Allink
v0.1
|
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) 00012 Resp = QUAD((1-x))*P1[y]+2*(1-x)*x*P2[y]+QUAD(x)*P3[y]; 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]; 00185 double Deltax2 = QUAD(P2[x]) - QUAD(P1[x]); 00186 double Dx = P3[x] - P1[x]; 00187 double Dx2 = QUAD(P3[x]) - QUAD(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]; 00228 double DxAB2 = QUAD(PA[x]) - QUAD(PB[x]); 00229 double DxAC2 = QUAD(PA[x]) - QUAD(PC[x]); 00230 double DyBA2 = QUAD(PB[y]) - QUAD(PA[y]); 00231 double DyCA2 = QUAD(PC[y]) - QUAD(PA[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; 00236 Cir.Rad = sqrt( QUAD((PA[x] - Cir.xC)) + QUAD((PA[y] - Cir.yC)) ); 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 00312 double Noto34 = DyBC*DxAD - DyAC*DxAD - DyDC*DxAB + DyAC*DxAB; 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++){ 00339 Sp->AddCoe( Coeff->Val(r,c)*Py[c] , 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++){ 00368 Sp->AddCoe( Coeff->Val(r,c)*Py[c] , r); 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]; 00386 Due += QUAD(Px[i]); 00387 ZeroUno += Py[i]; 00388 UnoUno += Px[i]*Py[i]; 00389 ZeroDue += QUAD(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]; 00426 Due += QUAD(Px[i]*Peso[i]); 00427 ZeroUno += Py[i]; 00428 UnoUno += Px[i]*Py[i]*Peso[i]; 00429 ZeroDue += QUAD(Py[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]; 00453 Due += QUAD(Px[i]); 00454 ZeroUno += log10(Py[i]); 00455 UnoUno += Px[i]*log10(Py[i]); 00456 ZeroDue += QUAD(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); 00504 Due += QUAD(*(Px+i)); 00505 Tre += *(Px+i)*QUAD(*(Px+i)); 00506 Quattro += QUAD(QUAD(*Px+i)); 00507 UnoUno += *(Px+i) * *(Py+i); 00508 DueUno += QUAD(*(Px+i)) * *(Py+i); 00509 ZeroUno += *(Py+i); 00510 ZeroDue += QUAD(*(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); 00639 // MaOut->Add(vo,vvo,Blendx*Blendy * PlIn->Val(vi,vvi)); 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 //