MatematicaInterp.cpp
```00001 #include "../include/Matematica.h"
00002
00003 double Matematica::QBezier(double *P1,double *P2,double *P3,double x,int y){
00004   double Resp;
00005 //   int Order = 3;
00006 //   double *Coeff;
00007 //   Coeff = (double *)malloc(Order*sizeof(double));
00008 //   for(int i=0;i<Order;i++){
00009 //     Coeff[i] = Binomial(Order,i)*P1[y];
00010 //   }//deCasteljau algorithm
00011 //   //P(x) = sum_i Coeff[i]*t^i*(t-1)^(n-1)
00013   return Resp;
00014 }
00015 SPLINE Matematica::Spline3Beg(double *P1,double *P2,double *P3,int x,int y){
00016   SPLINE Sp; Sp.a0 = 0.;
00017   if(x < 0 || x>3 || y < 0 || y>3){
00018     printf("x,y must be 0 <= %d , %d < 3\n",x,y);
00019     return Sp;
00020   }
00021   double Deltax = P2[x] - P1[x];
00022   Sp.a0 = P1[y];
00023   SPLINE Par = Parab2(P1,P2,P3,x,y);
00024   Sp.a2 = 0.;
00025   //  printf("%lf %lf\n",Sp.a2,Par.A);
00026   Sp.a3 = (Par.a2 - Sp.a2)/(3.*Deltax);
00027   Sp.a1 = (P2[y] - P1[y])/Deltax - Sp.a2*Deltax - Sp.a3*Deltax*Deltax;
00028   SpMem.a2 = Par.a2;
00029   SpMem.a1 = Sp.a1+2*Sp.a2*Deltax+3*Sp.a3*Deltax*Deltax;
00030   Sp.a4 = 0.;
00031   //printf("%lf %lf|%lf %lf %lf\n",Sp.a0,Sp.a1,SplineA1,Sp.a2,Sp.a3);
00032   return Sp;
00033 }
00034 SPLINE Matematica::Spline3(double *P1,double *P2,double *P3,int x,int y){
00035   SPLINE Sp; Sp.a0 = 0.;
00036   if(x < 0 || x>3 || y < 0 || y>3){
00037     printf("x,y must be 0 <= %d , %d < 3\n",x,y);
00038     return Sp;
00039   }
00040   double Deltax = P2[x] - P1[x];
00041   SPLINE Par = Parab2(P1,P2,P3,x,y);
00042   Sp.a0 = P1[y];
00043   Sp.a2 = SpMem.a2;
00044   Sp.a3 = (Par.a2 - Sp.a2)/(3.*Deltax);
00045   Sp.a1 = (P2[y] - P1[y])/Deltax - Sp.a2*Deltax - Sp.a3*Deltax*Deltax;
00046   SpMem.a2 = Par.a2;
00047   //printf("%lf %lf|%lf %lf|%lf %lf\n",Sp.a0,SplineA1,Sp.a1,Sp.a2,SplineA2,Sp.a3);
00048   SpMem.a1 = Sp.a1;
00049   Sp.a4 = 0.;
00050   return Sp;
00051 }
00052 SPLINE Matematica::Spline3End(double *P1,double *P2,int x,int y){
00053   SPLINE Sp; Sp.a0 = 0.;
00054   if(x < 0 || x>3 || y < 0 || y>3){
00055     printf("x,y must be 0 <= %d , %d < 3\n",x,y);
00056     return Sp;
00057   }
00058   double Deltax = P2[x] - P1[x];
00059   Sp.a0 = P1[y];
00060   Sp.a2 = SpMem.a2;
00061   Sp.a3 = (0. - SpMem.a2)/(3.*Deltax);
00062   Sp.a1 = (P2[y] - P1[y])/Deltax - Sp.a2*Deltax - Sp.a3*Deltax*Deltax;
00063   //  printf("%lf %lf %lf %lf\n",Sp.a0,Sp.a1,Sp.a2,Sp.a3);
00064   SpMem.a2 = 0.;
00065   Sp.a4 = 0.;
00066   return Sp;
00067 }
00068 SPLINE Matematica::Spline4Beg(double *P1,double *P2,double *P3,double *P4,int x,int y){
00069   SPLINE Sp; Sp.a0 = 0.;
00070   if(x < 0 || x>3 || y < 0 || y>3){
00071     printf("x,y must be 0 <= %d , %d < 3\n",x,y);
00072     return Sp;
00073   }
00074   double Deltax = P2[x] - P1[x];
00075   SPLINE Cub = Cubica(P1,P2,P3,P4,x,y);
00076   SPLINE Par = Parab2(P1,P2,P3,x,y);
00077   Sp.a0 = P1[y];
00078   Sp.a2 = 0.;
00079   Sp.a4 = 3.*Cub.a3/Deltax + (Sp.a2 - Cub.a2)/QUAD(Deltax);
00080   Sp.a4 /= 12.;
00081   Sp.a3 = Cub.a3 - 6.*Sp.a4*Deltax;
00082   Sp.a1 = (P2[y] - P1[y])/Deltax - Sp.a2*Deltax - Sp.a3*Deltax*Deltax - Sp.a4*Deltax*Deltax*Deltax;
00083   SpMem.a2 = Cub.a2;
00084   SpMem.a3 = Sp.a3;
00085   return Sp;
00086 }
00087 SPLINE Matematica::Spline4(double *P1,double *P2,double *P3,double *P4,int x,int y){
00088   SPLINE Sp; Sp.a0 = 0.;
00089   if(x < 0 || x>3 || y < 0 || y>3){
00090     printf("x,y must be 0 <= %d , %d < 3\n",x,y);
00091     return Sp;
00092   }
00093   double Deltax = P2[x] - P1[x];
00094   SPLINE Cub = Cubica(P1,P2,P3,P4,x,y);
00095   SPLINE Par = Parab2(P1,P2,P3,x,y);
00096   Sp.a0 = P1[y];
00097   Sp.a2 = SpMem.a2;
00098   Sp.a4 = 3.*Cub.a3/Deltax + (Sp.a2 - Cub.a2)/QUAD(Deltax);
00099   Sp.a4 /= 12.;
00100   Sp.a3 = Cub.a3 - 6.*Sp.a4*Deltax;
00101   Sp.a1 = (P2[y] - P1[y])/Deltax - Sp.a2*Deltax - Sp.a3*Deltax*Deltax - Sp.a4*Deltax*Deltax*Deltax;
00102   SpMem.a2 = Cub.a2;
00103   //  printf("%lf %lf|%lf %lf|%lf %lf %lf\n",Sp.a0,SpMem.a1,Sp.a1,Sp.a2,SpMem.a2,Sp.a3,Sp.a4);
00104   SpMem.a1 = Sp.a1;
00105   SpMem.a3 = Cub.a3;
00106   return Sp;
00107 }
00108 SPLINE Matematica::Spline4(double *P1,double *P2,double *P3,int x,int y){
00109   SPLINE Sp; Sp.a0 = 0.;
00110   if(x < 0 || x>3 || y < 0 || y>3){
00111     printf("x,y must be 0 <= %d , %d < 3\n",x,y);
00112     return Sp;
00113   }
00114   double Deltax = P2[x] - P1[x];
00115   SPLINE Par = Parab2(P1,P2,P3,x,y);
00116   Sp.a0 = P1[y];
00117   Sp.a2 = SpMem.a2;
00118   Sp.a1 = SpMem.a1;
00119 //   Sp.a1 = Par.B - 13.*Sp.a2*Deltax/2. - 3.*(P2[y] - P1[y])/Deltax - Par.A*Deltax/2.;
00120 //   Sp.a1 /= 13.;
00121   Sp.a4 = Par.a2 + 2.*Sp.a2 - 3.*(P2[y] - P1[y])/(Deltax*Deltax) + 3.*Sp.a1/Deltax;
00122   Sp.a4 /= 3.;
00123   Sp.a3 = (P2[y] - P1[y])/(Deltax*Deltax*Deltax) - Sp.a1/(Deltax*Deltax) - Sp.a2/Deltax - Sp.a4*Deltax;
00124   SpMem.a2 = Par.a2;
00125   //printf("%lf %lf|%lf %lf|%lf %lf %lf\n",Sp.a0,SpMem.a1,Sp.a1,Sp.a2,SpMem.a2,Sp.a3,Sp.a4);
00126   SpMem.a1 = Par.a1;
00127   return Sp;
00128 }
00129 SPLINE Matematica::Spline4PreEnd(double *P1,double *P2,double *P3,int x,int y){
00130   SPLINE Sp; Sp.a0 = 0.;
00131   if(x < 0 || x>3 || y < 0 || y>3){
00132     printf("x,y must be 0 <= %d , %d < 3\n",x,y);
00133     return Sp;
00134   }
00135   double Deltax = P2[x] - P1[x];
00136   SPLINE Par = Parab2(P1,P2,P3,x,y);
00137   Sp.a0 = P1[y];
00138   Sp.a2 = SpMem.a2;
00139   Sp.a3 = (Sp.a2 - Par.a2)/Deltax + 4.*SpMem.a3;
00140   Sp.a4 = (SpMem.a3 - Sp.a3)/(6.*Deltax);
00141   Sp.a1 = (P2[y] - P1[y])/Deltax - Sp.a2*Deltax - Sp.a3*Deltax*Deltax - Sp.a4*Deltax*Deltax*Deltax;
00142   SpMem.a2 = Par.a2;
00143   SpMem.a1 = Sp.a1;
00144   SpMem.a3 = Sp.a3;
00145   return Sp;
00146 }
00147 SPLINE Matematica::Spline4End(double *P1,double *P2,int x,int y){
00148   SPLINE Sp; Sp.a0 = 0.;
00149   if(x < 0 || x>3 || y < 0 || y>3){
00150     printf("x,y must be 0 <= %d , %d < 3\n",x,y);
00151     return Sp;
00152   }
00153   double Deltax = P2[x] - P1[x];
00154   Sp.a0 = P1[y];
00155   Sp.a2 = SpMem.a2;
00156   //Sp.a3 = (Sp.a2 - 0.)/Deltax + 4.*SpMem.a3;
00157   Sp.a3 = 0.;
00158   Sp.a4 = (SpMem.a3 - Sp.a3)/(6.*Deltax);
00159   Sp.a1 = (P2[y] - P1[y])/Deltax - Sp.a2*Deltax - Sp.a3*Deltax*Deltax - Sp.a4*Deltax*Deltax*Deltax;
00160   SpMem.a2 = 0.;
00161   SpMem.a1 = Sp.a1;
00162   SpMem.a3 = Sp.a3;
00163   return Sp;
00164 }
00165 //BSplines
00166 //   double w = x - floor(x) - 1.;
00167 //   xSp.a3 = (1./6.)*w * w * w;
00168 //   xSp.a0 = (1./6.) + (1./2.)*w*(w-1.)-xSp.a3;
00169 //   xSp.a2 = w + xSp.a0 - 2*xSp.a3;
00170 //   xSp.a1 = 1. - xSp.a0 - xSp.a2 - xSp.a3;
00171 //   Resp = xSp.a0*P1[0]+xSp.a1*P2[0]+xSp.a2*P3[0];
00172   //  w = y - 3.;
00173 //   ySp.a3 = (1./6.)*w * w * w;
00174 //   ySp.a0 = (1./6.) + (1./2.)*w*(w-1.)-ySp.a3;
00175 //   ySp.a2 = w + ySp.a0 - 2*ySp.a3;
00176 //   ySp.a1 = 1. - ySp.a0 - ySp.a2 - ySp.a3;
00177   //  Resp *= ySp.a0*P1[1]+ySp.a1*P2[1]+ySp.a2*P3[1];
00178 SPLINE Matematica::Parab(double *P1,double *P2,double *P3,int x,int y){
00179   SPLINE Par; Par.a0 = 0.;
00180   if(x < 0 || x>3 || y < 0 || y>3){
00181     printf("x,y must be 0 <= %d , %d < 3\n",x,y);
00182     return Par;
00183   }
00184   double Deltax = P2[x] - P1[x];
00186   double Dx = P3[x] - P1[x];
00188   double Deltay = P2[y] - P1[y];
00189   double Dy = P3[y] - P1[y];
00190   double b1= Dy - Deltay/Deltax2*Dx2;
00191   double b2= Dx - Deltax/Deltax2*Dx2;
00192   Par.a1 = b1/b2;
00193   Par.a2 = Deltay/Deltax2 - Par.a1*Deltax/Deltax2;
00194   Par.a0 = P1[y] - Par.a1*P1[x] - Par.a2*QUAD(P1[x]);
00195   Par.a3 = 0.;
00196   Par.a4 = 0.;
00197   return Par;
00198 }
00199 SPLINE Matematica::Parab2(double *PA,double *PB,double *PC,int x,int y){
00200   SPLINE Par; Par.a0 = 0.;
00201   if(x < 0 || x>3 || y < 0 || y>3){
00202     printf("x,y must be 0 <= %d , %d < 3\n",x,y);
00203     return Par;
00204   }
00205   double DxAB = PA[x] - PB[x];
00206   double DxCB = PC[x] - PB[x];
00207   double DyCB = PC[y] - PB[y];
00208   double DyAB = PA[y] - PB[y];
00209   double a1 = DyCB * DxAB - DyAB * DxCB;
00210   double a2 = DxCB*DxCB*DxAB - DxCB*DxAB*DxAB;
00211   Par.a2 = a1/a2;
00212   Par.a1 = DyAB/DxAB - Par.a2*DxAB;
00213   Par.a0 = PB[y];
00214   Par.a3 = 0.;
00215   Par.a4 = 0.;
00216   return Par;
00217 }
00218 CIRCLE Matematica::Osculante(double *PA,double *PB,double *PC,int x,int y){
00219   CIRCLE Cir; Cir.yC = 0.;
00220   if(x < 0 || x>3 || y < 0 || y>3){
00221     printf("x,y must be 0 <= %d , %d < 3\n",x,y);
00222     return Cir;
00223   }
00224   double DxCA = PC[x] - PA[x];
00225   double DxBA = PB[x] - PA[x];
00226   double DyAC = PA[y] - PC[y];
00227   double DyAB = PA[y] - PB[y];
00232   double a1 = DxCA*(DyBA2 - DxAB2) - DxBA*(DyCA2 - DxAB2);
00233   double a2 = 2.*(DyAC*DxBA - DyAB*DxCA);
00234   Cir.yC = a1/a2;
00235   Cir.xC = (DyBA2 - DxAB2 + 2.*Cir.yC*DyAB)/DxBA;
00237   return Cir;
00238 }
00239 SPLINE Matematica::Cubica(double *PA,double *PB,double *PC,double *PD,int x,int y){
00240   SPLINE Cub; Cub.a0 = 0.;
00241   if(x < 0 || x>3 || y < 0 || y>3){
00242     printf("x,y must be 0 <= %d , %d < 3\n",x,y);
00243     return Cub;
00244   }
00245   double DyDB = PD[y] - PB[y];
00246   double DyCB = PC[y] - PB[y];
00247   double DyAB = PA[y] - PB[y];
00248   double DxDB = PD[x] - PB[x];
00249   double DxCB = PC[x] - PB[x];
00250   double DxAB = PA[x] - PB[x];
00251   double DxCA = PC[x] - PA[x];
00252   double Num3 = DyDB/DxDB - DyAB/DxAB - (DyCB/DxCB - DyAB/DxAB)*(DxDB - DxAB)/DxCA;
00253   double Den3 = DxDB*DxDB - DxAB*DxAB - (DxCB*DxCB - DxAB*DxAB)*(DxDB - DxAB)/DxCA;
00254   Cub.a3 = Num3/Den3;
00255   Cub.a2 = DyCB/DxCB - DyAB/DxAB - Cub.a3*(DxCB*DxCB - DxAB*DxAB);
00256   Cub.a2 /= DxCA;
00257   Cub.a1 = DyAB/DxAB - Cub.a2*DxAB - Cub.a3*DxAB*DxAB;
00258   Cub.a0 = PB[y];
00259   Cub.a4 = 0.;
00260   return Cub;
00261 }
00262 SPLINE Matematica::Forth(double *PA,double *PB,double *PC,double *PD,double *PE,int x,int y){
00263   SPLINE Forth; Forth.a0=0.;
00264   if(x < 0 || x>3 || y < 0 || y>3){
00265     printf("x,y must be 0 <= %d , %d < 3\n",x,y);
00266     return Forth;
00267   }
00268   double Dx = ASS((PC[x] - PB[x]));
00269   Forth.a3 = PA[y]/12. - PB[y]/6. + PD[y]/6. - PE[y]/12.;
00270   Forth.a3 /= -Dx*Dx*Dx;
00271   Forth.a1 = - PA[y]/12. + 4.*PB[y]/6. - 4.*PD[y]/6. + PE[y]/12.;
00272   Forth.a1 /= -Dx;
00273   Forth.a4 = PA[y]/24. - PB[y]/6. + 1./4.*PC[y] - PD[y]/6. + PE[y]/24.;
00274   Forth.a4 /= Dx*Dx*Dx*Dx;
00275   Forth.a2 = -PA[y]/24. + 2./3.*PB[y] - 5./4.*PC[y] + 2./3.*PD[y] - PE[y]/24.;
00276   Forth.a2 /= Dx*Dx;
00277   Forth.a0 = PC[y];
00278   //printf("Dx %lf Coeff %lf %lf %lf %lf %lf \n",Dx,Forth.a0,Forth.a1,Forth.a2,Forth.a3,Forth.a4);
00279   return Forth;
00280
00281   double DxCA = PC[x] - PA[x];
00282   double DxCB = PC[x] - PB[x];
00283   double DxCD = PC[x] - PD[x];
00284   double DxCE = PC[x] - PE[x];
00285   double DxAB = PA[x] - PB[x];
00286   double DxAD = PA[x] - PD[x];
00287   double DyAC = (PA[y] - PC[y])/DxCA;
00288   double DyBC = (PB[y] - PC[y])/DxCB;
00289   double DyDC = (PD[y] - PC[y])/DxCD;
00290   double DyEC = (PE[y] - PC[y])/DxCE;
00291   //if(DxCB == DxCD)
00292   if (1==1)
00293     {
00294       double DxDB = PD[x] - PB[x];
00295       double DxEA = PE[x] - PA[x];
00296       double DyBD = (PB[y] - PD[y])/DxDB;
00297       double DyAE = (PA[y] - PE[y])/DxEA;
00298       //      Forth.a4 = - Noto34/Num34;
00299       Forth.a4 = DyBC - DyDC - DyEC + DyAC;
00300       Forth.a4 /= (DxCE*DxCE*DxCE - DxCA*DxCA*DxCA) - (DxCB*DxCB*DxCB - DxCD*DxCD*DxCD);
00301       Forth.a2 = DyBC - DyDC - Forth.a4*(DxCB*DxCB*DxCB - DxCD*DxCD*DxCD);
00302       Forth.a3 = DyBD - DyAE;
00303       Forth.a3 /= (DxCB*DxCB*DxCB - DxCD*DxCD*DxCD)/DxDB -
00304    (DxCA*DxCA*DxCA - DxCE*DxCE*DxCE)/DxEA;
00305       Forth.a1 = DyAE - Forth.a3*(DxCA*DxCA*DxCA - DxCE*DxCE*DxCE)/DxEA;
00306     }
00307   else
00308     {
00309       double DxAE = PA[x] - PE[x];
00310       double Den34 = DxAB*(DxCB*DxCB - DxCA*DxCA) - DxAB*(DxCD*DxCD - DxCA*DxCA);// 0.;
00311       double Num34 = DxAB*(DxCD*DxCD*DxCD - DxCA*DxCA*DxCA) - DxAD*(DxCB*DxCB*DxCB - DxCA*DxCA*DxCA);//21Dx
00313       double Num4 = -DxAB*(Noto34*(DxCE*DxCE - DxCA*DxCA) - DyEC + DyAC) +
00314    DxAE*(Noto34*(DxCB*DxCB - DxCA*DxCA) - DyBC + DyAC);
00315       double Den4 = DxAB*(Num34/Den34*(DxCE*DxCE - DxCA*DxCA) + (DxCE*DxCE*DxCE -  DxCA*DxCA*DxCA) ) -
00316    DxAE*(Num34/Den34*(DxCB*DxCB - DxCA*DxCA) + (DxCB*DxCB*DxCB - DxCA*DxCA*DxCA));
00317       //printf("Dx %lf %lf %lf %lf %lf %lf %lf\n",DxCB,Den34,Num34,Noto34,Num4,Den4,Num34/Den34);
00318       Forth.a4 = Num4/Den4;
00319       Forth.a3 = Num34/Den34*Forth.a4 + Noto34;
00320       Forth.a2 = DyBC - DyAC - Forth.a3*(DxCB*DxCB - DxCA*DxCA) - Forth.a4* (DxCB*DxCB*DxCB - DxCA*DxCA*DxCA);
00321       Forth.a2 /= DxAB;
00322       Forth.a1 = DyAC - Forth.a2*DxCA - Forth.a3*DxCA*DxCA - Forth.a4*DxCA*DxCA*DxCA;
00323     }
00324   Forth.a0 = PC[y];
00325   //printf("%lf %lf %lf %lf %lf \n",Forth.a0,Forth.a1,Forth.a2,Forth.a3,Forth.a4);
00326   return Forth;
00327 }
00328 int Matematica::Polinomio(double *Px,double *Py,int NMass,Spline *Sp){
00329   Matrice *Coeff = new Matrice(NMass,NMass);
00330   for(int r=0;r<NMass;r++){
00331     Sp->SetCoe( 0. , r);
00332     for(int c=0;c<NMass;c++){
00333       Coeff->Set(r,c,Elevato(Px[r],c));
00334     }
00335   }
00336   Coeff->Invert();
00337   for(int c=0;c<NMass;c++){
00338     for(int r=0;r<NMass;r++){
00340     }
00341   }
00342   delete  Coeff;
00343   return 1;
00344 }
00345 int Matematica::DerMatrix(double *Px,double *Py,int NMass,SPLINE Wg,Spline *Sp){
00346   Matrice *Coeff = new Matrice(NMass,NMass);
00347   double MenoDue = - .75*Wg.a3 + 1.5*Wg.a4;
00348   //MenoDue += .125*Wg.a1 - .125*Wg.a2; //O(h^4)
00349   double MenoUno = -.5*Wg.a1 + Wg.a2 + 1.5*Wg.a3 - 6.*Wg.a4;
00350   //MenoUno += -.125*Wg.a1 - .5*Wg.a2;
00351   double Zero = Wg.a0 - 2.*Wg.a2 + 9.*Wg.a4;
00352   //Zero += 0.75*Wg.a2;
00353   double PiuUno = .5*Wg.a1 + Wg.a2 - 1.5*Wg.a3 - 6.*Wg.a4;
00354   //PiuUno += .25*Wg.a1 + .5*Wg.a2;
00355   double PiuDue = .75*Wg.a3 + 1.5*Wg.a4;
00356   //PiuDue += -.125*Wg.a1 -.125*Wg.a2;
00357   for(int r=0;r<NMass;r++){
00358     if(r> 1) Coeff->Set(r,r-2,MenoDue*Px[r-2]);
00359     if(r> 0) Coeff->Set(r,r-1,MenoUno*Px[r-1]);
00360     if(r< NMass-1) Coeff->Set(r,r+1,PiuUno*Px[r+1]);
00361     if(r< NMass-2) Coeff->Set(r,r+2,PiuDue*Px[r+2]);
00362     Coeff->Set(r,r,Zero*Px[r]);
00363   }
00364   //Coeff->Print();
00365   Coeff->Invert();
00366   for(int r=0;r<NMass;r++){
00367     for(int c=0;c<NMass;c++){
00369     }
00370   }
00371   delete Coeff;
00372   return 0;
00373 }
00374 double Matematica::LinInterp(double Px1,double Px2,double Py1,double Py2,double x){
00375   double m = (Py2-Py1)/(Px2-Px1);
00376   double q = Py1 - m*Px1;
00377   return m*x+q;
00378 }
00379 RETTA Matematica::InterRett(double *Px,double *Py,int NMass){
00380   RETTA r1;
00381   double Uno=0.;double Due=0.;double Tre=0.;double Quattro=0.;
00382   double UnoUno=0.;double DueUno=0.;double ZeroUno=0.;double ZeroDue=0.;
00383   for(int i=0;i<NMass;i++){
00384     //printf("%lf %lf\n",Px[i],Py[i]);
00385     Uno += Px[i];
00387     ZeroUno += Py[i];
00388     UnoUno += Px[i]*Py[i];
00390   }
00391   double Mediax = Uno / (double) NMass;
00392   double Mediay = ZeroUno / (double) NMass;
00393   double Scartox = (Due - NMass*Uno*Uno)/(double)(NMass-0);
00394   double Scartoy = (ZeroDue - NMass*ZeroUno*ZeroUno)/(double)(NMass-0);
00395   r1.m = (NMass*UnoUno - Uno*ZeroUno) / (NMass*Due - Uno*Uno);
00396   r1.q = (ZeroUno - r1.m*Uno)/NMass;
00397   double Posteriori = 0.;
00398   r1.Cov = UnoUno/(double)NMass - Mediax*Mediay;
00399   for(int i=0;i<NMass;i++){
00400     Posteriori += QUAD(( Px[i]*r1.m + r1.q - Py[i] ));
00401   }
00402   r1.Corr = r1.Cov / ( sqrt(Scartox*Scartoy) );
00403   r1.ErrY = sqrt(Posteriori/(double)(NMass-2));
00404   r1.ErrM = r1.ErrY * sqrt( NMass / (NMass * Due - Uno*Uno));
00405   r1.ErrQ = r1.ErrY * sqrt( Due / (NMass*Due - Uno*Uno) );
00406   //  printf("m %lf q %lf r %lf sigma %lf Mediax %lf Mediay %lf\n",r1.m,r1.q,r1.r,r1.Corr,Mediax,Mediay);
00407   return r1;
00408 }
00409 // Vettore Matematica::Directive(double *Px,int NMass){
00410 //   Vettore v1(3);
00411
00412
00413 // }
00414 RETTA Matematica::InterRett(double *Px,double *Py,double *Peso,int NMass){
00415   RETTA r1;
00416   double Mediax = 0.;
00417   double Mediay = 0.;
00418   double Scartox = 0.;
00419   double Scartoy = 0.;
00420   double Pesox = 0.;
00421   double Uno=0.;double Due=0.;double Tre=0.;double Quattro=0.;
00422   double UnoUno=0.;double DueUno=0.;double ZeroUno=0.;double ZeroDue=0.;
00423   for(int i=0;i<NMass;i++){
00424     Pesox += Peso[i];
00425     Uno += Px[i]*Peso[i];
00427     ZeroUno += Py[i];
00428     UnoUno += Px[i]*Py[i]*Peso[i];
00430   }
00431   Mediax = Uno / (double) NMass/Pesox;
00432   Mediay = ZeroUno / (double) NMass;
00433   Scartox = (Due - Uno*Uno/QUAD(Pesox))/(double)NMass;
00434   Scartoy = (ZeroDue - ZeroUno*ZeroUno)/(double)NMass;
00435   r1.Corr = (UnoUno - Uno*ZeroUno) / (double) NMass/Pesox ;
00436   r1.r = r1.Corr / ( sqrt(Scartox*Scartoy) );
00437   r1.m = (NMass*UnoUno - Uno*ZeroUno) / (NMass*Due - Uno*Uno)/Pesox;
00438   r1.ErrM = Scartoy * sqrt( Due / (NMass * Due - Uno*Uno)/QUAD(Pesox));
00439   r1.q = (ZeroUno - r1.m*Uno)/NMass;
00440   r1.ErrQ = Scartoy * sqrt( NMass / (NMass * Due - Uno*Uno)/QUAD(Pesox) );
00441   return r1;
00442 }
00443 RETTA Matematica::InterExp(double *Px,double *Py,int NMass){
00444   RETTA r1;
00445   double Uno=0.;double Due=0.;double Tre=0.;double Quattro=0.;
00446   double UnoUno=0.;double DueUno=0.;double ZeroUno=0.;double ZeroDue=0.;
00447   for(int i=0;i<NMass;i++){
00448     if(Py[i] <= 0.){
00449       //printf("Negative number not allowed for an exponential interpolation %lf\n",Py[i]);
00450       continue;
00451     }
00452     Uno += Px[i];
00454     ZeroUno += log10(Py[i]);
00455     UnoUno += Px[i]*log10(Py[i]);
00457   }
00458   double Mediax = Uno / (double) NMass;
00459   double Mediay = ZeroUno / (double) NMass;
00460   double Scartox = (Due - NMass*Uno*Uno)/(double)(NMass-1);
00461   double Scartoy = (ZeroDue - NMass*ZeroUno*ZeroUno)/(double)(NMass-1);
00462   r1.m = (NMass*UnoUno - Uno*ZeroUno) / (NMass*Due - Uno*Uno);
00463   r1.q = (ZeroUno - r1.m*Uno)/NMass;
00464   double Posteriori = 0.;
00465   r1.Cov = UnoUno/(double)NMass - Mediax*Mediay;
00466   for(int i=0;i<NMass;i++){
00467     Posteriori += QUAD(( Px[i]*r1.m + r1.q - log10(Py[i]) ));
00468   }
00469   r1.Corr = r1.Cov / ( sqrt(Scartox*Scartoy) );
00470   r1.ErrY = sqrt(Posteriori/(double)(NMass-2));
00471   r1.ErrM = r1.ErrY * sqrt( NMass / (NMass * Due - Uno*Uno));
00472   r1.ErrQ = r1.ErrY * sqrt( Due / (NMass*Due - Uno*Uno) );
00473   //  printf("m %lf q %lf r %lf sigma %lf Mediax %lf Mediay %lf\n",r1.m,r1.q,r1.r,r1.Corr,Mediax,Mediay);
00474   return r1;
00475 }
00476 MOMENTI Matematica::InterGauss(double *Px,double *Py,int NMax){
00477   MOMENTI m1; m1.Uno=0.; m1.Due=0.; m1.Tre=0.;m1.Delta=0.;m1.Num=0;
00478   m1.Min = Px[0];
00479   m1.Max = Px[0];
00480   double Count = 0.;
00481   double Sum2 = 0.;
00482   for(int i=0;i<NMax;i++){
00483 #ifdef MAT_DEBUG
00484     if(Py[i] < 0.){printf("Invalid weight %lf < 0\n",Py[i]);}
00485 #endif
00486     if(Px[i]<m1.Min) m1.Min = Px[i];
00487     if(Px[i]>m1.Max) m1.Max = Px[i];
00488     m1.Uno += Px[i]*Py[i];
00489     Sum2 += SQR(Px[i])*Py[i];
00490     Count += Py[i];
00491   }
00492   m1.Uno /= Count;
00493   m1.Due = sqrt(Sum2 - m1.Uno*NMax)/(double)(NMax-1);
00494   return m1;
00495 }
00496 PARABOLA Matematica::MinimoParabola(double a, double b,double *Px,double *Py,int NMass){
00497   PARABOLA Par;
00498   double Uno=0.;double Due=0.;double Tre=0.;double Quattro=0.;
00499   double UnoUno=0.;double DueUno=0.;double ZeroUno=0.;double ZeroDue=0.;
00500   for(int i=0;*(Px+i)<b || i<NMass;i++){
00501     //    printf("%.0f\t",*(Px+i));
00502     if(*(Px+i)>=a){
00503       Uno += *(Px+i);
00507       UnoUno += *(Px+i) * *(Py+i);
00508       DueUno += QUAD(*(Px+i)) * *(Py+i);
00509       ZeroUno += *(Py+i);
00511       //      printf("%.2g\t%.2g\t%.2g\t%.2g\n",Uno,Due,Tre,Quattro);
00512     }
00513   }
00514   double X2 = NMass*Due - Uno*Uno;
00515   double X3 = NMass*Tre - Due*Uno;
00516   double X4 = NMass*Quattro - Due*Due;
00517   double XY = NMass *UnoUno - Uno*ZeroUno;
00518   double X2Y = NMass*DueUno - Due*ZeroUno;
00519   Par.a2 = (X4*X2 - X3*X3)/(X2Y*X2 - X3*XY);
00520   Par.a1 = (Par.a2*XY - X3)/X2;
00521   Par.a0 = (ZeroUno - Par.a1*Uno - Par.a2*Due)/(double)NMass;
00522   Par.Minimo = -Par.a1/(2*Par.a2);
00523   Par.MinimoY = Par.a2*SQR(Par.Minimo) + Par.a1*Par.Minimo + Par.a0;
00524   //printf("y = %.2f*x^2 + %.2f*x + %.2f\t a=0.075,b=-50.93,c=8623\n",Par.a2,Par.a1,Par.a0);
00525   return Par;
00526 }
00527 PARABOLA Matematica::MinimoParabola(double *Px,double *Py,int NMass){
00528   PARABOLA Par;
00529   double Uno=0.;double Due=0.;double Tre=0.;double Quattro=0.;
00530   double UnoUno=0.;double DueUno=0.;double ZeroUno=0.;double ZeroDue=0.;
00531   for(int i=0;i<NMass;i++){
00532     Uno += Px[i];
00533     Due += Px[i]*Px[i];
00534     Tre += Px[i]*Px[i]*Px[i];
00535     Quattro += Px[i]*Px[i]*Px[i]*Px[i];
00536     UnoUno += Px[i]*Py[i];
00537     DueUno += Px[i]*Px[i]*Py[i];
00538     ZeroUno += Py[i];
00539     ZeroDue += Py[i]*Py[i];
00540   }
00541   double X2 = NMass*Due - Uno*Uno;
00542   double X3 = NMass*Tre - Due*Uno;
00543   double X4 = NMass*Quattro - Due*Due;
00544   double XY = NMass *UnoUno - Uno*ZeroUno;
00545   double X2Y = NMass*DueUno - Due*ZeroUno;
00546   Par.a2 = (X2Y*X2 - X3*XY)/(X4*X2 - X3*X3);
00547   Par.a1 = (XY - Par.a2*X3)/X2;
00548   Par.a0 = (ZeroUno - Par.a1*Uno - Par.a2*Due)/(double)NMass;
00549   Par.Minimo = -Par.a1/(2.*Par.a2);
00550   Par.MinimoY = Par.a2*SQR(Par.Minimo) + Par.a1*Par.Minimo + Par.a0;
00551   //printf("y = %.2f*x^2 + %.2f*x + %.2f\t a=0.075,b=-50.93,c=8623\n",Par.a2,Par.a1,Par.a0);
00552   return Par;
00553 }
00554 double Matematica::Blend(const double *dPoint,double x,int nPoint,int nOrder){
00555   // if( ( x < dPoint[nPoint]) || ( x >= dPoint[nPoint+1]))
00556   //   return 0.;
00557   if(nOrder == 1){
00558     if( ( x >= dPoint[nPoint])&&( x < dPoint[nPoint+1]))
00559       return 1.;
00560     else
00561       return 0.;
00562   }
00563   double Resp = 0.;
00564   double Primo = (x - dPoint[nPoint])/(dPoint[nPoint+nOrder] - dPoint[nPoint])*Blend(dPoint,x,nPoint,nOrder-1);
00565   double Secondo = (dPoint[nPoint+nOrder+1]-x)/(dPoint[nPoint+nOrder+1] - dPoint[nPoint+1])*Blend(dPoint,x,nPoint+1,nOrder-1);
00566   if( (dPoint[nPoint+nOrder-1] == dPoint[nPoint]) && (dPoint[nPoint+nOrder] == dPoint[nPoint+1]) )
00567     Resp = 0.;
00568   else if(dPoint[nPoint+nOrder] == dPoint[nPoint])
00569     Resp = Primo;
00570   else if(dPoint[nPoint+nOrder] == dPoint[nPoint+1])
00571     Resp = Secondo;
00572   else
00573     Resp = Primo + Secondo;
00574   //printf("%d %d - (%lf-%lf) = %lf\n",nPoint,nOrder,x,dPoint[nPoint],Resp);
00575   return Resp;
00576 }
00577 double Matematica::Blend(double *dPoint,size_t Incr,double x,int nP,int nO){
00578 /*********************************************************************
00579
00580  Simple b-spline curve algorithm
00581
00582  Copyright 1994 by Keith Vertanen (vertankd@cda.mrs.umn.edu)
00583
00584  Released to the public domain (your mileage may vary)
00585
00586 **********************************************************************/
00587   double Resp;
00588   if(nO == 1){
00589     if( (dPoint[nP*Incr] <= x)&&(x<dPoint[(nP+1)*Incr]))
00590       Resp = 1.;
00591     else
00592       Resp = 0.;
00593   }
00594   else {
00595     if( (dPoint[(nP+nO-1)*Incr] == dPoint[nP*Incr]) && (dPoint[(nP+nO)*Incr] == dPoint[(nP+1)*Incr]) )
00596       Resp = 0.;
00597     else if(dPoint[(nP+nO)*Incr] = dPoint[(nP)*Incr])
00598       Resp = (dPoint[(nP+nO)*Incr] - x)/(dPoint[(nP+nO)*Incr] - dPoint[(nP+1)*Incr])*Blend(dPoint,Incr,x,nP,nO-1);
00599     else if(dPoint[(nP+nO)*Incr] == dPoint[(nP+1)*Incr])
00600       Resp = (x - dPoint[(nP)*Incr])/(dPoint[(nP+nO-1)*Incr] - dPoint[(nP)*Incr])*Blend(dPoint,Incr,x,nP+1,nO-1);
00601     else
00602       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);
00603   }
00604   //printf("%d %d - (%lf-%lf) %d = %lf\n",nP,nO,x,dPoint[nP],Incr,Resp);
00605   return Resp;
00606 }
00607 // int Matematica::InterBSpline2D(Matrice *MaIn,Matrice *MaOut){
00608 //   double Max=1.;
00609 //   double Min=0.;
00610 //   int NOutShow=0;
00611 //   double DeltaIn=(Max-Min)/(double)(NIn-1);
00612 //   double DeltaOut=(Max-Min)/(double)(NOut-1);
00613 //   int NOrder = 3+1;
00614 //   double *dArray = (double *)calloc(NIn+NOrder+1,sizeof(double));
00615 //   for(int p=0;p<=NIn+NOrder;p++){
00616 //     if(p<NOrder){
00617 //       dArray[p] = Min;
00618 //     }
00619 //     else if( (NOrder<=p) && (p<=NIn) ){
00620 //       dArray[p] = (p-NOrder+1)*DeltaIn+Min;//Pm[p-NOrder].Pos[CLat1];//
00621 //     }
00622 //     else if( p>NIn){
00623 //       dArray[p] = (p-NOrder+2)*DeltaIn+Min;
00624 //     }
00625 //   }
00626 //   for(int vo=0;vo<NOut;vo++){
00627 //    for(int vvo=0;vvo<NOut;vvo++){
00628 //      MaOut->Set(vo,vvo,0.);
00629 //      double x = DeltaOut*vo+Min;
00630 //      //for(int vi=0;vi<NIn;vi++){
00631 //      for(int vi=vo-1;vi<vo+NOrder+1;vi++){
00632 //        if(vi < 0 || vi >= NIn) continue;
00633 //        double Blendx = Mat->Blend(dArray,x,vi,NOrder);
00634 //        double y = DeltaOut*vvo+Min;
00635 //        //for(int vvi=0;vvi<NIn;vvi++){
00636 //        for(int vvi=vvo-1;vvi<vvo+NOrder+1;vvi++){
00637 //     if(vvi < 0 || vvi >= NIn) continue;
00638 //     double Blendy = Mat->Blend(dArray,y,vvi,NOrder);
00640 //        }
00641 //      }
00642 //    }
00643 //   }
00644 //   for(int vo=0;vo<NOut;vo++){
00645 //     if(vo < NIn){//To arrange
00646 //       PlOut[vo][NOut-1] = PlIn[vo][NIn-1];
00647 //       PlOut[NOut - 1][vo] = PlIn[NIn-1][vo];
00648 //     }
00649 //   }
00650 //   NOutShow = NOut;
00651 //   free(dArray);
00652 //   return NOutShow;
00653 // }
00654 //
```