Allink  v0.1
VarDataInterp.cpp
1 #include "../include/VarData.h"
2 
3 int VarData::InterParab(PART *PmIn,PART *PmOut,int NIn,int NOut){
4  SPLINE Par1;
5  SPLINE Par2;
6  double Maxx=0.;
7  double Minx=0.;
8  int NShowOut=1;
9  for(int p=0;p<NIn;p++){
10  if(PmIn[p].Pos[0] < Minx)
11  Minx = PmIn[p].Pos[0];
12  if(PmIn[p].Pos[0] > Maxx)
13  Maxx = PmIn[p].Pos[0];
14  }
15  double DeltaxIn=(Maxx-Minx)/(double)(NIn-1);
16  double DeltaxOut=(Maxx-Minx)/(double)(NOut-1);
17  double x = Minx;
18  for(int p=0,pp=0;p<NIn;p++){
19  if(p < NIn-2){
20  Par1 = Mat->Parab(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,1);
21  Par2 = Mat->Parab(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,2);
22  }
23  for(;x < DeltaxIn*(p+1);x += DeltaxOut){
24  if(pp == NOut-1) return NOut;
25  PmOut[pp].Pos[0] = x;
26  PmOut[pp].Pos[1] = Par1.a2*x*x+Par1.a1*x+Par1.a0;
27  PmOut[pp].Pos[2] = Par2.a2*x*x+Par2.a1*x+Par2.a0;
28  // printf("%d) (%lf %lf %lf) \n",p,PmOut[pp].Pos[0],PmOut[pp].Pos[1],PmOut[pp].Pos[2]);
29  pp++;
30  NShowOut = pp;
31  PmOut[pp].Idx = PmIn[p].Idx;
32  }
33  }
34  return NShowOut;
35 }
36 int VarData::InterParab2(PART *PmIn,PART *PmOut,int NIn,int NOut){
37  SPLINE Par1;
38  SPLINE Par2;
39  double Maxx=PmIn[0].Pos[0];
40  double Minx=PmIn[0].Pos[0];
41  for(int p=0;p<NIn;p++){
42  if(PmIn[p].Pos[0] < Minx)
43  Minx = PmIn[p].Pos[0];
44  if(PmIn[p].Pos[0] > Maxx)
45  Maxx = PmIn[p].Pos[0];
46  }
47  double DeltaxIn=(double)NIn/((Maxx-Minx)*(double)NOut);
48  double DeltaxOut=0.;
49  int NOutShow=NOut;
50  double x = Minx;
51  for(int p=0,pp=0;p<NIn;p++){
52  if(p == 0){
53  Par1 = Mat->Parab2(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,1);
54  Par2 = Mat->Parab2(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,2);
55  DeltaxOut = (PmIn[p+1].Pos[0] - PmIn[p].Pos[0])*DeltaxIn;
56  for(double x = PmIn[p].Pos[0], Stepx = 0.;
57  x <= PmIn[p+1].Pos[0];
58  x += DeltaxOut,pp++){
59  Stepx = (x-PmIn[p].Pos[0]);
60  PmOut[pp].Pos[0] = x;
61  PmOut[pp].Pos[1] = Par1.a2*Stepx*Stepx+Par1.a1*Stepx+Par1.a0;
62  PmOut[pp].Pos[2] = Par2.a2*Stepx*Stepx+Par2.a1*Stepx+Par2.a0;
63  }
64  continue;
65  }
66  else if(p < NIn-1){
67  Par1 = Mat->Parab2(PmIn[p-1].Pos,PmIn[p].Pos,PmIn[p+1].Pos,0,1);
68  Par2 = Mat->Parab2(PmIn[p-1].Pos,PmIn[p].Pos,PmIn[p+1].Pos,0,2);
69  DeltaxOut = (PmIn[p+1].Pos[0] - PmIn[p].Pos[0])*DeltaxIn;
70  for(double x = PmIn[p].Pos[0], Stepx = 0.;
71  x <= PmIn[p+1].Pos[0];
72  x += DeltaxOut,pp++){
73  Stepx = (x-PmIn[p].Pos[0]);
74  if(pp == NOut-1) return NOut;
75  PmOut[pp].Pos[0] = x;
76  PmOut[pp].Pos[1] = Par1.a2*Stepx*Stepx+Par1.a1*Stepx+Par1.a0;
77  PmOut[pp].Pos[2] = Par2.a2*Stepx*Stepx+Par2.a1*Stepx+Par2.a0;
78  NOutShow = pp;
79  PmOut[pp].Idx = PmIn[p].Idx;
80  //printf("%d/%d) %lf (%lf %lf %lf) \n",p,pp,Maxx,PmOut[pp].Pos[0],PmOut[pp].Pos[1],PmOut[pp].Pos[2]);
81  }
82  }
83  }
84  return NOutShow;
85 }
86 int VarData::InterCubica(PART *PmIn,PART *PmOut,int NIn,int NOut){
87  SPLINE Cub1;
88  SPLINE Cub2;
89  double Maxx=PmIn[0].Pos[0];
90  double Minx=PmIn[0].Pos[0];
91  for(int p=0;p<NIn;p++){
92  if(PmIn[p].Pos[0] < Minx)
93  Minx = PmIn[p].Pos[0];
94  if(PmIn[p].Pos[0] > Maxx)
95  Maxx = PmIn[p].Pos[0];
96  }
97  double DeltaxIn=(double)NIn/((Maxx-Minx)*(double)NOut);
98  double DeltaxOut=0.;
99  int NOutShow=NOut;
100  double x = Minx;
101  for(int p=0,pp=0;p<NIn-2;p++){
102  Cub1 = Mat->Cubica(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,PmIn[p+3].Pos,0,1);
103  Cub2 = Mat->Cubica(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,PmIn[p+3].Pos,0,2);
104  DeltaxOut = (PmIn[p+1].Pos[0] - PmIn[p].Pos[0])*DeltaxIn;
105  for(double x = PmIn[p].Pos[0], Dx = 0.;
106  x <= PmIn[p+1].Pos[0];
107  x += DeltaxOut,pp++){
108  Dx = (x-PmIn[p+1].Pos[0]);
109  if(pp == NOut-1) return NOutShow;
110  PmOut[pp].Pos[0] = x;
111  PmOut[pp].Pos[1] = Cub1.a0+Cub1.a1*Dx+Cub1.a2*Dx*Dx+Cub1.a3*Dx*Dx*Dx;
112  PmOut[pp].Pos[2] = Cub2.a0+Cub2.a1*Dx+Cub2.a2*Dx*Dx+Cub2.a3*Dx*Dx*Dx;
113  NOutShow = pp;
114  PmOut[pp].Idx = PmIn[p].Idx;
115  //printf("%d/%d) %lf (%lf %lf %lf) \n",p,pp,Maxx,PmOut[pp].Pos[0],PmOut[pp].Pos[1],PmOut[pp].Pos[2]);
116  }
117  }
118  return NOutShow;
119 }
120 int VarData::InterForth(PART *PmIn,PART *PmOut,int NIn,int NOut){
121  SPLINE Forth1;
122  SPLINE Forth2;
123  double Maxx=PmIn[0].Pos[0];
124  double Minx=PmIn[0].Pos[0];
125  for(int p=0;p<NIn;p++){
126  if(PmIn[p].Pos[0] < Minx)
127  Minx = PmIn[p].Pos[0];
128  if(PmIn[p].Pos[0] > Maxx)
129  Maxx = PmIn[p].Pos[0];
130  }
131  double DeltaxIn=(double)NIn/((Maxx-Minx)*(double)NOut);
132  double DeltaxOut=0.;
133  int NOutShow=NOut;
134  double x = Minx;
135  double Ref;
136  for(int p=0,pp=0;p<NIn;p++){
137  if(p < 1){
138  Forth1 = Mat->Forth(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,PmIn[p+3].Pos,PmIn[p+4].Pos,0,1);
139  Forth2 = Mat->Forth(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,PmIn[p+3].Pos,PmIn[p+4].Pos,0,2);
140  Ref = PmIn[p+2].Pos[0];
141  }
142  else if(p < 2){
143  Forth1 = Mat->Forth(PmIn[p-1].Pos,PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,PmIn[p+3].Pos,0,1);
144  Forth2 = Mat->Forth(PmIn[p-1].Pos,PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,PmIn[p+3].Pos,0,2);
145  Ref = PmIn[p+1].Pos[0];
146  }
147  else if(p < NIn-2){
148  Forth1 = Mat->Forth(PmIn[p-2].Pos,PmIn[p-1].Pos,PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,1);
149  Forth2 = Mat->Forth(PmIn[p-2].Pos,PmIn[p-1].Pos,PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,2);
150  Ref = PmIn[p].Pos[0];
151  }
152  else if(p == NIn-2){
153  Forth1 = Mat->Forth(PmIn[p-3].Pos,PmIn[p-2].Pos,PmIn[p-1].Pos,PmIn[p].Pos,PmIn[p+1].Pos,0,1);
154  Forth2 = Mat->Forth(PmIn[p-3].Pos,PmIn[p-2].Pos,PmIn[p-1].Pos,PmIn[p].Pos,PmIn[p+1].Pos,0,2);
155  Ref = PmIn[p-1].Pos[0];
156  }
157  else if(p == NIn-1){
158  Forth1 = Mat->Forth(PmIn[p-4].Pos,PmIn[p-3].Pos,PmIn[p-2].Pos,PmIn[p-1].Pos,PmIn[p].Pos,0,1);
159  Forth2 = Mat->Forth(PmIn[p-4].Pos,PmIn[p-3].Pos,PmIn[p-2].Pos,PmIn[p-1].Pos,PmIn[p].Pos,0,2);
160  Ref = PmIn[p-2].Pos[0];
161  }
162  //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);
163  DeltaxOut = (PmIn[p+1].Pos[0] - PmIn[p].Pos[0])*DeltaxIn;
164  for(double x = PmIn[p].Pos[0], Dx = 0.;
165  x <= PmIn[p+1].Pos[0];
166  x += DeltaxOut,pp++){
167  Dx = (x-Ref);
168  if(pp == NOut-1) return NOut;
169  PmOut[pp].Pos[0] = x;
170  PmOut[pp].Pos[1] = Forth1.a0+Forth1.a1*Dx+Forth1.a2*Dx*Dx+Forth1.a3*Dx*Dx*Dx+Forth1.a4*Dx*Dx*Dx*Dx;
171  PmOut[pp].Pos[2] = Forth2.a0+Forth2.a1*Dx+Forth2.a2*Dx*Dx+Forth2.a3*Dx*Dx*Dx+Forth2.a4*Dx*Dx*Dx*Dx;;
172  NOutShow = pp;
173  PmOut[pp].Idx = PmIn[p].Idx;
174  //printf("%d/%d) %lf (%lf %lf %lf) \n",p,pp,Maxx,PmOut[pp].Pos[0],PmOut[pp].Pos[1],PmOut[pp].Pos[2]);
175  }
176  }
177  return NOutShow;
178 }int VarData::InterSpline3(PART *PmIn,PART *PmOut,int NIn,int NOut){
179  SPLINE Sp1;
180  SPLINE Sp2;
181  double Maxx=PmIn[0].Pos[0];
182  double Minx=PmIn[0].Pos[0];
183  for(int p=0;p<NIn;p++){
184  if(PmIn[p].Pos[0] < Minx)
185  Minx = PmIn[p].Pos[0];
186  if(PmIn[p].Pos[0] > Maxx)
187  Maxx = PmIn[p].Pos[0];
188  }
189  double DeltaxIn=(double)(NIn-1)/((Maxx-Minx)*(double)NOut);
190  double DeltaxOut=0.;
191  int NOutShow=0;
192  double RatioOutIn = DeltaxOut/DeltaxIn;
193  for(int p=0,pp=0;p<NIn;p+=1){
194  if( p == 0)
195  Sp1 = Mat->Spline3Beg(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,2);
196  else if( p < NIn-2 )
197  Sp1 = Mat->Spline3(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,2);
198  else if( p == NIn - 2 )
199  Sp1 = Mat->Spline3End(PmIn[p].Pos,PmIn[p+1].Pos,0,2);
200  DeltaxOut = (PmIn[p+1].Pos[0] - PmIn[p].Pos[0])*DeltaxIn;
201  for(double x = PmIn[p].Pos[0], Stepx = 0.;
202  x <= PmIn[p+1].Pos[0];
203  x += DeltaxOut,pp++){
204  //printf("%d %d %lf %lf\n",p,pp,x,DeltaxIn*(p+1));
205  if(pp == NOut-1) break;
206  Stepx = (x-PmIn[p].Pos[0]);
207  PmOut[pp].Pos[0] = x;
208  PmOut[pp].Pos[2] = Sp1.a3*Stepx*Stepx*Stepx +
209  Sp1.a2*Stepx*Stepx +
210  Sp1.a1*Stepx + Sp1.a0;
211  }
212  }
213  for(int p=0,pp=0;p<NIn;p+=1){
214  //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);
215  if( p == 0)
216  Sp2 = Mat->Spline3Beg(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,1);
217  else if( p < NIn-2 )
218  Sp2 = Mat->Spline3(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,1);
219  else if( p == NIn - 2 )
220  Sp2 = Mat->Spline3End(PmIn[p].Pos,PmIn[p+1].Pos,0,1);
221  DeltaxOut = (PmIn[p+1].Pos[0] - PmIn[p].Pos[0])*DeltaxIn;
222  for(double x = PmIn[p].Pos[0], Stepx = 0.;
223  x <= PmIn[p+1].Pos[0];
224  x += DeltaxOut,pp++){
225  if(pp == NOut-1) return NOutShow;
226  Stepx = (x-PmIn[p].Pos[0]);
227  PmOut[pp].Pos[1] = Sp2.a3*Stepx*Stepx*Stepx +
228  Sp2.a2*Stepx*Stepx +
229  Sp2.a1*Stepx + Sp2.a0;
230  //printf("%d/%d) %lf (%lf %lf %lf)\n",p,pp,Stepx,PmOut[pp].Pos[0],PmOut[pp].Pos[1],PmOut[pp].Pos[2]);
231  NOutShow = pp;
232  PmOut[pp].Idx = PmIn[p].Idx;
233  }
234  }
235  return NOutShow;
236 }
237 int VarData::InterSpline4(PART *PmIn,PART *PmOut,int NIn,int NOut){
238  SPLINE Sp1;
239  SPLINE Sp2;
240  double Maxx=PmIn[0].Pos[0];
241  double Minx=PmIn[0].Pos[0];
242  for(int p=0;p<NIn;p++){
243  if(PmIn[p].Pos[0] < Minx)
244  Minx = PmIn[p].Pos[0];
245  if(PmIn[p].Pos[0] > Maxx)
246  Maxx = PmIn[p].Pos[0];
247  }
248  double DeltaxIn=(double)(NIn-1)/((Maxx-Minx)*(double)NOut);
249  double DeltaxOut=0.;
250  int NOutShow=0;
251  double RatioOutIn = DeltaxOut/DeltaxIn;
252  for(int p=0,pp=0;p<NIn;p+=1){
253  if( p == 0)
254  Sp1 = Mat->Spline4Beg(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,PmIn[p+3].Pos,0,1);
255  //else if( p < NIn - 3 )
256  //Sp1 = Mat->Spline4(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,PmIn[p+3].Pos,0,1);
257  //else if( p == NIn - 3 )
258  //Sp1 = Mat->Spline3(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,1);
259  //Sp1 = Mat->Spline4PreEnd(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,1);
260  else if(p < NIn - 2)
261  Sp1 = Mat->Spline4(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,1);
262  else if( p == NIn - 2 )
263  Sp1 = Mat->Spline3End(PmIn[p].Pos,PmIn[p+1].Pos,0,1);
264  //Sp1 = Mat->Spline4End(PmIn[p].Pos,PmIn[p+1].Pos,0,1);
265  DeltaxOut = (PmIn[p+1].Pos[0] - PmIn[p].Pos[0])*DeltaxIn;
266  for(double x = PmIn[p].Pos[0], Dx = 0.;
267  x <= PmIn[p+1].Pos[0];
268  x += DeltaxOut,pp++){
269  //printf("%d %d %lf %lf\n",p,pp,x,DeltaxIn*(p+1));
270  if(pp == NOut-1) break;
271  Dx = (x-PmIn[p].Pos[0]);
272  PmOut[pp].Pos[0] = x;
273  PmOut[pp].Pos[1] = Sp1.a0+Sp1.a1*Dx+Sp1.a2*Dx*Dx+Sp1.a3*Dx*Dx*Dx+Sp1.a4*Dx*Dx*Dx*Dx;
274  }
275  }
276  for(int p=0,pp=0;p<NIn;p+=1){
277  //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);
278  if( p == 0)
279  Sp2 = Mat->Spline4Beg(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,PmIn[p+3].Pos,0,2);
280  //else if( p < NIn - 3 )
281  //Sp2 = Mat->Spline4(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,PmIn[p+3].Pos,0,2);
282  //else if( p == NIn - 3 )
283  //Sp2 = Mat->Spline3(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,2);
284  //Sp2 = Mat->Spline4PreEnd(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,2);
285  else if(p < NIn - 2){
286  Sp2 = Mat->Spline4(PmIn[p].Pos,PmIn[p+1].Pos,PmIn[p+2].Pos,0,2);
287  }
288  else if( p == NIn - 2 )
289  //Sp2 = Mat->Spline4End(PmIn[p].Pos,PmIn[p+1].Pos,0,2);
290  Sp2 = Mat->Spline3End(PmIn[p].Pos,PmIn[p+1].Pos,0,2);
291  DeltaxOut = (PmIn[p+1].Pos[0] - PmIn[p].Pos[0])*DeltaxIn;
292  for(double x = PmIn[p].Pos[0], Dx = 0.;
293  x <= PmIn[p+1].Pos[0];
294  x += DeltaxOut,pp++){
295  if(pp == NOut-1) return NOutShow;
296  Dx = (x-PmIn[p].Pos[0]);
297  PmOut[pp].Pos[2] = Sp2.a0+Sp2.a1*Dx+Sp2.a2*Dx*Dx+Sp2.a3*Dx*Dx*Dx+Sp2.a4*Dx*Dx*Dx*Dx;
298  //printf("%d/%d) %lf (%lf %lf %lf)\n",p,pp,Dx,PmOut[pp].Pos[0],PmOut[pp].Pos[1],PmOut[pp].Pos[2]);
299  NOutShow = pp;
300  PmOut[pp].Idx = PmIn[p].Idx;
301  }
302  }
303  return NOutShow;
304 }
305 int VarData::InterBSpline(PART *PmIn,PART *PmOut,int NIn,int NOut){
306  double Maxx=PmIn[0].Pos[CLat1];
307  double Minx=PmIn[0].Pos[CLat1];
308  int NOutShow=0;
309  for(int p=0;p<NIn;p++){
310  if(PmIn[p].Pos[CLat1] < Minx)
311  Minx = PmIn[p].Pos[CLat1];
312  if(PmIn[p].Pos[CLat1] > Maxx)
313  Maxx = PmIn[p].Pos[CLat1];
314  }
315  double DeltaxIn=(Maxx-Minx)/(double)(NIn-1);
316  double DeltaxOut=(Maxx-Minx)/(double)(NOut-1);
317  int NOrder = 3+1;
318  double *dArray = (double *)calloc(NIn+NOrder+1,sizeof(double));
319  for(int p=0;p<=NIn+NOrder;p++){
320  if(p<NOrder)
321  dArray[p] = Minx;
322  else if( (NOrder<=p) && (p<=NIn) )
323  dArray[p] = (p-NOrder+1)*DeltaxIn+Minx;
324  else if( p>NIn)
325  dArray[p] = (p-NOrder+2)*DeltaxIn+Minx;
326  }
327  for(int pp=0;pp<NOut-1;pp++){
328  for(int d=0;d<3;d++) PmOut[pp].Pos[d] = 0.;
329  double x = DeltaxOut*pp+Minx;
330  int vPos = (int)((x-Minx)/(Maxx-Minx)*NIn);
331  for(int p=vPos-1;p<vPos+NOrder;p++){
332  if(p >= NIn || p <0) continue;
333  double Blend = Mat->Blend(dArray,x,p,NOrder);
334  PmOut[pp].Pos[0] += Blend * PmIn[p].Pos[0];
335  PmOut[pp].Pos[1] += Blend * PmIn[p].Pos[1];
336  PmOut[pp].Pos[2] += Blend * PmIn[p].Pos[2];
337  }
338  }
339  PmOut[NOut-1].Pos[0] = PmIn[NIn-1].Pos[0];
340  PmOut[NOut-1].Pos[1] = PmIn[NIn-1].Pos[1];
341  PmOut[NOut-1].Pos[2] = PmIn[NIn-1].Pos[2];
342  NOutShow = NOut;
343  return NOutShow;
344 }
345 int VarData::InterBSpline2D(double **PlIn,double **PlOut,int NIn,int NOut){
346  double Max=Gen->Edge[CLat1];
347  double Min=0.;
348  int NOutShow=0;
349  double DeltaIn=(Max-Min)/(double)(NIn-1);
350  double DeltaOut=(Max-Min)/(double)(NOut-1);
351  double RatioInOut = (double)(NIn/(double)NOut);
352  int NOrder = 3+1;
353  double *dArray = (double *)calloc(NIn+NOrder+1,sizeof(double));
354  for(int p=0;p<=NIn+NOrder;p++){
355  if(p<NOrder){
356  dArray[p] = Min;
357  }
358  else if( (NOrder<=p) && (p<=NIn) ){
359  dArray[p] = (p-NOrder+1)*DeltaIn+Min;
360  }
361  else if(p>NIn){
362  dArray[p] = (p-NOrder+2)*DeltaIn+Min;
363  }
364  }
365  for(int vo=0;vo<NOut;vo++){
366  for(int vvo=0;vvo<NOut;vvo++){
367  PlOut[vo][vvo] = 0.;
368  double x = DeltaOut*vo+Min;
369  int vxPos = (int)(vo*RatioInOut);
370  for(int vi=vxPos-1,vn=0;vi<vxPos+NOrder+1;vi++){
371  vn = vi;
372  if(vi < 0) vn = NIn + vi;
373  if(vi >= NIn) vn = vi - NIn;
374  double Blendx = Mat->Blend(dArray,x,vn,NOrder);
375  double y = DeltaOut*vvo+Min;
376  int vyPos = (int)(vvo*RatioInOut);
377  for(int vvi=vyPos-1,vvn=0;vvi<vyPos+NOrder+1;vvi++){
378  vvn = vvi;
379  if(vvi < 0) vvn = NIn + vvi;
380  if(vvi >= NIn) vvn = vvi - NIn;
381  double Blendy = Mat->Blend(dArray,y,vvn,NOrder);
382  PlOut[vo][vvo] += Blendx*Blendy*PlIn[vn][vvn];
383  }
384  }
385  }
386  }
387  // for(int vo=0;vo<NOut;vo++){
388  // if(vo < NIn){//To be arranged
389  // PlOut[vo][NOut-1] = PlIn[vo][NIn-1];
390  // PlOut[NOut - 1][vo] = PlIn[NIn-1][vo];
391  // }
392  // }
393  //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]);
394  NOutShow = NOut;
395  free(dArray);
396  return NOutShow;
397 }
401 void VarData::SmoothGrid(int NSample){
402  double *PlotIn = (double *)calloc(SQR(NSample),sizeof(double));
403  double *PlotOut = (double *)calloc(SQR(NSample),sizeof(double));
404  for(int nx=0;nx<NSample;nx++){
405  for(int ny=0;ny<NSample;ny++){
406  int n = nx*NSample+ny;
407  PlotIn[n] = Pm[n].Pos[2];
408  }
409  }
410  InterBSpline2D(PlotIn,PlotOut,NSample,NSample);
411  for(int nx=0;nx<NSample;nx++){
412  for(int ny=0;ny<NSample;ny++){
413  int n = nx*NSample+ny;
414  if(pType(n) != 0) continue;
415  Pm[n].Pos[2] = PlotOut[n];
416  }
417  }
418  free(PlotIn);
419  free(PlotOut);
420 }
424 int VarData::InterBSpline2D(double *PlIn,double *PlOut,int NIn,int NOut){
425  double Max=Gen->Edge[CLat1];
426  double Min=0.;
427  int NOutShow=0;
428  double DeltaIn = (Max-Min)/(double)(NIn-1);
429  double DeltaOut = (Max-Min)/(double)(NOut-1);
430  double RatioInOut = (double)(NIn/(double)NOut);
431  int NOrder = 3+1;
432  double *dArray = (double *)calloc(NIn+NOrder+1,sizeof(double));
433  for(int p=0;p<=NIn+NOrder;p++){
434  if(p<NOrder){
435  dArray[p] = Min;
436  }
437  else if( (NOrder<=p) && (p<=NIn) ){
438  dArray[p] = (p-NOrder+1)*DeltaIn+Min;
439  }
440  else if(p>NIn){
441  dArray[p] = (p-NOrder+2)*DeltaIn+Min;
442  }
443  }
444  for(int vo=0;vo<NOut;vo++){
445  for(int vvo=0;vvo<NOut;vvo++){
446  PlOut[vo*NOut+vvo] = 0.;
447  double x = DeltaOut*vo+Min;
448  int vxPos = (int)(vo*RatioInOut);
449  for(int vi=vxPos-1,vn=0;vi<vxPos+NOrder+1;vi++){
450  vn = vi;
451  if(vi < 0) vn = NIn + vi;
452  if(vi >= NIn) vn = vi - NIn;
453  double Blendx = Mat->Blend(dArray,x,vn,NOrder);
454  double y = DeltaOut*vvo+Min;
455  int vyPos = (int)(vvo*RatioInOut);
456  for(int vvi=vyPos-1,vvn=0;vvi<vyPos+NOrder+1;vvi++){
457  vvn = vvi;
458  if(vvi < 0) vvn = NIn + vvi;
459  if(vvi >= NIn) vvn = vvi - NIn;
460  double Blendy = Mat->Blend(dArray,y,vvn,NOrder);
461  PlOut[vo*NOut+vvo] += Blendx*Blendy*PlIn[vn*NIn+vvn];
462  }
463  }
464  }
465  }
466  // for(int vo=0;vo<NOut;vo++){
467  // if(vo < NIn){//To be arranged
468  // PlOut[vo][NOut-1] = PlIn[vo][NIn-1];
469  // PlOut[NOut - 1][vo] = PlIn[NIn-1][vo];
470  // }
471  // }
472  //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]);
473  NOutShow = NOut;
474  free(dArray);
475  return NOutShow;
476 }
477 //minimum image convention does not seem to work
478 int VarData::InterBSpline1D(double *PlIn,double *PlOut,int NIn,int NOut){
479  double Max=Gen->Edge[CLat1];
480  double Min=0.;
481  int NOutShow=0;
482  double DeltaIn=(Max-Min)/(double)(NIn-1);
483  double DeltaOut=(Max-Min)/(double)(NOut-1);
484  int NOrder = 3+1;
485  double *dArray = (double *)calloc(NIn+2*NOrder+2,sizeof(double));
486  for(int vi=0;vi<=NIn+2*NOrder;vi++){
487  if(vi<NOrder){
488  dArray[vi] = Min;
489  }
490  else if(vi > NIn){
491  dArray[vi] = (vi-NOrder+2)*DeltaIn+Min;
492  }
493  else {
494  dArray[vi] = (vi-NOrder+1)*DeltaIn+Min;
495  }
496  }
497  for(int vo=0;vo<NOut;vo++){
498  //to avoid the minimum image convention
499  if(vo > NOut-NOrder) continue;
500  PlOut[vo] = 0.;
501  double x = DeltaOut*vo+Min;
502  for(int vi=vo-1,vn=0;vi<vo+NOrder+1;vi++){
503  vn = vi;
504  if(vi < 0){
505  vn = NIn + vi;
506  }
507  if(vi >= NIn){
508  vn = vi - NIn;
509  }
510  double Blendx = Mat->Blend(dArray,x,vn,NOrder);
511  PlOut[vo] += Blendx*PlIn[vn];
512  }
513  }
514  NOutShow = NOut;
515  free(dArray);
516  return NOutShow;
517 }
518 int VarData::InterPoly(PART *PmIn,PART *PmOut,int NIn,int NOut){
519  double *Pz = (double *)calloc(NIn,sizeof(double));
520  double *Px = (double *)calloc(NIn,sizeof(double));
521  double *Py = (double *)calloc(NIn,sizeof(double));
522  for(int p=0;p<NIn;p++){
523  Px[p] = PmIn[p].Pos[CLat1];
524  Py[p] = PmIn[p].Pos[CLat2];
525  Pz[p] = PmIn[p].Pos[CNorm];
526  }
527  Spline *Sp1 = new Spline(NIn);
528  Spline *Sp2 = new Spline(NIn);
529  Mat->Polinomio(Px,Py,NIn,Sp1);
530  Mat->Polinomio(Px,Pz,NIn,Sp2);
531  double Length = PmIn[NIn-1].Pos[CLat1] - PmIn[0].Pos[CLat1];
532  double Deltax = Length / (double)NOut;
533  double Pos = PmIn[0].Pos[CLat1];
534  for(int p=0;p<NOut;p++){
535  Pos += Deltax;
536  PmOut[p].Pos[CLat1] = Pos;
537  PmOut[p].Pos[CLat2] = 0.;
538  PmOut[p].Pos[CNorm] = 0.;
539  for(int n=0;n<NIn;n++){
540  PmOut[p].Pos[CLat2] += Sp1->GetCoe(n)*Mat->Elevato(Pos,n);
541  PmOut[p].Pos[CNorm] += Sp2->GetCoe(n)*Mat->Elevato(Pos,n);
542  }
543  //printf("(%lf %lf %lf)\n",PmOut[p].Pos[CLat1],PmOut[p].Pos[CLat2],PmOut[p].Pos[CNorm]);
544  }
545  delete Sp1;delete Sp2;
546  free(Pz);free(Py);free(Px);
547  return NOut;
548 }
549 // int VarData::InterDerMatrix(PART *Pm,int NMass,SPLINE Wg,double Offset){
550 // Matrice *Coeff = new Matrice(Wg,1);
551 // Matrice *Resp = new Matrice(NMass,NMass);
552 // for(int r=0;r<NMass;r++){
553 // if(Pm[r].Typ != 0){ Resp->Set(r,r,1.);continue;}
554 // if(r >= 2) Resp->Set(r,r-2,Coeff->Val(2,0));
555 // if(r >= 1) Resp->Set(r,r-1,Coeff->Val(2,1));
556 // if(r < NMass-1) Resp->Set(r,r+1,Coeff->Val(2,3));
557 // if(r < NMass-2) Resp->Set(r,r+2,Coeff->Val(2,4));
558 // Resp->Set(r,r,Coeff->Val(2,2));
559 // }
560 // double *Known = (double *) calloc(NMass,sizeof(double));
561 // double *UnKnown = (double *) calloc(NMass,sizeof(double));
562 // for(int p=0;p<NMass;p++)
563 // Known[p] = Pm[p].Pos[CNorm] - Offset;
564 // Resp->Solve(Known,UnKnown);
565 // for(int p=0;p<NMass;p++){
566 // if(Pm[p].Typ != 0){continue;}
567 // Pm[p].Pos[CNorm] = UnKnown[p] + Offset;
568 // //printf("%lf %lf\n",Known[p],UnKnown[p]);
569 // }
570 // free(Known);
571 // free(UnKnown);
572 // // Resp->invert();
573 // // for(int r=0;r<NMass;r++){
574 // // double Agg=0.;
575 // // for(int c=0;c<NMass;c++){
576 // // if(Pm[r].Typ == 0)
577 // // Agg += Resp->Val(r,c) * Pm[c].Pos[CNorm];
578 // // else
579 // // Agg += Resp->Val(r,c) * Pm[c].Pos[CNorm];
580 // // }
581 // // Pm[r].Pos[CNorm] = Agg;
582 // // }
583 // delete Coeff;
584 // delete Resp;
585 // return NMass;
586 // }
587 void VarData::SmoothGrid(int NSample,char *FWrite){
588  double InvNSample = 1./(double)NSample;
589  double **PlotIn = (double **)calloc(SQR(NSample),sizeof(double));
590  double **PlotOut = (double **)calloc(SQR(NSample),sizeof(double));
591  double **Count = (double **)calloc(SQR(NSample),sizeof(double));
592  double Round = 0.00001;
593  for(int v=0;v<NSample;v++){
594  PlotIn[v] = (double *)calloc(NSample,sizeof(double));
595  PlotOut[v] = (double *)calloc(NSample,sizeof(double));
596  Count[v] = (double *)calloc(NSample,sizeof(double));
597  }
598  int Nx = 0;
599  int Ny = 0;
600  for(int p=0;p<pNPart();p++){
601  if(Pm[p].Pos[CLat2] < Pm[p+1].Pos[CLat2]) Ny++;
602  else break;
603  }
604  Ny++;
605  Nx = (int)(pNPart()/(double)Ny);
606  for(int p=0;p<pNPart();p++){
607  int vx = (int)((Pm[p].Pos[CLat1]+Round)*pInvEdge(CLat1)*NSample);
608  int vy = (int)((Pm[p].Pos[CLat2]+Round)*pInvEdge(CLat2)*NSample);
609  if(vx < 0 || vx >= NSample)continue;
610  if(vy < 0 || vy >= NSample)continue;
611  PlotIn[vx][vy] += Pm[p].Pos[CNorm];
612  Count[vx][vy] += 1.;
613  }
614  for(int vx=0;vx<NSample;vx++)
615  for(int vy=0;vy<NSample;vy++)
616  PlotIn[vx][vy] *= Count[vx][vy] > 0. ? 1./Count[vx][vy] : 1.;
617  InterBSpline2D(PlotIn,PlotOut,NSample,NSample);
618  if(1==0){
619  FILE *F2Write = fopen(FWrite,"w");
620  for(int vx=0;vx<NSample;vx++){
621  for(int vy=0;vy<NSample;vy++){
622  double x = vx*pEdge(CLat1)*InvNSample;
623  double y = vy*pEdge(CLat2)*InvNSample;
624  fprintf(F2Write,"%lf %lf %lf\n",x,y,PlotOut[vx][vy]);
625  }
626  }
627  fclose(F2Write);
628  }
629  for(int v=0;v<NSample;v++){
630  free(PlotIn[v]);
631  free(PlotOut[v]);
632  free(Count[v]);
633  }
634  free(PlotIn);
635  free(PlotOut);
636  free(Count);
637 }
int InterSpline3(PART *PmIn, PART *PmOut, int NIn, int NOut)
third order spline
SPLINE Spline4Beg(double *P1, double *P2, double *P3, double *P4, int x, int y)
Four order spline first boundary.
double a4
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
double GetCoe(int n)
Print the n coefficient.
int InterBSpline1D(double *PlIn, double *PmOut, int NIn, int NOut)
1-d BSpline
int InterSpline4(PART *PmIn, PART *PmOut, int NIn, int NOut)
forth order spline
coeficients of a N order spline
double Edge[4]
xyzr edges of the simulation box
Definition: VarData.h:309
SPLINE Spline3Beg(double *P1, double *P2, double *P3, int x, int y)
Three order spline first boundary.
double pInvEdge(int d)
Inverted xyzr edges of the simulation box.
Definition: VarData.h:920
int pType(int p)
Return the type.
int CLat2
lateral coordinate
Definition: VarData.h:1078
double Pos[3]
xyz Position of the particle
Definition: VarData.h:216
double pEdge(int d)
xyzr edges of the simulation box
Definition: VarData.h:918
SPLINE Parab2(double *PA, double *PB, double *PC, int x, int y)
Three points parabolic interpolation.
int Idx
Particle identifier.
Definition: VarData.h:222
double a3
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
SPLINE Spline4(double *P1, double *P2, double *P3, double *P4, int x, int y)
Four order spline.
int InterBSpline(PART *PmIn, PART *PmOut, int NIn, int NOut)
BSpline.
int InterParab(PART *PmIn, PART *PmOut, int NIn, int nOut)
Discontinous parabolas.
int InterParab2(PART *PmIn, PART *PmOut, int NIn, int NOut)
Discontinous parabolas.
int InterPoly(PART *PmIn, PART *PmOut, int NIn, int nOut)
NIn-polynomian.
double a2
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
double a1
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
int InterForth(PART *PmIn, PART *PmOut, int NIn, int NOut)
Discontinous forth degree.
int Polinomio(double *P1, double *P2, int NMass, Spline *Sp)
Polinimial interpolation of NMass order.
SPLINE Spline3End(double *P1, double *P2, int x, int y)
Three order spline last boundary.
double Elevato(double x, int Volte)
Integer power.
void SmoothGrid(int NSample, char *FWrite)
Smooth a grid with BSplines.
Matematica * Mat
Implementation of all usefull algorythms.
Definition: VarData.h:527
SPLINE Parab(double *P1, double *P2, double *P3, int x, int y)
Three points parabolic interpolation.
int CLat1
lateral coordinate
Definition: VarData.h:1076
Information of every particle.
Definition: VarData.h:214
SPLINE Spline3(double *P1, double *P2, double *P3, int x, int y)
Three order spline.
int InterCubica(PART *PmIn, PART *PmOut, int NIn, int NOut)
Discontinous cubic.
double a0
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
Coefficient of a spline.
SPLINE Forth(double *PA, double *PB, double *PC, double *PD, double *PE, int x, int y)
Five points four order interpolation.
SPLINE Cubica(double *PA, double *PB, double *PC, double *PD, int x, int y)
Four point cubic interpolation.
double Blend(const double *dPoint, double x, int nPoint, int nOrder)
For the BSpline.
PART * Pm
Particle information of all particle.
Definition: VarData.h:1046
int InterBSpline2D(double **PlIn, double **PmOut, int NIn, int NOut)
2-d BSpline
int CNorm
Normal coordinate.
Definition: VarData.h:1074
int pNPart()
Number of particle.