|
Allink
v0.1
|
00001 #include "../include/VarData.h" 00002 00003 int VarData::InterParab(PART *PmIn,PART *PmOut,int NIn,int NOut){ 00004 SPLINE Par1; 00005 SPLINE Par2; 00006 double Maxx=0.; 00007 double Minx=0.; 00008 int NShowOut=1; 00009 for(int p=0;p<NIn;p++){ 00010 if(PmIn[p].Pos[0] < Minx) 00011 Minx = PmIn[p].Pos[0]; 00012 if(PmIn[p].Pos[0] > Maxx) 00013 Maxx = PmIn[p].Pos[0]; 00014 } 00015 double DeltaxIn=(Maxx-Minx)/(double)(NIn-1); 00016 double DeltaxOut=(Maxx-Minx)/(double)(NOut-1); 00017 double x = Minx; 00018 for(int p=0,pp=0;p<NIn;p++){ 00019 if(p < NIn-2){ 00020 Par1 = Mat->Parab(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,1); 00021 Par2 = Mat->Parab(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,2); 00022 } 00023 for(;x < DeltaxIn*(p+1);x += DeltaxOut){ 00024 if(pp == NOut-1) return NOut; 00025 PmOut[pp].Pos[0] = x; 00026 PmOut[pp].Pos[1] = Par1.a2*x*x+Par1.a1*x+Par1.a0; 00027 PmOut[pp].Pos[2] = Par2.a2*x*x+Par2.a1*x+Par2.a0; 00028 // printf("%d) (%lf %lf %lf) \n",p,PmOut[pp].Pos[0],PmOut[pp].Pos[1],PmOut[pp].Pos[2]); 00029 pp++; 00030 NShowOut = pp; 00031 PmOut[pp].Idx = PmIn[p].Idx; 00032 } 00033 } 00034 return NShowOut; 00035 } 00036 int VarData::InterParab2(PART *PmIn,PART *PmOut,int NIn,int NOut){ 00037 SPLINE Par1; 00038 SPLINE Par2; 00039 double Maxx=PmIn[0].Pos[0]; 00040 double Minx=PmIn[0].Pos[0]; 00041 for(int p=0;p<NIn;p++){ 00042 if(PmIn[p].Pos[0] < Minx) 00043 Minx = PmIn[p].Pos[0]; 00044 if(PmIn[p].Pos[0] > Maxx) 00045 Maxx = PmIn[p].Pos[0]; 00046 } 00047 double DeltaxIn=(double)NIn/((Maxx-Minx)*(double)NOut); 00048 double DeltaxOut=0.; 00049 int NOutShow=NOut; 00050 double x = Minx; 00051 for(int p=0,pp=0;p<NIn;p++){ 00052 if(p == 0){ 00053 Par1 = Mat->Parab2(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,1); 00054 Par2 = Mat->Parab2(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,2); 00055 DeltaxOut = (PmIn[p+1].Pos[0] - PmIn[p].Pos[0])*DeltaxIn; 00056 for(double x = PmIn[p].Pos[0], Stepx = 0.; 00057 x <= PmIn[p+1].Pos[0]; 00058 x += DeltaxOut,pp++){ 00059 Stepx = (x-PmIn[p].Pos[0]); 00060 PmOut[pp].Pos[0] = x; 00061 PmOut[pp].Pos[1] = Par1.a2*Stepx*Stepx+Par1.a1*Stepx+Par1.a0; 00062 PmOut[pp].Pos[2] = Par2.a2*Stepx*Stepx+Par2.a1*Stepx+Par2.a0; 00063 } 00064 continue; 00065 } 00066 else if(p < NIn-1){ 00067 Par1 = Mat->Parab2(PmIn[p-1].Pos,PmIn[p].Pos,PmIn[p+1].Pos,0,1); 00068 Par2 = Mat->Parab2(PmIn[p-1].Pos,PmIn[p].Pos,PmIn[p+1].Pos,0,2); 00069 DeltaxOut = (PmIn[p+1].Pos[0] - PmIn[p].Pos[0])*DeltaxIn; 00070 for(double x = PmIn[p].Pos[0], Stepx = 0.; 00071 x <= PmIn[p+1].Pos[0]; 00072 x += DeltaxOut,pp++){ 00073 Stepx = (x-PmIn[p].Pos[0]); 00074 if(pp == NOut-1) return NOut; 00075 PmOut[pp].Pos[0] = x; 00076 PmOut[pp].Pos[1] = Par1.a2*Stepx*Stepx+Par1.a1*Stepx+Par1.a0; 00077 PmOut[pp].Pos[2] = Par2.a2*Stepx*Stepx+Par2.a1*Stepx+Par2.a0; 00078 NOutShow = pp; 00079 PmOut[pp].Idx = PmIn[p].Idx; 00080 //printf("%d/%d) %lf (%lf %lf %lf) \n",p,pp,Maxx,PmOut[pp].Pos[0],PmOut[pp].Pos[1],PmOut[pp].Pos[2]); 00081 } 00082 } 00083 } 00084 return NOutShow; 00085 } 00086 int VarData::InterCubica(PART *PmIn,PART *PmOut,int NIn,int NOut){ 00087 SPLINE Cub1; 00088 SPLINE Cub2; 00089 double Maxx=PmIn[0].Pos[0]; 00090 double Minx=PmIn[0].Pos[0]; 00091 for(int p=0;p<NIn;p++){ 00092 if(PmIn[p].Pos[0] < Minx) 00093 Minx = PmIn[p].Pos[0]; 00094 if(PmIn[p].Pos[0] > Maxx) 00095 Maxx = PmIn[p].Pos[0]; 00096 } 00097 double DeltaxIn=(double)NIn/((Maxx-Minx)*(double)NOut); 00098 double DeltaxOut=0.; 00099 int NOutShow=NOut; 00100 double x = Minx; 00101 for(int p=0,pp=0;p<NIn-2;p++){ 00102 Cub1 = Mat->Cubica(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,PmIn[p+3].Pos,0,1); 00103 Cub2 = Mat->Cubica(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,PmIn[p+3].Pos,0,2); 00104 DeltaxOut = (PmIn[p+1].Pos[0] - PmIn[p].Pos[0])*DeltaxIn; 00105 for(double x = PmIn[p].Pos[0], Dx = 0.; 00106 x <= PmIn[p+1].Pos[0]; 00107 x += DeltaxOut,pp++){ 00108 Dx = (x-PmIn[p+1].Pos[0]); 00109 if(pp == NOut-1) return NOutShow; 00110 PmOut[pp].Pos[0] = x; 00111 PmOut[pp].Pos[1] = Cub1.a0+Cub1.a1*Dx+Cub1.a2*Dx*Dx+Cub1.a3*Dx*Dx*Dx; 00112 PmOut[pp].Pos[2] = Cub2.a0+Cub2.a1*Dx+Cub2.a2*Dx*Dx+Cub2.a3*Dx*Dx*Dx; 00113 NOutShow = pp; 00114 PmOut[pp].Idx = PmIn[p].Idx; 00115 //printf("%d/%d) %lf (%lf %lf %lf) \n",p,pp,Maxx,PmOut[pp].Pos[0],PmOut[pp].Pos[1],PmOut[pp].Pos[2]); 00116 } 00117 } 00118 return NOutShow; 00119 } 00120 int VarData::InterForth(PART *PmIn,PART *PmOut,int NIn,int NOut){ 00121 SPLINE Forth1; 00122 SPLINE Forth2; 00123 double Maxx=PmIn[0].Pos[0]; 00124 double Minx=PmIn[0].Pos[0]; 00125 for(int p=0;p<NIn;p++){ 00126 if(PmIn[p].Pos[0] < Minx) 00127 Minx = PmIn[p].Pos[0]; 00128 if(PmIn[p].Pos[0] > Maxx) 00129 Maxx = PmIn[p].Pos[0]; 00130 } 00131 double DeltaxIn=(double)NIn/((Maxx-Minx)*(double)NOut); 00132 double DeltaxOut=0.; 00133 int NOutShow=NOut; 00134 double x = Minx; 00135 double Ref; 00136 for(int p=0,pp=0;p<NIn;p++){ 00137 if(p < 1){ 00138 Forth1 = Mat->Forth(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,PmIn[p+3].Pos,PmIn[p+4].Pos,0,1); 00139 Forth2 = Mat->Forth(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,PmIn[p+3].Pos,PmIn[p+4].Pos,0,2); 00140 Ref = PmIn[p+2].Pos[0]; 00141 } 00142 else if(p < 2){ 00143 Forth1 = Mat->Forth(PmIn[p-1].Pos,PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,PmIn[p+3].Pos,0,1); 00144 Forth2 = Mat->Forth(PmIn[p-1].Pos,PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,PmIn[p+3].Pos,0,2); 00145 Ref = PmIn[p+1].Pos[0]; 00146 } 00147 else if(p < NIn-2){ 00148 Forth1 = Mat->Forth(PmIn[p-2].Pos,PmIn[p-1].Pos,PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,1); 00149 Forth2 = Mat->Forth(PmIn[p-2].Pos,PmIn[p-1].Pos,PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,2); 00150 Ref = PmIn[p].Pos[0]; 00151 } 00152 else if(p == NIn-2){ 00153 Forth1 = Mat->Forth(PmIn[p-3].Pos,PmIn[p-2].Pos,PmIn[p-1].Pos,PmIn[p].Pos,PmIn[p+1].Pos,0,1); 00154 Forth2 = Mat->Forth(PmIn[p-3].Pos,PmIn[p-2].Pos,PmIn[p-1].Pos,PmIn[p].Pos,PmIn[p+1].Pos,0,2); 00155 Ref = PmIn[p-1].Pos[0]; 00156 } 00157 else if(p == NIn-1){ 00158 Forth1 = Mat->Forth(PmIn[p-4].Pos,PmIn[p-3].Pos,PmIn[p-2].Pos,PmIn[p-1].Pos,PmIn[p].Pos,0,1); 00159 Forth2 = Mat->Forth(PmIn[p-4].Pos,PmIn[p-3].Pos,PmIn[p-2].Pos,PmIn[p-1].Pos,PmIn[p].Pos,0,2); 00160 Ref = PmIn[p-2].Pos[0]; 00161 } 00162 //printf("%d/%d (%lf %lf %lf) Sp (%lf %lf %lf %lf)\n",p,pp,PmIn[p].Pos[0],PmIn[p].Pos[1],PmIn[p].Pos[2],Forth2.a0,Forth2.a1,Forth2.a2,Forth2.a3); 00163 DeltaxOut = (PmIn[p+1].Pos[0] - PmIn[p].Pos[0])*DeltaxIn; 00164 for(double x = PmIn[p].Pos[0], Dx = 0.; 00165 x <= PmIn[p+1].Pos[0]; 00166 x += DeltaxOut,pp++){ 00167 Dx = (x-Ref); 00168 if(pp == NOut-1) return NOut; 00169 PmOut[pp].Pos[0] = x; 00170 PmOut[pp].Pos[1] = Forth1.a0+Forth1.a1*Dx+Forth1.a2*Dx*Dx+Forth1.a3*Dx*Dx*Dx+Forth1.a4*Dx*Dx*Dx*Dx; 00171 PmOut[pp].Pos[2] = Forth2.a0+Forth2.a1*Dx+Forth2.a2*Dx*Dx+Forth2.a3*Dx*Dx*Dx+Forth2.a4*Dx*Dx*Dx*Dx;; 00172 NOutShow = pp; 00173 PmOut[pp].Idx = PmIn[p].Idx; 00174 //printf("%d/%d) %lf (%lf %lf %lf) \n",p,pp,Maxx,PmOut[pp].Pos[0],PmOut[pp].Pos[1],PmOut[pp].Pos[2]); 00175 } 00176 } 00177 return NOutShow; 00178 }int VarData::InterSpline3(PART *PmIn,PART *PmOut,int NIn,int NOut){ 00179 SPLINE Sp1; 00180 SPLINE Sp2; 00181 double Maxx=PmIn[0].Pos[0]; 00182 double Minx=PmIn[0].Pos[0]; 00183 for(int p=0;p<NIn;p++){ 00184 if(PmIn[p].Pos[0] < Minx) 00185 Minx = PmIn[p].Pos[0]; 00186 if(PmIn[p].Pos[0] > Maxx) 00187 Maxx = PmIn[p].Pos[0]; 00188 } 00189 double DeltaxIn=(double)(NIn-1)/((Maxx-Minx)*(double)NOut); 00190 double DeltaxOut=0.; 00191 int NOutShow=0; 00192 double RatioOutIn = DeltaxOut/DeltaxIn; 00193 for(int p=0,pp=0;p<NIn;p+=1){ 00194 if( p == 0) 00195 Sp1 = Mat->Spline3Beg(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,2); 00196 else if( p < NIn-2 ) 00197 Sp1 = Mat->Spline3(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,2); 00198 else if( p == NIn - 2 ) 00199 Sp1 = Mat->Spline3End(PmIn[p].Pos,PmIn[p+1].Pos,0,2); 00200 DeltaxOut = (PmIn[p+1].Pos[0] - PmIn[p].Pos[0])*DeltaxIn; 00201 for(double x = PmIn[p].Pos[0], Stepx = 0.; 00202 x <= PmIn[p+1].Pos[0]; 00203 x += DeltaxOut,pp++){ 00204 //printf("%d %d %lf %lf\n",p,pp,x,DeltaxIn*(p+1)); 00205 if(pp == NOut-1) break; 00206 Stepx = (x-PmIn[p].Pos[0]); 00207 PmOut[pp].Pos[0] = x; 00208 PmOut[pp].Pos[2] = Sp1.a3*Stepx*Stepx*Stepx + 00209 Sp1.a2*Stepx*Stepx + 00210 Sp1.a1*Stepx + Sp1.a0; 00211 } 00212 } 00213 for(int p=0,pp=0;p<NIn;p+=1){ 00214 //printf("%d/%d (%lf %lf %lf) Sp (%lf %lf %lf %lf)\n",p,pp,PmIn[p].Pos[0],PmIn[p].Pos[1],PmIn[p].Pos[2],Sp2.a0,Sp2.a1,Sp2.a2,Sp2.a3); 00215 if( p == 0) 00216 Sp2 = Mat->Spline3Beg(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,1); 00217 else if( p < NIn-2 ) 00218 Sp2 = Mat->Spline3(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,1); 00219 else if( p == NIn - 2 ) 00220 Sp2 = Mat->Spline3End(PmIn[p].Pos,PmIn[p+1].Pos,0,1); 00221 DeltaxOut = (PmIn[p+1].Pos[0] - PmIn[p].Pos[0])*DeltaxIn; 00222 for(double x = PmIn[p].Pos[0], Stepx = 0.; 00223 x <= PmIn[p+1].Pos[0]; 00224 x += DeltaxOut,pp++){ 00225 if(pp == NOut-1) return NOutShow; 00226 Stepx = (x-PmIn[p].Pos[0]); 00227 PmOut[pp].Pos[1] = Sp2.a3*Stepx*Stepx*Stepx + 00228 Sp2.a2*Stepx*Stepx + 00229 Sp2.a1*Stepx + Sp2.a0; 00230 //printf("%d/%d) %lf (%lf %lf %lf)\n",p,pp,Stepx,PmOut[pp].Pos[0],PmOut[pp].Pos[1],PmOut[pp].Pos[2]); 00231 NOutShow = pp; 00232 PmOut[pp].Idx = PmIn[p].Idx; 00233 } 00234 } 00235 return NOutShow; 00236 } 00237 int VarData::InterSpline4(PART *PmIn,PART *PmOut,int NIn,int NOut){ 00238 SPLINE Sp1; 00239 SPLINE Sp2; 00240 double Maxx=PmIn[0].Pos[0]; 00241 double Minx=PmIn[0].Pos[0]; 00242 for(int p=0;p<NIn;p++){ 00243 if(PmIn[p].Pos[0] < Minx) 00244 Minx = PmIn[p].Pos[0]; 00245 if(PmIn[p].Pos[0] > Maxx) 00246 Maxx = PmIn[p].Pos[0]; 00247 } 00248 double DeltaxIn=(double)(NIn-1)/((Maxx-Minx)*(double)NOut); 00249 double DeltaxOut=0.; 00250 int NOutShow=0; 00251 double RatioOutIn = DeltaxOut/DeltaxIn; 00252 for(int p=0,pp=0;p<NIn;p+=1){ 00253 if( p == 0) 00254 Sp1 = Mat->Spline4Beg(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,PmIn[p+3].Pos,0,1); 00255 //else if( p < NIn - 3 ) 00256 //Sp1 = Mat->Spline4(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,PmIn[p+3].Pos,0,1); 00257 //else if( p == NIn - 3 ) 00258 //Sp1 = Mat->Spline3(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,1); 00259 //Sp1 = Mat->Spline4PreEnd(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,1); 00260 else if(p < NIn - 2) 00261 Sp1 = Mat->Spline4(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,1); 00262 else if( p == NIn - 2 ) 00263 Sp1 = Mat->Spline3End(PmIn[p].Pos,PmIn[p+1].Pos,0,1); 00264 //Sp1 = Mat->Spline4End(PmIn[p].Pos,PmIn[p+1].Pos,0,1); 00265 DeltaxOut = (PmIn[p+1].Pos[0] - PmIn[p].Pos[0])*DeltaxIn; 00266 for(double x = PmIn[p].Pos[0], Dx = 0.; 00267 x <= PmIn[p+1].Pos[0]; 00268 x += DeltaxOut,pp++){ 00269 //printf("%d %d %lf %lf\n",p,pp,x,DeltaxIn*(p+1)); 00270 if(pp == NOut-1) break; 00271 Dx = (x-PmIn[p].Pos[0]); 00272 PmOut[pp].Pos[0] = x; 00273 PmOut[pp].Pos[1] = Sp1.a0+Sp1.a1*Dx+Sp1.a2*Dx*Dx+Sp1.a3*Dx*Dx*Dx+Sp1.a4*Dx*Dx*Dx*Dx; 00274 } 00275 } 00276 for(int p=0,pp=0;p<NIn;p+=1){ 00277 //printf("%d/%d (%lf %lf %lf) Sp (%lf %lf %lf %lf)\n",p,pp,PmIn[p].Pos[0],PmIn[p].Pos[1],PmIn[p].Pos[2],Sp2.a0,Sp2.a1,Sp2.a2,Sp2.a3); 00278 if( p == 0) 00279 Sp2 = Mat->Spline4Beg(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,PmIn[p+3].Pos,0,2); 00280 //else if( p < NIn - 3 ) 00281 //Sp2 = Mat->Spline4(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,PmIn[p+3].Pos,0,2); 00282 //else if( p == NIn - 3 ) 00283 //Sp2 = Mat->Spline3(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,2); 00284 //Sp2 = Mat->Spline4PreEnd(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,2); 00285 else if(p < NIn - 2){ 00286 Sp2 = Mat->Spline4(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,2); 00287 } 00288 else if( p == NIn - 2 ) 00289 //Sp2 = Mat->Spline4End(PmIn[p].Pos,PmIn[p+1].Pos,0,2); 00290 Sp2 = Mat->Spline3End(PmIn[p].Pos,PmIn[p+1].Pos,0,2); 00291 DeltaxOut = (PmIn[p+1].Pos[0] - PmIn[p].Pos[0])*DeltaxIn; 00292 for(double x = PmIn[p].Pos[0], Dx = 0.; 00293 x <= PmIn[p+1].Pos[0]; 00294 x += DeltaxOut,pp++){ 00295 if(pp == NOut-1) return NOutShow; 00296 Dx = (x-PmIn[p].Pos[0]); 00297 PmOut[pp].Pos[2] = Sp2.a0+Sp2.a1*Dx+Sp2.a2*Dx*Dx+Sp2.a3*Dx*Dx*Dx+Sp2.a4*Dx*Dx*Dx*Dx; 00298 //printf("%d/%d) %lf (%lf %lf %lf)\n",p,pp,Dx,PmOut[pp].Pos[0],PmOut[pp].Pos[1],PmOut[pp].Pos[2]); 00299 NOutShow = pp; 00300 PmOut[pp].Idx = PmIn[p].Idx; 00301 } 00302 } 00303 return NOutShow; 00304 } 00305 int VarData::InterBSpline(PART *PmIn,PART *PmOut,int NIn,int NOut){ 00306 double Maxx=PmIn[0].Pos[CLat1]; 00307 double Minx=PmIn[0].Pos[CLat1]; 00308 int NOutShow=0; 00309 for(int p=0;p<NIn;p++){ 00310 if(PmIn[p].Pos[CLat1] < Minx) 00311 Minx = PmIn[p].Pos[CLat1]; 00312 if(PmIn[p].Pos[CLat1] > Maxx) 00313 Maxx = PmIn[p].Pos[CLat1]; 00314 } 00315 double DeltaxIn=(Maxx-Minx)/(double)(NIn-1); 00316 double DeltaxOut=(Maxx-Minx)/(double)(NOut-1); 00317 int NOrder = 3+1; 00318 double *dArray = (double *)calloc(NIn+NOrder+1,sizeof(double)); 00319 for(int p=0;p<=NIn+NOrder;p++){ 00320 if(p<NOrder) 00321 dArray[p] = Minx; 00322 else if( (NOrder<=p) && (p<=NIn) ) 00323 dArray[p] = (p-NOrder+1)*DeltaxIn+Minx; 00324 else if( p>NIn) 00325 dArray[p] = (p-NOrder+2)*DeltaxIn+Minx; 00326 } 00327 for(int pp=0;pp<NOut-1;pp++){ 00328 for(int d=0;d<3;d++) PmOut[pp].Pos[d] = 0.; 00329 double x = DeltaxOut*pp+Minx; 00330 int vPos = (int)((x-Minx)/(Maxx-Minx)*NIn); 00331 for(int p=vPos-1;p<vPos+NOrder;p++){ 00332 if(p >= NIn || p <0) continue; 00333 double Blend = Mat->Blend(dArray,x,p,NOrder); 00334 PmOut[pp].Pos[0] += Blend * PmIn[p].Pos[0]; 00335 PmOut[pp].Pos[1] += Blend * PmIn[p].Pos[1]; 00336 PmOut[pp].Pos[2] += Blend * PmIn[p].Pos[2]; 00337 } 00338 } 00339 PmOut[NOut-1].Pos[0] = PmIn[NIn-1].Pos[0]; 00340 PmOut[NOut-1].Pos[1] = PmIn[NIn-1].Pos[1]; 00341 PmOut[NOut-1].Pos[2] = PmIn[NIn-1].Pos[2]; 00342 NOutShow = NOut; 00343 return NOutShow; 00344 } 00345 int VarData::InterBSpline2D(double **PlIn,double **PlOut,int NIn,int NOut){ 00346 double Max=Gen->Edge[CLat1]; 00347 double Min=0.; 00348 int NOutShow=0; 00349 double DeltaIn=(Max-Min)/(double)(NIn-1); 00350 double DeltaOut=(Max-Min)/(double)(NOut-1); 00351 double RatioInOut = (double)(NIn/(double)NOut); 00352 int NOrder = 3+1; 00353 double *dArray = (double *)calloc(NIn+NOrder+1,sizeof(double)); 00354 for(int p=0;p<=NIn+NOrder;p++){ 00355 if(p<NOrder){ 00356 dArray[p] = Min; 00357 } 00358 else if( (NOrder<=p) && (p<=NIn) ){ 00359 dArray[p] = (p-NOrder+1)*DeltaIn+Min; 00360 } 00361 else if(p>NIn){ 00362 dArray[p] = (p-NOrder+2)*DeltaIn+Min; 00363 } 00364 } 00365 for(int vo=0;vo<NOut;vo++){ 00366 for(int vvo=0;vvo<NOut;vvo++){ 00367 PlOut[vo][vvo] = 0.; 00368 double x = DeltaOut*vo+Min; 00369 int vxPos = (int)(vo*RatioInOut); 00370 for(int vi=vxPos-1,vn=0;vi<vxPos+NOrder+1;vi++){ 00371 vn = vi; 00372 if(vi < 0) vn = NIn + vi; 00373 if(vi >= NIn) vn = vi - NIn; 00374 double Blendx = Mat->Blend(dArray,x,vn,NOrder); 00375 double y = DeltaOut*vvo+Min; 00376 int vyPos = (int)(vvo*RatioInOut); 00377 for(int vvi=vyPos-1,vvn=0;vvi<vyPos+NOrder+1;vvi++){ 00378 vvn = vvi; 00379 if(vvi < 0) vvn = NIn + vvi; 00380 if(vvi >= NIn) vvn = vvi - NIn; 00381 double Blendy = Mat->Blend(dArray,y,vvn,NOrder); 00382 PlOut[vo][vvo] += Blendx*Blendy*PlIn[vn][vvn]; 00383 } 00384 } 00385 } 00386 } 00387 // for(int vo=0;vo<NOut;vo++){ 00388 // if(vo < NIn){//To be arranged 00389 // PlOut[vo][NOut-1] = PlIn[vo][NIn-1]; 00390 // PlOut[NOut - 1][vo] = PlIn[NIn-1][vo]; 00391 // } 00392 // } 00393 //for(int vo=0;vo<NOut;vo++)for(int vvo=0;vvo<NOut;vvo++)printf("%d %d %lf %lf \n",vo,vvo,PlOut[vo][vvo],PlIn[vo][vvo]); 00394 NOutShow = NOut; 00395 free(dArray); 00396 return NOutShow; 00397 } 00401 void VarData::SmoothGrid(int NSample){ 00402 double *PlotIn = (double *)calloc(SQR(NSample),sizeof(double)); 00403 double *PlotOut = (double *)calloc(SQR(NSample),sizeof(double)); 00404 for(int nx=0;nx<NSample;nx++){ 00405 for(int ny=0;ny<NSample;ny++){ 00406 int n = nx*NSample+ny; 00407 PlotIn[n] = Pm[n].Pos[2]; 00408 } 00409 } 00410 InterBSpline2D(PlotIn,PlotOut,NSample,NSample); 00411 for(int nx=0;nx<NSample;nx++){ 00412 for(int ny=0;ny<NSample;ny++){ 00413 int n = nx*NSample+ny; 00414 if(pType(n) != 0) continue; 00415 Pm[n].Pos[2] = PlotOut[n]; 00416 } 00417 } 00418 free(PlotIn); 00419 free(PlotOut); 00420 } 00424 int VarData::InterBSpline2D(double *PlIn,double *PlOut,int NIn,int NOut){ 00425 double Max=Gen->Edge[CLat1]; 00426 double Min=0.; 00427 int NOutShow=0; 00428 double DeltaIn = (Max-Min)/(double)(NIn-1); 00429 double DeltaOut = (Max-Min)/(double)(NOut-1); 00430 double RatioInOut = (double)(NIn/(double)NOut); 00431 int NOrder = 3+1; 00432 double *dArray = (double *)calloc(NIn+NOrder+1,sizeof(double)); 00433 for(int p=0;p<=NIn+NOrder;p++){ 00434 if(p<NOrder){ 00435 dArray[p] = Min; 00436 } 00437 else if( (NOrder<=p) && (p<=NIn) ){ 00438 dArray[p] = (p-NOrder+1)*DeltaIn+Min; 00439 } 00440 else if(p>NIn){ 00441 dArray[p] = (p-NOrder+2)*DeltaIn+Min; 00442 } 00443 } 00444 for(int vo=0;vo<NOut;vo++){ 00445 for(int vvo=0;vvo<NOut;vvo++){ 00446 PlOut[vo*NOut+vvo] = 0.; 00447 double x = DeltaOut*vo+Min; 00448 int vxPos = (int)(vo*RatioInOut); 00449 for(int vi=vxPos-1,vn=0;vi<vxPos+NOrder+1;vi++){ 00450 vn = vi; 00451 if(vi < 0) vn = NIn + vi; 00452 if(vi >= NIn) vn = vi - NIn; 00453 double Blendx = Mat->Blend(dArray,x,vn,NOrder); 00454 double y = DeltaOut*vvo+Min; 00455 int vyPos = (int)(vvo*RatioInOut); 00456 for(int vvi=vyPos-1,vvn=0;vvi<vyPos+NOrder+1;vvi++){ 00457 vvn = vvi; 00458 if(vvi < 0) vvn = NIn + vvi; 00459 if(vvi >= NIn) vvn = vvi - NIn; 00460 double Blendy = Mat->Blend(dArray,y,vvn,NOrder); 00461 PlOut[vo*NOut+vvo] += Blendx*Blendy*PlIn[vn*NIn+vvn]; 00462 } 00463 } 00464 } 00465 } 00466 // for(int vo=0;vo<NOut;vo++){ 00467 // if(vo < NIn){//To be arranged 00468 // PlOut[vo][NOut-1] = PlIn[vo][NIn-1]; 00469 // PlOut[NOut - 1][vo] = PlIn[NIn-1][vo]; 00470 // } 00471 // } 00472 //for(int vo=0;vo<NOut;vo++)for(int vvo=0;vvo<NOut;vvo++)printf("%d %d %lf %lf \n",vo,vvo,PlOut[vo][vvo],PlIn[vo][vvo]); 00473 NOutShow = NOut; 00474 free(dArray); 00475 return NOutShow; 00476 } 00477 //minimum image convention does not seem to work 00478 int VarData::InterBSpline1D(double *PlIn,double *PlOut,int NIn,int NOut){ 00479 double Max=Gen->Edge[CLat1]; 00480 double Min=0.; 00481 int NOutShow=0; 00482 double DeltaIn=(Max-Min)/(double)(NIn-1); 00483 double DeltaOut=(Max-Min)/(double)(NOut-1); 00484 int NOrder = 3+1; 00485 double *dArray = (double *)calloc(NIn+2*NOrder+2,sizeof(double)); 00486 for(int vi=0;vi<=NIn+2*NOrder;vi++){ 00487 if(vi<NOrder){ 00488 dArray[vi] = Min; 00489 } 00490 else if(vi > NIn){ 00491 dArray[vi] = (vi-NOrder+2)*DeltaIn+Min; 00492 } 00493 else { 00494 dArray[vi] = (vi-NOrder+1)*DeltaIn+Min; 00495 } 00496 } 00497 for(int vo=0;vo<NOut;vo++){ 00498 //to avoid the minimum image convention 00499 if(vo > NOut-NOrder) continue; 00500 PlOut[vo] = 0.; 00501 double x = DeltaOut*vo+Min; 00502 for(int vi=vo-1,vn=0;vi<vo+NOrder+1;vi++){ 00503 vn = vi; 00504 if(vi < 0){ 00505 vn = NIn + vi; 00506 } 00507 if(vi >= NIn){ 00508 vn = vi - NIn; 00509 } 00510 double Blendx = Mat->Blend(dArray,x,vn,NOrder); 00511 PlOut[vo] += Blendx*PlIn[vn]; 00512 } 00513 } 00514 NOutShow = NOut; 00515 free(dArray); 00516 return NOutShow; 00517 } 00518 int VarData::InterPoly(PART *PmIn,PART *PmOut,int NIn,int NOut){ 00519 double *Pz = (double *)calloc(NIn,sizeof(double)); 00520 double *Px = (double *)calloc(NIn,sizeof(double)); 00521 double *Py = (double *)calloc(NIn,sizeof(double)); 00522 for(int p=0;p<NIn;p++){ 00523 Px[p] = PmIn[p].Pos[CLat1]; 00524 Py[p] = PmIn[p].Pos[CLat2]; 00525 Pz[p] = PmIn[p].Pos[CNorm]; 00526 } 00527 Spline *Sp1 = new Spline(NIn); 00528 Spline *Sp2 = new Spline(NIn); 00529 Mat->Polinomio(Px,Py,NIn,Sp1); 00530 Mat->Polinomio(Px,Pz,NIn,Sp2); 00531 double Length = PmIn[NIn-1].Pos[CLat1] - PmIn[0].Pos[CLat1]; 00532 double Deltax = Length / (double)NOut; 00533 double Pos = PmIn[0].Pos[CLat1]; 00534 for(int p=0;p<NOut;p++){ 00535 Pos += Deltax; 00536 PmOut[p].Pos[CLat1] = Pos; 00537 PmOut[p].Pos[CLat2] = 0.; 00538 PmOut[p].Pos[CNorm] = 0.; 00539 for(int n=0;n<NIn;n++){ 00540 PmOut[p].Pos[CLat2] += Sp1->GetCoe(n)*Mat->Elevato(Pos,n); 00541 PmOut[p].Pos[CNorm] += Sp2->GetCoe(n)*Mat->Elevato(Pos,n); 00542 } 00543 //printf("(%lf %lf %lf)\n",PmOut[p].Pos[CLat1],PmOut[p].Pos[CLat2],PmOut[p].Pos[CNorm]); 00544 } 00545 delete Sp1;delete Sp2; 00546 free(Pz);free(Py);free(Px); 00547 return NOut; 00548 } 00549 // int VarData::InterDerMatrix(PART *Pm,int NMass,SPLINE Wg,double Offset){ 00550 // Matrice *Coeff = new Matrice(Wg,1); 00551 // Matrice *Resp = new Matrice(NMass,NMass); 00552 // for(int r=0;r<NMass;r++){ 00553 // if(Pm[r].Typ != 0){ Resp->Set(r,r,1.);continue;} 00554 // if(r >= 2) Resp->Set(r,r-2,Coeff->Val(2,0)); 00555 // if(r >= 1) Resp->Set(r,r-1,Coeff->Val(2,1)); 00556 // if(r < NMass-1) Resp->Set(r,r+1,Coeff->Val(2,3)); 00557 // if(r < NMass-2) Resp->Set(r,r+2,Coeff->Val(2,4)); 00558 // Resp->Set(r,r,Coeff->Val(2,2)); 00559 // } 00560 // double *Known = (double *) calloc(NMass,sizeof(double)); 00561 // double *UnKnown = (double *) calloc(NMass,sizeof(double)); 00562 // for(int p=0;p<NMass;p++) 00563 // Known[p] = Pm[p].Pos[CNorm] - Offset; 00564 // Resp->Solve(Known,UnKnown); 00565 // for(int p=0;p<NMass;p++){ 00566 // if(Pm[p].Typ != 0){continue;} 00567 // Pm[p].Pos[CNorm] = UnKnown[p] + Offset; 00568 // //printf("%lf %lf\n",Known[p],UnKnown[p]); 00569 // } 00570 // free(Known); 00571 // free(UnKnown); 00572 // // Resp->invert(); 00573 // // for(int r=0;r<NMass;r++){ 00574 // // double Agg=0.; 00575 // // for(int c=0;c<NMass;c++){ 00576 // // if(Pm[r].Typ == 0) 00577 // // Agg += Resp->Val(r,c) * Pm[c].Pos[CNorm]; 00578 // // else 00579 // // Agg += Resp->Val(r,c) * Pm[c].Pos[CNorm]; 00580 // // } 00581 // // Pm[r].Pos[CNorm] = Agg; 00582 // // } 00583 // delete Coeff; 00584 // delete Resp; 00585 // return NMass; 00586 // } 00587 void VarData::SmoothGrid(int NSample,char *FWrite){ 00588 double InvNSample = 1./(double)NSample; 00589 double **PlotIn = (double **)calloc(SQR(NSample),sizeof(double)); 00590 double **PlotOut = (double **)calloc(SQR(NSample),sizeof(double)); 00591 double **Count = (double **)calloc(SQR(NSample),sizeof(double)); 00592 double Round = 0.00001; 00593 for(int v=0;v<NSample;v++){ 00594 PlotIn[v] = (double *)calloc(NSample,sizeof(double)); 00595 PlotOut[v] = (double *)calloc(NSample,sizeof(double)); 00596 Count[v] = (double *)calloc(NSample,sizeof(double)); 00597 } 00598 int Nx = 0; 00599 int Ny = 0; 00600 for(int p=0;p<pNPart();p++){ 00601 if(Pm[p].Pos[CLat2] < Pm[p+1].Pos[CLat2]) Ny++; 00602 else break; 00603 } 00604 Ny++; 00605 Nx = (int)(pNPart()/(double)Ny); 00606 for(int p=0;p<pNPart();p++){ 00607 int vx = (int)((Pm[p].Pos[CLat1]+Round)*pInvEdge(CLat1)*NSample); 00608 int vy = (int)((Pm[p].Pos[CLat2]+Round)*pInvEdge(CLat2)*NSample); 00609 if(vx < 0 || vx >= NSample)continue; 00610 if(vy < 0 || vy >= NSample)continue; 00611 PlotIn[vx][vy] += Pm[p].Pos[CNorm]; 00612 Count[vx][vy] += 1.; 00613 } 00614 for(int vx=0;vx<NSample;vx++) 00615 for(int vy=0;vy<NSample;vy++) 00616 PlotIn[vx][vy] *= Count[vx][vy] > 0. ? 1./Count[vx][vy] : 1.; 00617 InterBSpline2D(PlotIn,PlotOut,NSample,NSample); 00618 if(1==0){ 00619 FILE *F2Write = fopen(FWrite,"w"); 00620 for(int vx=0;vx<NSample;vx++){ 00621 for(int vy=0;vy<NSample;vy++){ 00622 double x = vx*pEdge(CLat1)*InvNSample; 00623 double y = vy*pEdge(CLat2)*InvNSample; 00624 fprintf(F2Write,"%lf %lf %lf\n",x,y,PlotOut[vx][vy]); 00625 } 00626 } 00627 fclose(F2Write); 00628 } 00629 for(int v=0;v<NSample;v++){ 00630 free(PlotIn[v]); 00631 free(PlotOut[v]); 00632 free(Count[v]); 00633 } 00634 free(PlotIn); 00635 free(PlotOut); 00636 free(Count); 00637 }
1.7.6.1