Allink  v0.1
MatematicaInterp.cpp
1 #include "../include/Matematica.h"
2 
3 double Matematica::QBezier(double *P1,double *P2,double *P3,double x,int y){
4  double Resp;
5 // int Order = 3;
6 // double *Coeff;
7 // Coeff = (double *)malloc(Order*sizeof(double));
8 // for(int i=0;i<Order;i++){
9 // Coeff[i] = Binomial(Order,i)*P1[y];
10 // }//deCasteljau algorithm
11 // //P(x) = sum_i Coeff[i]*t^i*(t-1)^(n-1)
12  Resp = QUAD((1-x))*P1[y]+2*(1-x)*x*P2[y]+QUAD(x)*P3[y];
13  return Resp;
14 }
15 SPLINE Matematica::Spline3Beg(double *P1,double *P2,double *P3,int x,int y){
16  SPLINE Sp; Sp.a0 = 0.;
17  if(x < 0 || x>3 || y < 0 || y>3){
18  printf("x,y must be 0 <= %d , %d < 3\n",x,y);
19  return Sp;
20  }
21  double Deltax = P2[x] - P1[x];
22  Sp.a0 = P1[y];
23  SPLINE Par = Parab2(P1,P2,P3,x,y);
24  Sp.a2 = 0.;
25  // printf("%lf %lf\n",Sp.a2,Par.A);
26  Sp.a3 = (Par.a2 - Sp.a2)/(3.*Deltax);
27  Sp.a1 = (P2[y] - P1[y])/Deltax - Sp.a2*Deltax - Sp.a3*Deltax*Deltax;
28  SpMem.a2 = Par.a2;
29  SpMem.a1 = Sp.a1+2*Sp.a2*Deltax+3*Sp.a3*Deltax*Deltax;
30  Sp.a4 = 0.;
31  //printf("%lf %lf|%lf %lf %lf\n",Sp.a0,Sp.a1,SplineA1,Sp.a2,Sp.a3);
32  return Sp;
33 }
34 SPLINE Matematica::Spline3(double *P1,double *P2,double *P3,int x,int y){
35  SPLINE Sp; Sp.a0 = 0.;
36  if(x < 0 || x>3 || y < 0 || y>3){
37  printf("x,y must be 0 <= %d , %d < 3\n",x,y);
38  return Sp;
39  }
40  double Deltax = P2[x] - P1[x];
41  SPLINE Par = Parab2(P1,P2,P3,x,y);
42  Sp.a0 = P1[y];
43  Sp.a2 = SpMem.a2;
44  Sp.a3 = (Par.a2 - Sp.a2)/(3.*Deltax);
45  Sp.a1 = (P2[y] - P1[y])/Deltax - Sp.a2*Deltax - Sp.a3*Deltax*Deltax;
46  SpMem.a2 = Par.a2;
47  //printf("%lf %lf|%lf %lf|%lf %lf\n",Sp.a0,SplineA1,Sp.a1,Sp.a2,SplineA2,Sp.a3);
48  SpMem.a1 = Sp.a1;
49  Sp.a4 = 0.;
50  return Sp;
51 }
52 SPLINE Matematica::Spline3End(double *P1,double *P2,int x,int y){
53  SPLINE Sp; Sp.a0 = 0.;
54  if(x < 0 || x>3 || y < 0 || y>3){
55  printf("x,y must be 0 <= %d , %d < 3\n",x,y);
56  return Sp;
57  }
58  double Deltax = P2[x] - P1[x];
59  Sp.a0 = P1[y];
60  Sp.a2 = SpMem.a2;
61  Sp.a3 = (0. - SpMem.a2)/(3.*Deltax);
62  Sp.a1 = (P2[y] - P1[y])/Deltax - Sp.a2*Deltax - Sp.a3*Deltax*Deltax;
63  // printf("%lf %lf %lf %lf\n",Sp.a0,Sp.a1,Sp.a2,Sp.a3);
64  SpMem.a2 = 0.;
65  Sp.a4 = 0.;
66  return Sp;
67 }
68 SPLINE Matematica::Spline4Beg(double *P1,double *P2,double *P3,double *P4,int x,int y){
69  SPLINE Sp; Sp.a0 = 0.;
70  if(x < 0 || x>3 || y < 0 || y>3){
71  printf("x,y must be 0 <= %d , %d < 3\n",x,y);
72  return Sp;
73  }
74  double Deltax = P2[x] - P1[x];
75  SPLINE Cub = Cubica(P1,P2,P3,P4,x,y);
76  SPLINE Par = Parab2(P1,P2,P3,x,y);
77  Sp.a0 = P1[y];
78  Sp.a2 = 0.;
79  Sp.a4 = 3.*Cub.a3/Deltax + (Sp.a2 - Cub.a2)/QUAD(Deltax);
80  Sp.a4 /= 12.;
81  Sp.a3 = Cub.a3 - 6.*Sp.a4*Deltax;
82  Sp.a1 = (P2[y] - P1[y])/Deltax - Sp.a2*Deltax - Sp.a3*Deltax*Deltax - Sp.a4*Deltax*Deltax*Deltax;
83  SpMem.a2 = Cub.a2;
84  SpMem.a3 = Sp.a3;
85  return Sp;
86 }
87 SPLINE Matematica::Spline4(double *P1,double *P2,double *P3,double *P4,int x,int y){
88  SPLINE Sp; Sp.a0 = 0.;
89  if(x < 0 || x>3 || y < 0 || y>3){
90  printf("x,y must be 0 <= %d , %d < 3\n",x,y);
91  return Sp;
92  }
93  double Deltax = P2[x] - P1[x];
94  SPLINE Cub = Cubica(P1,P2,P3,P4,x,y);
95  SPLINE Par = Parab2(P1,P2,P3,x,y);
96  Sp.a0 = P1[y];
97  Sp.a2 = SpMem.a2;
98  Sp.a4 = 3.*Cub.a3/Deltax + (Sp.a2 - Cub.a2)/QUAD(Deltax);
99  Sp.a4 /= 12.;
100  Sp.a3 = Cub.a3 - 6.*Sp.a4*Deltax;
101  Sp.a1 = (P2[y] - P1[y])/Deltax - Sp.a2*Deltax - Sp.a3*Deltax*Deltax - Sp.a4*Deltax*Deltax*Deltax;
102  SpMem.a2 = Cub.a2;
103  // printf("%lf %lf|%lf %lf|%lf %lf %lf\n",Sp.a0,SpMem.a1,Sp.a1,Sp.a2,SpMem.a2,Sp.a3,Sp.a4);
104  SpMem.a1 = Sp.a1;
105  SpMem.a3 = Cub.a3;
106  return Sp;
107 }
108 SPLINE Matematica::Spline4(double *P1,double *P2,double *P3,int x,int y){
109  SPLINE Sp; Sp.a0 = 0.;
110  if(x < 0 || x>3 || y < 0 || y>3){
111  printf("x,y must be 0 <= %d , %d < 3\n",x,y);
112  return Sp;
113  }
114  double Deltax = P2[x] - P1[x];
115  SPLINE Par = Parab2(P1,P2,P3,x,y);
116  Sp.a0 = P1[y];
117  Sp.a2 = SpMem.a2;
118  Sp.a1 = SpMem.a1;
119 // Sp.a1 = Par.B - 13.*Sp.a2*Deltax/2. - 3.*(P2[y] - P1[y])/Deltax - Par.A*Deltax/2.;
120 // Sp.a1 /= 13.;
121  Sp.a4 = Par.a2 + 2.*Sp.a2 - 3.*(P2[y] - P1[y])/(Deltax*Deltax) + 3.*Sp.a1/Deltax;
122  Sp.a4 /= 3.;
123  Sp.a3 = (P2[y] - P1[y])/(Deltax*Deltax*Deltax) - Sp.a1/(Deltax*Deltax) - Sp.a2/Deltax - Sp.a4*Deltax;
124  SpMem.a2 = Par.a2;
125  //printf("%lf %lf|%lf %lf|%lf %lf %lf\n",Sp.a0,SpMem.a1,Sp.a1,Sp.a2,SpMem.a2,Sp.a3,Sp.a4);
126  SpMem.a1 = Par.a1;
127  return Sp;
128 }
129 SPLINE Matematica::Spline4PreEnd(double *P1,double *P2,double *P3,int x,int y){
130  SPLINE Sp; Sp.a0 = 0.;
131  if(x < 0 || x>3 || y < 0 || y>3){
132  printf("x,y must be 0 <= %d , %d < 3\n",x,y);
133  return Sp;
134  }
135  double Deltax = P2[x] - P1[x];
136  SPLINE Par = Parab2(P1,P2,P3,x,y);
137  Sp.a0 = P1[y];
138  Sp.a2 = SpMem.a2;
139  Sp.a3 = (Sp.a2 - Par.a2)/Deltax + 4.*SpMem.a3;
140  Sp.a4 = (SpMem.a3 - Sp.a3)/(6.*Deltax);
141  Sp.a1 = (P2[y] - P1[y])/Deltax - Sp.a2*Deltax - Sp.a3*Deltax*Deltax - Sp.a4*Deltax*Deltax*Deltax;
142  SpMem.a2 = Par.a2;
143  SpMem.a1 = Sp.a1;
144  SpMem.a3 = Sp.a3;
145  return Sp;
146 }
147 SPLINE Matematica::Spline4End(double *P1,double *P2,int x,int y){
148  SPLINE Sp; Sp.a0 = 0.;
149  if(x < 0 || x>3 || y < 0 || y>3){
150  printf("x,y must be 0 <= %d , %d < 3\n",x,y);
151  return Sp;
152  }
153  double Deltax = P2[x] - P1[x];
154  Sp.a0 = P1[y];
155  Sp.a2 = SpMem.a2;
156  //Sp.a3 = (Sp.a2 - 0.)/Deltax + 4.*SpMem.a3;
157  Sp.a3 = 0.;
158  Sp.a4 = (SpMem.a3 - Sp.a3)/(6.*Deltax);
159  Sp.a1 = (P2[y] - P1[y])/Deltax - Sp.a2*Deltax - Sp.a3*Deltax*Deltax - Sp.a4*Deltax*Deltax*Deltax;
160  SpMem.a2 = 0.;
161  SpMem.a1 = Sp.a1;
162  SpMem.a3 = Sp.a3;
163  return Sp;
164 }
165 //BSplines
166 // double w = x - floor(x) - 1.;
167 // xSp.a3 = (1./6.)*w * w * w;
168 // xSp.a0 = (1./6.) + (1./2.)*w*(w-1.)-xSp.a3;
169 // xSp.a2 = w + xSp.a0 - 2*xSp.a3;
170 // xSp.a1 = 1. - xSp.a0 - xSp.a2 - xSp.a3;
171 // Resp = xSp.a0*P1[0]+xSp.a1*P2[0]+xSp.a2*P3[0];
172  // w = y - 3.;
173 // ySp.a3 = (1./6.)*w * w * w;
174 // ySp.a0 = (1./6.) + (1./2.)*w*(w-1.)-ySp.a3;
175 // ySp.a2 = w + ySp.a0 - 2*ySp.a3;
176 // ySp.a1 = 1. - ySp.a0 - ySp.a2 - ySp.a3;
177  // Resp *= ySp.a0*P1[1]+ySp.a1*P2[1]+ySp.a2*P3[1];
178 SPLINE Matematica::Parab(double *P1,double *P2,double *P3,int x,int y){
179  SPLINE Par; Par.a0 = 0.;
180  if(x < 0 || x>3 || y < 0 || y>3){
181  printf("x,y must be 0 <= %d , %d < 3\n",x,y);
182  return Par;
183  }
184  double Deltax = P2[x] - P1[x];
185  double Deltax2 = QUAD(P2[x]) - QUAD(P1[x]);
186  double Dx = P3[x] - P1[x];
187  double Dx2 = QUAD(P3[x]) - QUAD(P1[x]);
188  double Deltay = P2[y] - P1[y];
189  double Dy = P3[y] - P1[y];
190  double b1= Dy - Deltay/Deltax2*Dx2;
191  double b2= Dx - Deltax/Deltax2*Dx2;
192  Par.a1 = b1/b2;
193  Par.a2 = Deltay/Deltax2 - Par.a1*Deltax/Deltax2;
194  Par.a0 = P1[y] - Par.a1*P1[x] - Par.a2*QUAD(P1[x]);
195  Par.a3 = 0.;
196  Par.a4 = 0.;
197  return Par;
198 }
199 SPLINE Matematica::Parab2(double *PA,double *PB,double *PC,int x,int y){
200  SPLINE Par; Par.a0 = 0.;
201  if(x < 0 || x>3 || y < 0 || y>3){
202  printf("x,y must be 0 <= %d , %d < 3\n",x,y);
203  return Par;
204  }
205  double DxAB = PA[x] - PB[x];
206  double DxCB = PC[x] - PB[x];
207  double DyCB = PC[y] - PB[y];
208  double DyAB = PA[y] - PB[y];
209  double a1 = DyCB * DxAB - DyAB * DxCB;
210  double a2 = DxCB*DxCB*DxAB - DxCB*DxAB*DxAB;
211  Par.a2 = a1/a2;
212  Par.a1 = DyAB/DxAB - Par.a2*DxAB;
213  Par.a0 = PB[y];
214  Par.a3 = 0.;
215  Par.a4 = 0.;
216  return Par;
217 }
218 CIRCLE Matematica::Osculante(double *PA,double *PB,double *PC,int x,int y){
219  CIRCLE Cir; Cir.yC = 0.;
220  if(x < 0 || x>3 || y < 0 || y>3){
221  printf("x,y must be 0 <= %d , %d < 3\n",x,y);
222  return Cir;
223  }
224  double DxCA = PC[x] - PA[x];
225  double DxBA = PB[x] - PA[x];
226  double DyAC = PA[y] - PC[y];
227  double DyAB = PA[y] - PB[y];
228  double DxAB2 = QUAD(PA[x]) - QUAD(PB[x]);
229  double DxAC2 = QUAD(PA[x]) - QUAD(PC[x]);
230  double DyBA2 = QUAD(PB[y]) - QUAD(PA[y]);
231  double DyCA2 = QUAD(PC[y]) - QUAD(PA[y]);
232  double a1 = DxCA*(DyBA2 - DxAB2) - DxBA*(DyCA2 - DxAB2);
233  double a2 = 2.*(DyAC*DxBA - DyAB*DxCA);
234  Cir.yC = a1/a2;
235  Cir.xC = (DyBA2 - DxAB2 + 2.*Cir.yC*DyAB)/DxBA;
236  Cir.Rad = sqrt( QUAD((PA[x] - Cir.xC)) + QUAD((PA[y] - Cir.yC)) );
237  return Cir;
238 }
239 SPLINE Matematica::Cubica(double *PA,double *PB,double *PC,double *PD,int x,int y){
240  SPLINE Cub; Cub.a0 = 0.;
241  if(x < 0 || x>3 || y < 0 || y>3){
242  printf("x,y must be 0 <= %d , %d < 3\n",x,y);
243  return Cub;
244  }
245  double DyDB = PD[y] - PB[y];
246  double DyCB = PC[y] - PB[y];
247  double DyAB = PA[y] - PB[y];
248  double DxDB = PD[x] - PB[x];
249  double DxCB = PC[x] - PB[x];
250  double DxAB = PA[x] - PB[x];
251  double DxCA = PC[x] - PA[x];
252  double Num3 = DyDB/DxDB - DyAB/DxAB - (DyCB/DxCB - DyAB/DxAB)*(DxDB - DxAB)/DxCA;
253  double Den3 = DxDB*DxDB - DxAB*DxAB - (DxCB*DxCB - DxAB*DxAB)*(DxDB - DxAB)/DxCA;
254  Cub.a3 = Num3/Den3;
255  Cub.a2 = DyCB/DxCB - DyAB/DxAB - Cub.a3*(DxCB*DxCB - DxAB*DxAB);
256  Cub.a2 /= DxCA;
257  Cub.a1 = DyAB/DxAB - Cub.a2*DxAB - Cub.a3*DxAB*DxAB;
258  Cub.a0 = PB[y];
259  Cub.a4 = 0.;
260  return Cub;
261 }
262 SPLINE Matematica::Forth(double *PA,double *PB,double *PC,double *PD,double *PE,int x,int y){
263  SPLINE Forth; Forth.a0=0.;
264  if(x < 0 || x>3 || y < 0 || y>3){
265  printf("x,y must be 0 <= %d , %d < 3\n",x,y);
266  return Forth;
267  }
268  double Dx = ASS((PC[x] - PB[x]));
269  Forth.a3 = PA[y]/12. - PB[y]/6. + PD[y]/6. - PE[y]/12.;
270  Forth.a3 /= -Dx*Dx*Dx;
271  Forth.a1 = - PA[y]/12. + 4.*PB[y]/6. - 4.*PD[y]/6. + PE[y]/12.;
272  Forth.a1 /= -Dx;
273  Forth.a4 = PA[y]/24. - PB[y]/6. + 1./4.*PC[y] - PD[y]/6. + PE[y]/24.;
274  Forth.a4 /= Dx*Dx*Dx*Dx;
275  Forth.a2 = -PA[y]/24. + 2./3.*PB[y] - 5./4.*PC[y] + 2./3.*PD[y] - PE[y]/24.;
276  Forth.a2 /= Dx*Dx;
277  Forth.a0 = PC[y];
278  //printf("Dx %lf Coeff %lf %lf %lf %lf %lf \n",Dx,Forth.a0,Forth.a1,Forth.a2,Forth.a3,Forth.a4);
279  return Forth;
280 
281  double DxCA = PC[x] - PA[x];
282  double DxCB = PC[x] - PB[x];
283  double DxCD = PC[x] - PD[x];
284  double DxCE = PC[x] - PE[x];
285  double DxAB = PA[x] - PB[x];
286  double DxAD = PA[x] - PD[x];
287  double DyAC = (PA[y] - PC[y])/DxCA;
288  double DyBC = (PB[y] - PC[y])/DxCB;
289  double DyDC = (PD[y] - PC[y])/DxCD;
290  double DyEC = (PE[y] - PC[y])/DxCE;
291  //if(DxCB == DxCD)
292  if (1==1)
293  {
294  double DxDB = PD[x] - PB[x];
295  double DxEA = PE[x] - PA[x];
296  double DyBD = (PB[y] - PD[y])/DxDB;
297  double DyAE = (PA[y] - PE[y])/DxEA;
298  // Forth.a4 = - Noto34/Num34;
299  Forth.a4 = DyBC - DyDC - DyEC + DyAC;
300  Forth.a4 /= (DxCE*DxCE*DxCE - DxCA*DxCA*DxCA) - (DxCB*DxCB*DxCB - DxCD*DxCD*DxCD);
301  Forth.a2 = DyBC - DyDC - Forth.a4*(DxCB*DxCB*DxCB - DxCD*DxCD*DxCD);
302  Forth.a3 = DyBD - DyAE;
303  Forth.a3 /= (DxCB*DxCB*DxCB - DxCD*DxCD*DxCD)/DxDB -
304  (DxCA*DxCA*DxCA - DxCE*DxCE*DxCE)/DxEA;
305  Forth.a1 = DyAE - Forth.a3*(DxCA*DxCA*DxCA - DxCE*DxCE*DxCE)/DxEA;
306  }
307  else
308  {
309  double DxAE = PA[x] - PE[x];
310  double Den34 = DxAB*(DxCB*DxCB - DxCA*DxCA) - DxAB*(DxCD*DxCD - DxCA*DxCA);// 0.;
311  double Num34 = DxAB*(DxCD*DxCD*DxCD - DxCA*DxCA*DxCA) - DxAD*(DxCB*DxCB*DxCB - DxCA*DxCA*DxCA);//21Dx
312  double Noto34 = DyBC*DxAD - DyAC*DxAD - DyDC*DxAB + DyAC*DxAB;
313  double Num4 = -DxAB*(Noto34*(DxCE*DxCE - DxCA*DxCA) - DyEC + DyAC) +
314  DxAE*(Noto34*(DxCB*DxCB - DxCA*DxCA) - DyBC + DyAC);
315  double Den4 = DxAB*(Num34/Den34*(DxCE*DxCE - DxCA*DxCA) + (DxCE*DxCE*DxCE - DxCA*DxCA*DxCA) ) -
316  DxAE*(Num34/Den34*(DxCB*DxCB - DxCA*DxCA) + (DxCB*DxCB*DxCB - DxCA*DxCA*DxCA));
317  //printf("Dx %lf %lf %lf %lf %lf %lf %lf\n",DxCB,Den34,Num34,Noto34,Num4,Den4,Num34/Den34);
318  Forth.a4 = Num4/Den4;
319  Forth.a3 = Num34/Den34*Forth.a4 + Noto34;
320  Forth.a2 = DyBC - DyAC - Forth.a3*(DxCB*DxCB - DxCA*DxCA) - Forth.a4* (DxCB*DxCB*DxCB - DxCA*DxCA*DxCA);
321  Forth.a2 /= DxAB;
322  Forth.a1 = DyAC - Forth.a2*DxCA - Forth.a3*DxCA*DxCA - Forth.a4*DxCA*DxCA*DxCA;
323  }
324  Forth.a0 = PC[y];
325  //printf("%lf %lf %lf %lf %lf \n",Forth.a0,Forth.a1,Forth.a2,Forth.a3,Forth.a4);
326  return Forth;
327 }
328 int Matematica::Polinomio(double *Px,double *Py,int NMass,Spline *Sp){
329  Matrice *Coeff = new Matrice(NMass,NMass);
330  for(int r=0;r<NMass;r++){
331  Sp->SetCoe( 0. , r);
332  for(int c=0;c<NMass;c++){
333  Coeff->Set(r,c,Elevato(Px[r],c));
334  }
335  }
336  Coeff->Invert();
337  for(int c=0;c<NMass;c++){
338  for(int r=0;r<NMass;r++){
339  Sp->AddCoe( Coeff->Val(r,c)*Py[c] , r);
340  }
341  }
342  delete Coeff;
343  return 1;
344 }
345 int Matematica::DerMatrix(double *Px,double *Py,int NMass,SPLINE Wg,Spline *Sp){
346  Matrice *Coeff = new Matrice(NMass,NMass);
347  double MenoDue = - .75*Wg.a3 + 1.5*Wg.a4;
348  //MenoDue += .125*Wg.a1 - .125*Wg.a2; //O(h^4)
349  double MenoUno = -.5*Wg.a1 + Wg.a2 + 1.5*Wg.a3 - 6.*Wg.a4;
350  //MenoUno += -.125*Wg.a1 - .5*Wg.a2;
351  double Zero = Wg.a0 - 2.*Wg.a2 + 9.*Wg.a4;
352  //Zero += 0.75*Wg.a2;
353  double PiuUno = .5*Wg.a1 + Wg.a2 - 1.5*Wg.a3 - 6.*Wg.a4;
354  //PiuUno += .25*Wg.a1 + .5*Wg.a2;
355  double PiuDue = .75*Wg.a3 + 1.5*Wg.a4;
356  //PiuDue += -.125*Wg.a1 -.125*Wg.a2;
357  for(int r=0;r<NMass;r++){
358  if(r> 1) Coeff->Set(r,r-2,MenoDue*Px[r-2]);
359  if(r> 0) Coeff->Set(r,r-1,MenoUno*Px[r-1]);
360  if(r< NMass-1) Coeff->Set(r,r+1,PiuUno*Px[r+1]);
361  if(r< NMass-2) Coeff->Set(r,r+2,PiuDue*Px[r+2]);
362  Coeff->Set(r,r,Zero*Px[r]);
363  }
364  //Coeff->Print();
365  Coeff->Invert();
366  for(int r=0;r<NMass;r++){
367  for(int c=0;c<NMass;c++){
368  Sp->AddCoe( Coeff->Val(r,c)*Py[c] , r);
369  }
370  }
371  delete Coeff;
372  return 0;
373 }
374 double Matematica::LinInterp(double Px1,double Px2,double Py1,double Py2,double x){
375  double m = (Py2-Py1)/(Px2-Px1);
376  double q = Py1 - m*Px1;
377  return m*x+q;
378 }
379 RETTA Matematica::InterRett(double *Px,double *Py,int NMass){
380  RETTA r1;
381  double Uno=0.;double Due=0.;double Tre=0.;double Quattro=0.;
382  double UnoUno=0.;double DueUno=0.;double ZeroUno=0.;double ZeroDue=0.;
383  for(int i=0;i<NMass;i++){
384  //printf("%lf %lf\n",Px[i],Py[i]);
385  Uno += Px[i];
386  Due += QUAD(Px[i]);
387  ZeroUno += Py[i];
388  UnoUno += Px[i]*Py[i];
389  ZeroDue += QUAD(Py[i]);
390  }
391  double Mediax = Uno / (double) NMass;
392  double Mediay = ZeroUno / (double) NMass;
393  double Scartox = (Due - NMass*Uno*Uno)/(double)(NMass-0);
394  double Scartoy = (ZeroDue - NMass*ZeroUno*ZeroUno)/(double)(NMass-0);
395  r1.m = (NMass*UnoUno - Uno*ZeroUno) / (NMass*Due - Uno*Uno);
396  r1.q = (ZeroUno - r1.m*Uno)/NMass;
397  double Posteriori = 0.;
398  r1.Cov = UnoUno/(double)NMass - Mediax*Mediay;
399  for(int i=0;i<NMass;i++){
400  Posteriori += QUAD(( Px[i]*r1.m + r1.q - Py[i] ));
401  }
402  r1.Corr = r1.Cov / ( sqrt(Scartox*Scartoy) );
403  r1.ErrY = sqrt(Posteriori/(double)(NMass-2));
404  r1.ErrM = r1.ErrY * sqrt( NMass / (NMass * Due - Uno*Uno));
405  r1.ErrQ = r1.ErrY * sqrt( Due / (NMass*Due - Uno*Uno) );
406  // printf("m %lf q %lf r %lf sigma %lf Mediax %lf Mediay %lf\n",r1.m,r1.q,r1.r,r1.Corr,Mediax,Mediay);
407  return r1;
408 }
409 // Vettore Matematica::Directive(double *Px,int NMass){
410 // Vettore v1(3);
411 
412 
413 // }
414 RETTA Matematica::InterRett(double *Px,double *Py,double *Peso,int NMass){
415  RETTA r1;
416  double Mediax = 0.;
417  double Mediay = 0.;
418  double Scartox = 0.;
419  double Scartoy = 0.;
420  double Pesox = 0.;
421  double Uno=0.;double Due=0.;double Tre=0.;double Quattro=0.;
422  double UnoUno=0.;double DueUno=0.;double ZeroUno=0.;double ZeroDue=0.;
423  for(int i=0;i<NMass;i++){
424  Pesox += Peso[i];
425  Uno += Px[i]*Peso[i];
426  Due += QUAD(Px[i]*Peso[i]);
427  ZeroUno += Py[i];
428  UnoUno += Px[i]*Py[i]*Peso[i];
429  ZeroDue += QUAD(Py[i]);
430  }
431  Mediax = Uno / (double) NMass/Pesox;
432  Mediay = ZeroUno / (double) NMass;
433  Scartox = (Due - Uno*Uno/QUAD(Pesox))/(double)NMass;
434  Scartoy = (ZeroDue - ZeroUno*ZeroUno)/(double)NMass;
435  r1.Corr = (UnoUno - Uno*ZeroUno) / (double) NMass/Pesox ;
436  r1.r = r1.Corr / ( sqrt(Scartox*Scartoy) );
437  r1.m = (NMass*UnoUno - Uno*ZeroUno) / (NMass*Due - Uno*Uno)/Pesox;
438  r1.ErrM = Scartoy * sqrt( Due / (NMass * Due - Uno*Uno)/QUAD(Pesox));
439  r1.q = (ZeroUno - r1.m*Uno)/NMass;
440  r1.ErrQ = Scartoy * sqrt( NMass / (NMass * Due - Uno*Uno)/QUAD(Pesox) );
441  return r1;
442 }
443 RETTA Matematica::InterExp(double *Px,double *Py,int NMass){
444  RETTA r1;
445  double Uno=0.;double Due=0.;double Tre=0.;double Quattro=0.;
446  double UnoUno=0.;double DueUno=0.;double ZeroUno=0.;double ZeroDue=0.;
447  for(int i=0;i<NMass;i++){
448  if(Py[i] <= 0.){
449  //printf("Negative number not allowed for an exponential interpolation %lf\n",Py[i]);
450  continue;
451  }
452  Uno += Px[i];
453  Due += QUAD(Px[i]);
454  ZeroUno += log10(Py[i]);
455  UnoUno += Px[i]*log10(Py[i]);
456  ZeroDue += QUAD(log10(Py[i]));
457  }
458  double Mediax = Uno / (double) NMass;
459  double Mediay = ZeroUno / (double) NMass;
460  double Scartox = (Due - NMass*Uno*Uno)/(double)(NMass-1);
461  double Scartoy = (ZeroDue - NMass*ZeroUno*ZeroUno)/(double)(NMass-1);
462  r1.m = (NMass*UnoUno - Uno*ZeroUno) / (NMass*Due - Uno*Uno);
463  r1.q = (ZeroUno - r1.m*Uno)/NMass;
464  double Posteriori = 0.;
465  r1.Cov = UnoUno/(double)NMass - Mediax*Mediay;
466  for(int i=0;i<NMass;i++){
467  Posteriori += QUAD(( Px[i]*r1.m + r1.q - log10(Py[i]) ));
468  }
469  r1.Corr = r1.Cov / ( sqrt(Scartox*Scartoy) );
470  r1.ErrY = sqrt(Posteriori/(double)(NMass-2));
471  r1.ErrM = r1.ErrY * sqrt( NMass / (NMass * Due - Uno*Uno));
472  r1.ErrQ = r1.ErrY * sqrt( Due / (NMass*Due - Uno*Uno) );
473  // printf("m %lf q %lf r %lf sigma %lf Mediax %lf Mediay %lf\n",r1.m,r1.q,r1.r,r1.Corr,Mediax,Mediay);
474  return r1;
475 }
476 MOMENTI Matematica::InterGauss(double *Px,double *Py,int NMax){
477  MOMENTI m1; m1.Uno=0.; m1.Due=0.; m1.Tre=0.;m1.Delta=0.;m1.Num=0;
478  m1.Min = Px[0];
479  m1.Max = Px[0];
480  double Count = 0.;
481  double Sum2 = 0.;
482  for(int i=0;i<NMax;i++){
483 #ifdef MAT_DEBUG
484  if(Py[i] < 0.){printf("Invalid weight %lf < 0\n",Py[i]);}
485 #endif
486  if(Px[i]<m1.Min) m1.Min = Px[i];
487  if(Px[i]>m1.Max) m1.Max = Px[i];
488  m1.Uno += Px[i]*Py[i];
489  Sum2 += SQR(Px[i])*Py[i];
490  Count += Py[i];
491  }
492  m1.Uno /= Count;
493  m1.Due = sqrt(Sum2 - m1.Uno*NMax)/(double)(NMax-1);
494  return m1;
495 }
496 PARABOLA Matematica::MinimoParabola(double a, double b,double *Px,double *Py,int NMass){
497  PARABOLA Par;
498  double Uno=0.;double Due=0.;double Tre=0.;double Quattro=0.;
499  double UnoUno=0.;double DueUno=0.;double ZeroUno=0.;double ZeroDue=0.;
500  for(int i=0;*(Px+i)<b || i<NMass;i++){
501  // printf("%.0f\t",*(Px+i));
502  if(*(Px+i)>=a){
503  Uno += *(Px+i);
504  Due += QUAD(*(Px+i));
505  Tre += *(Px+i)*QUAD(*(Px+i));
506  Quattro += QUAD(QUAD(*Px+i));
507  UnoUno += *(Px+i) * *(Py+i);
508  DueUno += QUAD(*(Px+i)) * *(Py+i);
509  ZeroUno += *(Py+i);
510  ZeroDue += QUAD(*(Py+i));
511  // printf("%.2g\t%.2g\t%.2g\t%.2g\n",Uno,Due,Tre,Quattro);
512  }
513  }
514  double X2 = NMass*Due - Uno*Uno;
515  double X3 = NMass*Tre - Due*Uno;
516  double X4 = NMass*Quattro - Due*Due;
517  double XY = NMass *UnoUno - Uno*ZeroUno;
518  double X2Y = NMass*DueUno - Due*ZeroUno;
519  Par.a2 = (X4*X2 - X3*X3)/(X2Y*X2 - X3*XY);
520  Par.a1 = (Par.a2*XY - X3)/X2;
521  Par.a0 = (ZeroUno - Par.a1*Uno - Par.a2*Due)/(double)NMass;
522  Par.Minimo = -Par.a1/(2*Par.a2);
523  Par.MinimoY = Par.a2*SQR(Par.Minimo) + Par.a1*Par.Minimo + Par.a0;
524  //printf("y = %.2f*x^2 + %.2f*x + %.2f\t a=0.075,b=-50.93,c=8623\n",Par.a2,Par.a1,Par.a0);
525  return Par;
526 }
527 PARABOLA Matematica::MinimoParabola(double *Px,double *Py,int NMass){
528  PARABOLA Par;
529  double Uno=0.;double Due=0.;double Tre=0.;double Quattro=0.;
530  double UnoUno=0.;double DueUno=0.;double ZeroUno=0.;double ZeroDue=0.;
531  for(int i=0;i<NMass;i++){
532  Uno += Px[i];
533  Due += Px[i]*Px[i];
534  Tre += Px[i]*Px[i]*Px[i];
535  Quattro += Px[i]*Px[i]*Px[i]*Px[i];
536  UnoUno += Px[i]*Py[i];
537  DueUno += Px[i]*Px[i]*Py[i];
538  ZeroUno += Py[i];
539  ZeroDue += Py[i]*Py[i];
540  }
541  double X2 = NMass*Due - Uno*Uno;
542  double X3 = NMass*Tre - Due*Uno;
543  double X4 = NMass*Quattro - Due*Due;
544  double XY = NMass *UnoUno - Uno*ZeroUno;
545  double X2Y = NMass*DueUno - Due*ZeroUno;
546  Par.a2 = (X2Y*X2 - X3*XY)/(X4*X2 - X3*X3);
547  Par.a1 = (XY - Par.a2*X3)/X2;
548  Par.a0 = (ZeroUno - Par.a1*Uno - Par.a2*Due)/(double)NMass;
549  Par.Minimo = -Par.a1/(2.*Par.a2);
550  Par.MinimoY = Par.a2*SQR(Par.Minimo) + Par.a1*Par.Minimo + Par.a0;
551  //printf("y = %.2f*x^2 + %.2f*x + %.2f\t a=0.075,b=-50.93,c=8623\n",Par.a2,Par.a1,Par.a0);
552  return Par;
553 }
554 double Matematica::Blend(const double *dPoint,double x,int nPoint,int nOrder){
555  // if( ( x < dPoint[nPoint]) || ( x >= dPoint[nPoint+1]))
556  // return 0.;
557  if(nOrder == 1){
558  if( ( x >= dPoint[nPoint])&&( x < dPoint[nPoint+1]))
559  return 1.;
560  else
561  return 0.;
562  }
563  double Resp = 0.;
564  double Primo = (x - dPoint[nPoint])/(dPoint[nPoint+nOrder] - dPoint[nPoint])*Blend(dPoint,x,nPoint,nOrder-1);
565  double Secondo = (dPoint[nPoint+nOrder+1]-x)/(dPoint[nPoint+nOrder+1] - dPoint[nPoint+1])*Blend(dPoint,x,nPoint+1,nOrder-1);
566  if( (dPoint[nPoint+nOrder-1] == dPoint[nPoint]) && (dPoint[nPoint+nOrder] == dPoint[nPoint+1]) )
567  Resp = 0.;
568  else if(dPoint[nPoint+nOrder] == dPoint[nPoint])
569  Resp = Primo;
570  else if(dPoint[nPoint+nOrder] == dPoint[nPoint+1])
571  Resp = Secondo;
572  else
573  Resp = Primo + Secondo;
574  //printf("%d %d - (%lf-%lf) = %lf\n",nPoint,nOrder,x,dPoint[nPoint],Resp);
575  return Resp;
576 }
577 double Matematica::Blend(double *dPoint,size_t Incr,double x,int nP,int nO){
578 /*********************************************************************
579 
580  Simple b-spline curve algorithm
581 
582  Copyright 1994 by Keith Vertanen (vertankd@cda.mrs.umn.edu)
583 
584  Released to the public domain (your mileage may vary)
585 
586 **********************************************************************/
587  double Resp;
588  if(nO == 1){
589  if( (dPoint[nP*Incr] <= x)&&(x<dPoint[(nP+1)*Incr]))
590  Resp = 1.;
591  else
592  Resp = 0.;
593  }
594  else {
595  if( (dPoint[(nP+nO-1)*Incr] == dPoint[nP*Incr]) && (dPoint[(nP+nO)*Incr] == dPoint[(nP+1)*Incr]) )
596  Resp = 0.;
597  else if(dPoint[(nP+nO)*Incr] = dPoint[(nP)*Incr])
598  Resp = (dPoint[(nP+nO)*Incr] - x)/(dPoint[(nP+nO)*Incr] - dPoint[(nP+1)*Incr])*Blend(dPoint,Incr,x,nP,nO-1);
599  else if(dPoint[(nP+nO)*Incr] == dPoint[(nP+1)*Incr])
600  Resp = (x - dPoint[(nP)*Incr])/(dPoint[(nP+nO-1)*Incr] - dPoint[(nP)*Incr])*Blend(dPoint,Incr,x,nP+1,nO-1);
601  else
602  Resp = (dPoint[(nP+nO)*Incr] - x)/(dPoint[(nP+nO)*Incr] - dPoint[(nP+1)*Incr])*Blend(dPoint,Incr,x,nP,nO-1) + (x - dPoint[nP*Incr])/(dPoint[(nP+nO-1)*Incr] - dPoint[nP*Incr])*Blend(dPoint,Incr,x,nP+1,nO-1);
603  }
604  //printf("%d %d - (%lf-%lf) %d = %lf\n",nP,nO,x,dPoint[nP],Incr,Resp);
605  return Resp;
606 }
607 // int Matematica::InterBSpline2D(Matrice *MaIn,Matrice *MaOut){
608 // double Max=1.;
609 // double Min=0.;
610 // int NOutShow=0;
611 // double DeltaIn=(Max-Min)/(double)(NIn-1);
612 // double DeltaOut=(Max-Min)/(double)(NOut-1);
613 // int NOrder = 3+1;
614 // double *dArray = (double *)calloc(NIn+NOrder+1,sizeof(double));
615 // for(int p=0;p<=NIn+NOrder;p++){
616 // if(p<NOrder){
617 // dArray[p] = Min;
618 // }
619 // else if( (NOrder<=p) && (p<=NIn) ){
620 // dArray[p] = (p-NOrder+1)*DeltaIn+Min;//Pm[p-NOrder].Pos[CLat1];//
621 // }
622 // else if( p>NIn){
623 // dArray[p] = (p-NOrder+2)*DeltaIn+Min;
624 // }
625 // }
626 // for(int vo=0;vo<NOut;vo++){
627 // for(int vvo=0;vvo<NOut;vvo++){
628 // MaOut->Set(vo,vvo,0.);
629 // double x = DeltaOut*vo+Min;
630 // //for(int vi=0;vi<NIn;vi++){
631 // for(int vi=vo-1;vi<vo+NOrder+1;vi++){
632 // if(vi < 0 || vi >= NIn) continue;
633 // double Blendx = Mat->Blend(dArray,x,vi,NOrder);
634 // double y = DeltaOut*vvo+Min;
635 // //for(int vvi=0;vvi<NIn;vvi++){
636 // for(int vvi=vvo-1;vvi<vvo+NOrder+1;vvi++){
637 // if(vvi < 0 || vvi >= NIn) continue;
638 // double Blendy = Mat->Blend(dArray,y,vvi,NOrder);
639 // MaOut->Add(vo,vvo,Blendx*Blendy * PlIn->Val(vi,vvi));
640 // }
641 // }
642 // }
643 // }
644 // for(int vo=0;vo<NOut;vo++){
645 // if(vo < NIn){//To arrange
646 // PlOut[vo][NOut-1] = PlIn[vo][NIn-1];
647 // PlOut[NOut - 1][vo] = PlIn[NIn-1][vo];
648 // }
649 // }
650 // NOutShow = NOut;
651 // free(dArray);
652 // return NOutShow;
653 // }
654 //
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 a1
a0 + a1 x + a2 x^2
coeficients of a N order spline
int DerMatrix(double *Px, double *Py, int NMass, SPLINE Wg, Spline *Sp)
Boh.
double Min
Minimum value.
SPLINE Spline3Beg(double *P1, double *P2, double *P3, int x, int y)
Three order spline first boundary.
Parabolas coefficients.
SPLINE Spline4End(double *P1, double *P2, int x, int y)
Four order spline last boundary.
void AddCoe(double Val, int n)
Add to the n coefficient.
double xC
x position of the center
double Due
Second moment.
double Tre
Third moment.
PARABOLA MinimoParabola(double a, double b, double *Px, double *Py, int NMass)
Minimum of the Parabola between.
double a2
a0 + a1 x + a2 x^2
double LinInterp(double Px1, double Px2, double Py1, double Py2, double x)
Linear interpolation between two points.
SPLINE Spline4PreEnd(double *P1, double *P2, double *P3, int x, int y)
Four order spline just before the end.
double MinimoY
Minimum of the parabola.
double m
y = m*x + q
void Invert()
Invert.
double ErrQ
Error on the intercept.
Linear interpolation.
double Cov
Covariance.
SPLINE Parab2(double *PA, double *PB, double *PC, int x, int y)
Three points parabolic interpolation.
double a3
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
Moments of a distribution.
SPLINE Spline4(double *P1, double *P2, double *P3, double *P4, int x, int y)
Four order spline.
RETTA InterExp(double *Px, double *Py, int NMass)
Exponential interpolation.
double Rad
Radius of the circle.
double Val(int row)
Returns a value in 1d.
double r
r factor
double ErrM
Error on the slope.
double a2
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
double a0
a0 + a1 x + a2 x^2
double a1
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
double ErrY
Error a posteriori.
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.
Matrice computes the algebric operations on matrices.
double Elevato(double x, int Volte)
Integer power.
bool Set(int row, int column, double Val)
Set a coefficient.
CIRCLE Osculante(double *PA, double *PB, double *PC, int x, int y)
Osculant circle.
double Minimo
Minimum of the parabola.
double Max
Maximum value.
double Uno
First moment.
MOMENTI InterGauss(double *Px, double *Py, int NMass)
Gaussian interpolation.
SPLINE Parab(double *P1, double *P2, double *P3, int x, int y)
Three points parabolic interpolation.
int Num
Number of points considered.
SPLINE Spline3(double *P1, double *P2, double *P3, int x, int y)
Three order spline.
double a0
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
Radius and center of a circle.
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.
double QBezier(double *P1, double *P2, double *P3, double x, int y)
QBezier curve of three points.
double yC
y position of the center
double q
y = m*x + q
double Delta
Step interval.
void SetCoe(double Val, int n)
Set the n coefficient.
double Corr
Correlation.
RETTA InterRett(double *Px, double *Py, int NMass)
Linear interpolation.