Allink: src/Matematica/MatematicaFunc.cpp Source File

00001 #include "../include/Matematica.h"
00002 
00003 double Matematica::ContactAngle(double x){
00004   double Ris;
00005   //Ris = PreFact*pow(3./(2.+3.*sin(x)-CUB(sin(x)) ), 4./3. )*QUAD(QUAD(cos(x))) - Ypsilon;
00006   Ris = PreFact*pow(1-cos(x), -4./3. )*pow(sin(x),2.) - Ypsilon;
00007   return Ris;
00008 }
00009 double Matematica::fProva(double x){
00010   double Ris;
00011   //Ris = pow(x,5.)/((exp(x)-1)*(1-exp(-x)));
00012   // Ris = pow(x,5.)*exp(x)/QUAD(exp(x)-1);
00013   // Ris = 1.*pow(x-150,4.)-1.*pow(x-150,5.)+x;
00014   // Ris = pow(x-150,2.)+x-150;
00015   Ris = 1./(13.9332529115775*x*x - 4.*13.9332529115775*0.13710248135043*x + 7.86306085419497) + .5/(-0.0805855688278211*x + 10.6925178203198);
00016   //Ris = 1./(12.5342*x*x) + .5/(-0.207807999999993*x + 17.4283);
00017   return Ris;
00018 }
00019 double Matematica::F(double T,double TD){
00020   return pow((T/TD),5.)*Integrazione(.00001,TD/T);
00021 }
00022 double Matematica::Integrazione(double a,double b){
00023   double Delta = (b-a)/(double)NPassi;
00024   double Risp = 0.;
00025   for(int i=0;a+((double)i)*Delta<b;i+=3){//+3 se si sovrappongono
00026     Risp += 3.*3.*Delta*( Evalx(a+Delta*i)+3.*Evalx(a+Delta*(i+1))+3.*Evalx(a+Delta*(i+2))+Evalx(a+Delta*(i+3)) )/8.;
00027   }
00028   return Risp;
00029 }
00030 double Matematica::Integrazione(double *st,double *sw,int NMass){
00031   double NMassInv = 1/((double)NMass);
00032   sw[0] = st[0];
00033   for(int i=1;i<NMass;i++){
00034     sw[i] = st[i] + sw[i-1];
00035   }
00036   return sw[NMass-1];
00037 }
00038 double Matematica::Df(double x,double Delta){
00039   return ( Evalx(x+Delta)-Evalx(x) )/Delta;
00040 }
00041 void Matematica::Derivata(double *st,double *sw,int NMass){
00042   for(int i=0;i<NMass-1;i++){
00043     sw[i] = (st[i+1]-st[i])/2.;
00044   }
00045 }
00046 void Matematica::DerO4(double *st,double *sw,int NMass){
00047   sw[0] = 0.;
00048   sw[NMass-1] = 0.;
00049   for(int i=1;i<NMass-1;i++){
00050     if(i<2)
00051       sw[i] = (st[i+1] - st[i-1])/2.;
00052     else if (i < NMass - 2)
00053       sw[i] = st[i-2] - 8.*st[i-1] + 8.*st[i+1] - st[i+2];
00054     else if (i < NMass - 1)
00055       sw[i] = (st[i+1] - st[i-1])/2.;
00056   }
00057 }
00058 int Matematica::Zeri(double a,double b,double *Radici,int NRadici){
00059   double Uno;double Due;double Delta;
00060   int rad=0;
00061   RADICE Rad;
00062   if( a > b ){
00063     Uno = b ; Due = a;
00064   }
00065   else{
00066     Uno = a ; Due = b;
00067   }
00068   Delta = (Due - Uno)/2.;
00069   //  for(int p=0;p<4;p++){
00070   for(int i=0;i<NRadici*2;i++){
00071       Delta = (Due - Uno)/(2*NRadici);
00072       Rad = RegulaFalsi(Uno+i*Delta,Due-(2*NRadici-i-1)*Delta);
00073       //printf("%d %d) %g|%g %g|%g %g|%g\n",NRadici,i,Uno+i*Delta,Rad.iLim,Due-(2*NRadici-i-1)*Delta,Rad.sLim,Delta,Rad.sLim-Rad.iLim);
00074       if(Rad.IfRis == 1){
00075    printf("Found a root in %lf\n",Rad.Zero);
00076    Radici[rad] = Rad.Zero;
00077    rad++;
00078       }
00079     }
00080     //  }
00081   return rad;
00082 }
00083 double Matematica::Estremo(double a,double b){
00084   double Uno;double Due;double Tre;
00085   //  a<b : Delta=(b-a)/NPassi ? Delta=(a-b)/NPassi;
00086   double Delta=0.;
00087   int NLim = 1000;
00088   Delta=(b-a)/(double)NLim;
00089   Uno=a;Due=b;Tre=0.;
00090   for(int i=0;i<NLim;i++){
00091     if( ASS((Evalx(Tre)-0.)) < PrecMinimo){
00092       break;
00093     }
00094     Tre = Due - (Due - Uno)/(Df(Due,Delta)-Df(Uno,Delta))*Df(Due,Delta);
00095     Due = Tre;
00096     Uno = Due;
00097   }
00098   return Tre;
00099 }
00100 RADICE Matematica::RegulaFalsi(double a,double b){
00101   RADICE Rad; 
00102   double Uno=a;double Due=b;double Tre=0.;
00103   double dIncr = 100.;
00104   double Delta = (b-a)/dIncr;
00105   Uno=a;Due=b;Tre=0.;
00106   FILE *CONTROLLA;
00107   CONTROLLA = fopen("RegulaFalsi.dat","w");
00108   for(int i=0;i<NPassi;i++){
00109     if( ASS(Evalx(Tre)) < PrecMinimo){
00110       break;
00111     }
00112     else if( ASS(Evalx(Due)) < PrecMinimo){
00113       Tre = Due;
00114       break;
00115     }
00116     else if( ASS(Evalx(Uno)) < PrecMinimo){
00117       Tre = Uno;
00118       break;
00119     }
00120     if( Evalx(Due) < 0. && Evalx(Uno) > 0.){
00121       if( Evalx(Due) < 0. &&  Evalx(Uno + Delta) > 0.){
00122    Uno += Delta;
00123       }
00124       if(Evalx(Due - Delta) < 0. && Evalx(Uno) > 0.){
00125    Due -= Delta;
00126       }
00127       else {
00128    dIncr *= 10.;
00129       }
00130     }
00131     else if( Evalx(Due) > 0. && Evalx(Uno) < 0.){
00132       if( Evalx(Due) > 0. && Evalx(Uno+Delta) < 0.){
00133    Uno += Delta;
00134       }
00135       if(Evalx(Due - Delta) > 0.  && Evalx(Uno) < 0.){
00136    Due -= Delta;
00137       }
00138       else {
00139    dIncr *= 10.;
00140       }
00141     }
00142     else{ 
00143       Uno += (b-a) / dIncr;
00144     }
00145     if(Uno > Due){
00146       Uno = a;
00147       dIncr *= 10.;
00148     }
00149     Delta = (b-a)/dIncr;
00150     Tre = (Due + Uno)/2.;
00151     //printf("%d) Uno %g|%g  Due %g|%g  Evalx(Tre) %g|%g Delta %g Incr %g\n",i,Uno,Evalx(Uno),Due,Evalx(Due),Tre,Evalx(Tre),Delta,dIncr);
00152     //fprintf(CONTROLLA,"%lf %lf\n",Tre,Evalx(Tre));
00153   }
00154   if( !(ASS((Evalx(Tre)-0.)) < PrecMinimo)){
00155     Rad.IfRis = 0;
00156   }
00157   else 
00158     Rad.IfRis = 1;
00159   Rad.Zero = Tre;
00160   Rad.iLim = Uno;
00161   Rad.sLim = Due;
00162   fclose(CONTROLLA);
00163   return Rad;
00164 }
00165 RADICE Matematica::Newton(double a){
00166   RADICE Rad;
00167   double Uno=a,Due =0.,Tre=0.;
00168   double m=0.,q=0.;
00169   FILE *CONTROLLA;
00170   CONTROLLA = fopen("Newton.dat","w");
00171   for(int i=0;i<NPassi;i++){
00172     if( ASS((Evalx(Tre))) < PrecMinimo){
00173       break;
00174     }
00175     Due = Uno + 1e-5;
00176     m = (Evalx(Due) - Evalx(Uno))/(Due - Uno);
00177     q = Evalx(Uno) - m*Uno;
00178     Tre = -q/m;
00179     printf("%d) %lf  %lf  %lf\n",i,Uno,Due,Evalx(Tre));
00180     fprintf(CONTROLLA,"%lf %lf\n",Tre,Evalx(Tre));
00181     Uno = Tre;
00182   }
00183   if( !(ASS((Evalx(Tre)-0.)) < PrecMinimo)){
00184     printf("Calculation failed\n");
00185     Rad.IfRis = 0;
00186   }
00187   else 
00188     Rad.IfRis = 1;
00189   Rad.Zero = Tre;
00190   Rad.iLim = Uno;
00191   Rad.sLim = Due;
00192   fclose(CONTROLLA);
00193   return Rad;
00194 }
00195 double Matematica::Gauss(double Media,double Scarto,double x){
00196   return 1./(Scarto*sqrt(DUE_PI))*exp(- .5*SQR((x-Media)/Scarto) );
00197 }
00198 double Matematica::IntegrazioneGauss(double a,double b,double Scarto){
00199   double Delta = (b-a)/(double)NPassi;
00200   double Risp = 0.;
00201   for(int i=0;a+((double)i)*Delta<b;i+=3){//+3 se si sovrappongono
00202     Risp += 3*Delta*( Gauss(0.,Scarto,a+Delta*i)+3*Gauss(0.,Scarto,a+Delta*(i+1))+3*Gauss(0.,Scarto,a+Delta*(i+2))+Gauss(0.,Scarto,a+Delta*(i+3)) )/8;
00203   }
00204   return Risp;
00205 }
00206 void Matematica::SquareGradient(double *st,double *sw,int NMass){
00207   DerO4(st,sw,NMass);
00208 }
00209 void Matematica::NormalizeVect(double *st,int NMass){
00210   double Norm = 0.;
00211   for(int n=0;n<NMass;n++){
00212     Norm += SQR(st[n]);
00213   }
00214   Norm = 1./sqrt(Norm);
00215   for(int n=0;n<NMass;n++){
00216     st[n] = st[n]*Norm;
00217   }  
00218 }
00219 double Matematica::Norm(double *st,int NMass){
00220   double Norm = 0.;
00221   for(int n=0;n<NMass;n++){
00222     Norm += SQR(st[n]);
00223   }
00224   return sqrt(Norm);
00225 }
00226 void Matematica::Modulo(double *st,double *sw,int NMass){
00227   for(int n=0;n<NMass;n++){
00228     sw[n] = POS(st[n]);
00229   }
00230 }
00231 double Matematica::Fattoriale(int n){
00232   if(n < 0) { printf("Il fattoriale di numeri negativi non ha senso\n"); return 0;}
00233   if( n == 0 ) return 1;
00234   int Ris=1;
00235   for(int i=n;i>0;i--){
00236     //    printf("Fatt %d\n",i);
00237     Ris *= i;
00238   }
00239   return (double) Ris;
00240 }
00241 double Matematica::Gamma(int n){
00242   return Fattoriale(n-1);
00243 }
00244 double Matematica::Elevato(double x,int Volte){
00245   double Risp=1.;  
00246   double Moltx = Volte >= 0 ? x : 1./x;
00247   for(int v=0;v<Volte;v++)
00248     Risp *= Moltx;
00249   return Risp;
00250 }
00251 double Matematica::Bessel(double Val,int Ord){
00252   int NOrd = POS(Ord);
00253   int NMax = 10;
00254   double Risp=0.;
00255   for(int n=0;n<NMax;n++){
00256     //    printf("  %d\n",n);
00257     Risp += Elevato(-1.,n)/(Fattoriale(n)*Gamma(NOrd+n+1))*Elevato(.5*Val,2*n+NOrd);
00258   }
00259   return Ord>=0 ? Risp : Elevato(-1.,Ord)*Risp;
00260 }
00261 double Matematica::Neumann(double Val,int Ord){
00262   double Angolo = Ord*.5*DUE_PI+.0001;
00263   return Bessel(Val,Ord)*(cos(Angolo) - Elevato(-1.,Ord))/sin(Angolo);
00264 }
00265 double Matematica::QuasiBessel(double Val,int Ord){
00266   if( Val > (double)Ord +1.)
00267     return sqrt((2./(PI*Val))) *cos( Val- (2.*Ord+1.)*PI*.25);
00268   else
00269     return Elevato(.5*Val,Ord)/(Fattoriale(Ord));
00270 }
00271 double Matematica::QuasiNeumann(double Val,int Ord){
00272   if(Ord == 0)
00273     return 2./PI*log(Val);
00274   if(Val > (double) Ord+1.)
00275     return sqrt(2./(PI*Val))*sin(Val-(2.*Ord+1.)*PI*.25);
00276   else
00277     return Elevato(2,Ord)*Fattoriale(Ord-1)/PI*Elevato(Val,-Ord);
00278 }
00279 double Matematica::Segno(int n){
00280   return (n%2)==1 ? -1. : 1.;
00281 }
00282 double Matematica::WeightFunction(double x,double a){
00283   double Risp=2.*x*x*x - 3.*(a+1.)*x*x-3.*a*a+1.;
00284   Risp /= CUB((1.-a));
00285   return Risp;
00286 }
00287 double Matematica::WeightFunction2(double x,double a){
00288   double Risp= -a*x*x + 1. + 4./3.*PI*a;
00289   Risp /= CUB((1.-a));
00290   return Risp;
00291 }
00292 double Matematica::LJHamaker(double Rad,double RadNp,double Theta){
00293   //return 1./CUB(Rad-RadNp);
00294   double Num = 2.*DUE_PI*QUAD(RadNp)*sin(Theta);
00295   double Den = QUAD(RadNp) + QUAD(Rad+RadNp) - 2.*(Rad+RadNp)*RadNp*cos(Theta);
00296   return Num/CUB(Den);// - Num/CUB(CUB(Den));
00297 }
00298 double Matematica::LJ39(double r,double r_np){
00299   return pow(1./(r-r_np),9.) - pow(1./(r-r_np),3.);
00300 }
00301 double Matematica::LJHamaker(double Rad,double RadNp){
00302   double ThetaMax = PI;
00303   double ThetaMin = 0.;
00304   double ThetaDelta = (ThetaMax-ThetaMin)/100.;
00305   double Risp = 0.;
00306   for(double Theta = ThetaMin;Theta<ThetaMax;Theta+=ThetaDelta){
00307     Risp += ThetaDelta*.5*(LJHamaker(Rad,RadNp,Theta) + LJHamaker(Rad,RadNp,Theta+ThetaDelta) );
00308   }
00309   return Risp;
00310 }
00311 double Matematica::LJHamakerCum(double Rad,double RadNpMin,double RadNpMax){
00312   double RadNpDelta = (RadNpMax-RadNpMin)/100.;
00313   double Risp = 0.;
00314   for(double RadNp = RadNpMin;RadNp<RadNpMax;RadNp+=RadNpDelta){
00315     //Risp += RadNpDelta*LJHamaker(Rad,RadNp+RadNpDelta*0.5);
00316     Risp += RadNpDelta*.5*(LJHamaker(Rad,RadNp) + LJHamaker(Rad,RadNp+RadNpDelta) );
00317     //Risp += 3.*RadNpDelta*(LJHamaker(Rad,RadNp) + 3.*LJHamaker(Rad,RadNp+RadNpDelta) + 3.*LJHamaker(Rad,RadNp+2.*RadNpDelta) + 3.*LJHamaker(Rad,RadNp+3.*RadNpDelta) )/8.;
00318   }
00319   return Risp;
00320 }
00321 double Potenziale(double Dist, double RadNp){
00322   double Rad  = Dist + RadNp;
00323   double Pre1 = 2./(12.*Rad);
00324   double Post1 = RadNp / CUB(Rad+RadNp) + RadNp/CUB(Rad-RadNp);
00325   double Pre2 = 1./(12.*Rad);
00326   double Post2 = 1./QUAD(Rad+RadNp) - 1./QUAD(Rad+RadNp);
00327   return PI*(Pre1*Post1 + Pre2*Post2);
00328 }
00329 double Potenziale2(double Dist, double RadNp){
00330   double Rad = Dist + RadNp;
00331   double Pre1 = -( QUAD(RadNp) - QUAD(Rad) )/(4.*Rad);
00332   double Post1 = 1./QUAD(QUAD(Rad+RadNp)) - 1./QUAD(QUAD(Rad-RadNp));
00333   double Pre2  = -2./(3.);
00334   double Post2 = 1./CUB(Rad+RadNp) - 1./CUB(Rad+RadNp);
00335   double Pre3  = 1./(2.*Rad);
00336   double Post3 = 1./QUAD(Rad+RadNp) - 1./QUAD(Rad+RadNp);
00337   return PI*(Pre1*Post1 + Pre2*Post2 + Pre3*Post3);
00338 }
00339 void Matematica::IntegraA3(){
00340   double RadNpMin = 0.;
00341   double RadNpMax = .001;
00342   double RadNpDelta = (RadNpMax - RadNpMin)/20.;
00343   double RadMin = 0.1;
00344   double RadMax = 3.;
00345   double RadDelta = (RadMax-RadMin)/100.;
00346   char *FileName = (char *)calloc(60,sizeof(char));
00347   double RadNp = RadNpMax;
00348   //for(double RadNp = RadNpMin;RadNp <= RadNpMax;RadNp += RadNpDelta)
00349     {
00350     sprintf(FileName,"Potential%0.2f.dat",RadNp);
00351     FILE *POT = fopen(FileName,"w");
00352     for(double Rad = RadMin;Rad < RadMax;Rad += RadDelta){
00353       fprintf(POT,"%lf %g %g %g \n",Rad,LJHamaker(Rad,RadNp),LJHamakerCum(Rad,RadNpMin,RadNp),Potenziale(Rad,RadNp) );
00354     }
00355     fclose(POT);
00356   }
00357 }