Allink  v0.1
MatematicaFunc.cpp
1 #include "../include/Matematica.h"
2 
3 double Matematica::ContactAngle(double x){
4  double Ris;
5  //Ris = PreFact*pow(3./(2.+3.*sin(x)-CUB(sin(x)) ), 4./3. )*QUAD(QUAD(cos(x))) - Ypsilon;
6  Ris = PreFact*pow(1-cos(x), -4./3. )*pow(sin(x),2.) - Ypsilon;
7  return Ris;
8 }
9 double Matematica::fProva(double x){
10  double Ris;
11  //Ris = pow(x,5.)/((exp(x)-1)*(1-exp(-x)));
12  // Ris = pow(x,5.)*exp(x)/QUAD(exp(x)-1);
13  // Ris = 1.*pow(x-150,4.)-1.*pow(x-150,5.)+x;
14  // Ris = pow(x-150,2.)+x-150;
15  Ris = 1./(13.9332529115775*x*x - 4.*13.9332529115775*0.13710248135043*x + 7.86306085419497) + .5/(-0.0805855688278211*x + 10.6925178203198);
16  //Ris = 1./(12.5342*x*x) + .5/(-0.207807999999993*x + 17.4283);
17  return Ris;
18 }
19 double Matematica::F(double T,double TD){
20  return pow((T/TD),5.)*Integrazione(.00001,TD/T);
21 }
22 double Matematica::Integrazione(double a,double b){
23  double Delta = (b-a)/(double)NPassi;
24  double Risp = 0.;
25  for(int i=0;a+((double)i)*Delta<b;i+=3){//+3 se si sovrappongono
26  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.;
27  }
28  return Risp;
29 }
30 double Matematica::Integrazione(double *st,double *sw,int NMass){
31  double NMassInv = 1/((double)NMass);
32  sw[0] = st[0];
33  for(int i=1;i<NMass;i++){
34  sw[i] = st[i] + sw[i-1];
35  }
36  return sw[NMass-1];
37 }
38 double Matematica::Df(double x,double Delta){
39  return ( Evalx(x+Delta)-Evalx(x) )/Delta;
40 }
41 void Matematica::Derivata(double *st,double *sw,int NMass){
42  for(int i=0;i<NMass-1;i++){
43  sw[i] = (st[i+1]-st[i])/2.;
44  }
45 }
46 void Matematica::DerO4(double *st,double *sw,int NMass){
47  sw[0] = 0.;
48  sw[NMass-1] = 0.;
49  for(int i=1;i<NMass-1;i++){
50  if(i<2)
51  sw[i] = (st[i+1] - st[i-1])/2.;
52  else if (i < NMass - 2)
53  sw[i] = st[i-2] - 8.*st[i-1] + 8.*st[i+1] - st[i+2];
54  else if (i < NMass - 1)
55  sw[i] = (st[i+1] - st[i-1])/2.;
56  }
57 }
58 int Matematica::Zeri(double a,double b,double *Radici,int NRadici){
59  double Uno;double Due;double Delta;
60  int rad=0;
61  RADICE Rad;
62  if( a > b ){
63  Uno = b ; Due = a;
64  }
65  else{
66  Uno = a ; Due = b;
67  }
68  Delta = (Due - Uno)/2.;
69  // for(int p=0;p<4;p++){
70  for(int i=0;i<NRadici*2;i++){
71  Delta = (Due - Uno)/(2*NRadici);
72  Rad = RegulaFalsi(Uno+i*Delta,Due-(2*NRadici-i-1)*Delta);
73  //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);
74  if(Rad.IfRis == 1){
75  printf("Found a root in %lf\n",Rad.Zero);
76  Radici[rad] = Rad.Zero;
77  rad++;
78  }
79  }
80  // }
81  return rad;
82 }
83 double Matematica::Estremo(double a,double b){
84  double Uno;double Due;double Tre;
85  // a<b : Delta=(b-a)/NPassi ? Delta=(a-b)/NPassi;
86  double Delta=0.;
87  int NLim = 1000;
88  Delta=(b-a)/(double)NLim;
89  Uno=a;Due=b;Tre=0.;
90  for(int i=0;i<NLim;i++){
91  if( ASS((Evalx(Tre)-0.)) < PrecMinimo){
92  break;
93  }
94  Tre = Due - (Due - Uno)/(Df(Due,Delta)-Df(Uno,Delta))*Df(Due,Delta);
95  Due = Tre;
96  Uno = Due;
97  }
98  return Tre;
99 }
100 RADICE Matematica::RegulaFalsi(double a,double b){
101  RADICE Rad;
102  double Uno=a;double Due=b;double Tre=0.;
103  double dIncr = 100.;
104  double Delta = (b-a)/dIncr;
105  Uno=a;Due=b;Tre=0.;
106  FILE *CONTROLLA;
107  CONTROLLA = fopen("RegulaFalsi.dat","w");
108  for(int i=0;i<NPassi;i++){
109  if( ASS(Evalx(Tre)) < PrecMinimo){
110  break;
111  }
112  else if( ASS(Evalx(Due)) < PrecMinimo){
113  Tre = Due;
114  break;
115  }
116  else if( ASS(Evalx(Uno)) < PrecMinimo){
117  Tre = Uno;
118  break;
119  }
120  if( Evalx(Due) < 0. && Evalx(Uno) > 0.){
121  if( Evalx(Due) < 0. && Evalx(Uno + Delta) > 0.){
122  Uno += Delta;
123  }
124  if(Evalx(Due - Delta) < 0. && Evalx(Uno) > 0.){
125  Due -= Delta;
126  }
127  else {
128  dIncr *= 10.;
129  }
130  }
131  else if( Evalx(Due) > 0. && Evalx(Uno) < 0.){
132  if( Evalx(Due) > 0. && Evalx(Uno+Delta) < 0.){
133  Uno += Delta;
134  }
135  if(Evalx(Due - Delta) > 0. && Evalx(Uno) < 0.){
136  Due -= Delta;
137  }
138  else {
139  dIncr *= 10.;
140  }
141  }
142  else{
143  Uno += (b-a) / dIncr;
144  }
145  if(Uno > Due){
146  Uno = a;
147  dIncr *= 10.;
148  }
149  Delta = (b-a)/dIncr;
150  Tre = (Due + Uno)/2.;
151  //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);
152  //fprintf(CONTROLLA,"%lf %lf\n",Tre,Evalx(Tre));
153  }
154  if( !(ASS((Evalx(Tre)-0.)) < PrecMinimo)){
155  Rad.IfRis = 0;
156  }
157  else
158  Rad.IfRis = 1;
159  Rad.Zero = Tre;
160  Rad.iLim = Uno;
161  Rad.sLim = Due;
162  fclose(CONTROLLA);
163  return Rad;
164 }
166  RADICE Rad;
167  double Uno=a,Due =0.,Tre=0.;
168  double m=0.,q=0.;
169  FILE *CONTROLLA;
170  CONTROLLA = fopen("Newton.dat","w");
171  for(int i=0;i<NPassi;i++){
172  if( ASS((Evalx(Tre))) < PrecMinimo){
173  break;
174  }
175  Due = Uno + 1e-5;
176  m = (Evalx(Due) - Evalx(Uno))/(Due - Uno);
177  q = Evalx(Uno) - m*Uno;
178  Tre = -q/m;
179  printf("%d) %lf %lf %lf\n",i,Uno,Due,Evalx(Tre));
180  fprintf(CONTROLLA,"%lf %lf\n",Tre,Evalx(Tre));
181  Uno = Tre;
182  }
183  if( !(ASS((Evalx(Tre)-0.)) < PrecMinimo)){
184  printf("Calculation failed\n");
185  Rad.IfRis = 0;
186  }
187  else
188  Rad.IfRis = 1;
189  Rad.Zero = Tre;
190  Rad.iLim = Uno;
191  Rad.sLim = Due;
192  fclose(CONTROLLA);
193  return Rad;
194 }
195 double Matematica::Gauss(double Media,double Scarto,double x){
196  return 1./(Scarto*sqrt(DUE_PI))*exp(- .5*SQR((x-Media)/Scarto) );
197 }
198 double Matematica::IntegrazioneGauss(double a,double b,double Scarto){
199  double Delta = (b-a)/(double)NPassi;
200  double Risp = 0.;
201  for(int i=0;a+((double)i)*Delta<b;i+=3){//+3 se si sovrappongono
202  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;
203  }
204  return Risp;
205 }
206 void Matematica::SquareGradient(double *st,double *sw,int NMass){
207  DerO4(st,sw,NMass);
208 }
209 void Matematica::NormalizeVect(double *st,int NMass){
210  double Norm = 0.;
211  for(int n=0;n<NMass;n++){
212  Norm += SQR(st[n]);
213  }
214  Norm = 1./sqrt(Norm);
215  for(int n=0;n<NMass;n++){
216  st[n] = st[n]*Norm;
217  }
218 }
219 double Matematica::Norm(double *st,int NMass){
220  double Norm = 0.;
221  for(int n=0;n<NMass;n++){
222  Norm += SQR(st[n]);
223  }
224  return sqrt(Norm);
225 }
226 void Matematica::Modulo(double *st,double *sw,int NMass){
227  for(int n=0;n<NMass;n++){
228  sw[n] = POS(st[n]);
229  }
230 }
232  if(n < 0) { printf("Il fattoriale di numeri negativi non ha senso\n"); return 0;}
233  if( n == 0 ) return 1;
234  int Ris=1;
235  for(int i=n;i>0;i--){
236  // printf("Fatt %d\n",i);
237  Ris *= i;
238  }
239  return (double) Ris;
240 }
241 double Matematica::Gamma(int n){
242  return Fattoriale(n-1);
243 }
244 double Matematica::Elevato(double x,int Volte){
245  double Risp=1.;
246  double Moltx = Volte >= 0 ? x : 1./x;
247  for(int v=0;v<Volte;v++)
248  Risp *= Moltx;
249  return Risp;
250 }
251 double Matematica::Bessel(double Val,int Ord){
252  int NOrd = POS(Ord);
253  int NMax = 10;
254  double Risp=0.;
255  for(int n=0;n<NMax;n++){
256  // printf(" %d\n",n);
257  Risp += Elevato(-1.,n)/(Fattoriale(n)*Gamma(NOrd+n+1))*Elevato(.5*Val,2*n+NOrd);
258  }
259  return Ord>=0 ? Risp : Elevato(-1.,Ord)*Risp;
260 }
261 double Matematica::Neumann(double Val,int Ord){
262  double Angolo = Ord*.5*DUE_PI+.0001;
263  return Bessel(Val,Ord)*(cos(Angolo) - Elevato(-1.,Ord))/sin(Angolo);
264 }
265 double Matematica::QuasiBessel(double Val,int Ord){
266  if( Val > (double)Ord +1.)
267  return sqrt((2./(PI*Val))) *cos( Val- (2.*Ord+1.)*PI*.25);
268  else
269  return Elevato(.5*Val,Ord)/(Fattoriale(Ord));
270 }
271 double Matematica::QuasiNeumann(double Val,int Ord){
272  if(Ord == 0)
273  return 2./PI*log(Val);
274  if(Val > (double) Ord+1.)
275  return sqrt(2./(PI*Val))*sin(Val-(2.*Ord+1.)*PI*.25);
276  else
277  return Elevato(2,Ord)*Fattoriale(Ord-1)/PI*Elevato(Val,-Ord);
278 }
279 double Matematica::Segno(int n){
280  return (n%2)==1 ? -1. : 1.;
281 }
282 double Matematica::WeightFunction(double x,double a){
283  double Risp=2.*x*x*x - 3.*(a+1.)*x*x-3.*a*a+1.;
284  Risp /= CUB((1.-a));
285  return Risp;
286 }
287 double Matematica::WeightFunction2(double x,double a){
288  double Risp= -a*x*x + 1. + 4./3.*PI*a;
289  Risp /= CUB((1.-a));
290  return Risp;
291 }
292 double Matematica::LJHamaker(double Rad,double RadNp,double Theta){
293  //return 1./CUB(Rad-RadNp);
294  double Num = 2.*DUE_PI*QUAD(RadNp)*sin(Theta);
295  double Den = QUAD(RadNp) + QUAD(Rad+RadNp) - 2.*(Rad+RadNp)*RadNp*cos(Theta);
296  return Num/CUB(Den);// - Num/CUB(CUB(Den));
297 }
298 double Matematica::LJ39(double r,double r_np){
299  return pow(1./(r-r_np),9.) - pow(1./(r-r_np),3.);
300 }
301 double Matematica::LJHamaker(double Rad,double RadNp){
302  double ThetaMax = PI;
303  double ThetaMin = 0.;
304  double ThetaDelta = (ThetaMax-ThetaMin)/100.;
305  double Risp = 0.;
306  for(double Theta = ThetaMin;Theta<ThetaMax;Theta+=ThetaDelta){
307  Risp += ThetaDelta*.5*(LJHamaker(Rad,RadNp,Theta) + LJHamaker(Rad,RadNp,Theta+ThetaDelta) );
308  }
309  return Risp;
310 }
311 double Matematica::LJHamakerCum(double Rad,double RadNpMin,double RadNpMax){
312  double RadNpDelta = (RadNpMax-RadNpMin)/100.;
313  double Risp = 0.;
314  for(double RadNp = RadNpMin;RadNp<RadNpMax;RadNp+=RadNpDelta){
315  //Risp += RadNpDelta*LJHamaker(Rad,RadNp+RadNpDelta*0.5);
316  Risp += RadNpDelta*.5*(LJHamaker(Rad,RadNp) + LJHamaker(Rad,RadNp+RadNpDelta) );
317  //Risp += 3.*RadNpDelta*(LJHamaker(Rad,RadNp) + 3.*LJHamaker(Rad,RadNp+RadNpDelta) + 3.*LJHamaker(Rad,RadNp+2.*RadNpDelta) + 3.*LJHamaker(Rad,RadNp+3.*RadNpDelta) )/8.;
318  }
319  return Risp;
320 }
321 double Potenziale(double Dist, double RadNp){
322  double Rad = Dist + RadNp;
323  double Pre1 = 2./(12.*Rad);
324  double Post1 = RadNp / CUB(Rad+RadNp) + RadNp/CUB(Rad-RadNp);
325  double Pre2 = 1./(12.*Rad);
326  double Post2 = 1./QUAD(Rad+RadNp) - 1./QUAD(Rad+RadNp);
327  return PI*(Pre1*Post1 + Pre2*Post2);
328 }
329 double Potenziale2(double Dist, double RadNp){
330  double Rad = Dist + RadNp;
331  double Pre1 = -( QUAD(RadNp) - QUAD(Rad) )/(4.*Rad);
332  double Post1 = 1./QUAD(QUAD(Rad+RadNp)) - 1./QUAD(QUAD(Rad-RadNp));
333  double Pre2 = -2./(3.);
334  double Post2 = 1./CUB(Rad+RadNp) - 1./CUB(Rad+RadNp);
335  double Pre3 = 1./(2.*Rad);
336  double Post3 = 1./QUAD(Rad+RadNp) - 1./QUAD(Rad+RadNp);
337  return PI*(Pre1*Post1 + Pre2*Post2 + Pre3*Post3);
338 }
340  double RadNpMin = 0.;
341  double RadNpMax = .001;
342  double RadNpDelta = (RadNpMax - RadNpMin)/20.;
343  double RadMin = 0.1;
344  double RadMax = 3.;
345  double RadDelta = (RadMax-RadMin)/100.;
346  char *FileName = (char *)calloc(60,sizeof(char));
347  double RadNp = RadNpMax;
348  //for(double RadNp = RadNpMin;RadNp <= RadNpMax;RadNp += RadNpDelta)
349  {
350  sprintf(FileName,"Potential%0.2f.dat",RadNp);
351  FILE *POT = fopen(FileName,"w");
352  for(double Rad = RadMin;Rad < RadMax;Rad += RadDelta){
353  fprintf(POT,"%lf %g %g %g \n",Rad,LJHamaker(Rad,RadNp),LJHamakerCum(Rad,RadNpMin,RadNp),Potenziale(Rad,RadNp) );
354  }
355  fclose(POT);
356  }
357 }
double Ypsilon
External parameter to calculate the contact angle.
Definition: Matematica.h:330
void SquareGradient(double *st, double *sw, int NMass)
Square of the gradient.
void DerO4(double *st, double *sw, int NMass)
Derivate O(4) of.
void NormalizeVect(double *st, int NMass)
Normalize.
double Integrazione(double *Punti, double *sw, int NMass)
Integral of.
int Zeri(double a, double b, double *Radici, int NRadici)
Find the.
double sLim
Superior limit.
double fProva(double x)
Trial function.
double QuasiBessel(double Val, int Ord)
A faster Bessel.
void Derivata(double *st, double *sw, int NMass)
Derivate of.
double Bessel(double Val, int Ord)
Bessel function.
int IfRis
If the zero was founds.
void Modulo(double *st, double *sw, int NMass)
Compute the modulus.
double QuasiNeumann(double Val, int Ord)
A faster Neumann.
double Df(double x, double Delta)
Boh.
double iLim
Iferior limit.
double Zero
Point of the zero.
double LJ39(double r, double r_np)
Integration of the LJ 6 term.
double ContactAngle(double x)
Definition of the contact angle.
double WeightFunction2(double x, double a)
Definition of a weighting function.
double PreFact
External parameter in the definition of the contact angle.
Definition: Matematica.h:332
RADICE Newton(double a)
Use Newton to find the roots.
void IntegraA3()
Perform a integration of a LJ6 Potential.
double Gauss(double Media, double Scarto, double x)
Gaussian.
double Elevato(double x, int Volte)
Integer power.
double Estremo(double a, double b)
Other algorithm to find the roots.
Where a root was searched.
double F(double TD, double T)
Boh.
double Gamma(int n)
Euler&#39;s gamma.
double Evalx(double x)
Pointer to a generic function.
Definition: Matematica.h:113
double Norm(double *st, int NMass)
Norm of an array.
double LJHamaker(double r, double r_np, double theta)
Integration of the LJ 6 term.
double Fattoriale(int n)
Compute the factorial.
double IntegrazioneGauss(double a, double b, double Scarto)
Itegrate a Gaussian.
RADICE RegulaFalsi(double a, double b)
Use regula falsi algorithm to find the roots.
double Neumann(double Val, int Ord)
Neumann function.
double Segno(int n)
Sign of -^n.
double WeightFunction(double x, double a)
Definition of a weighting function.
double LJHamakerCum(double Rad, double RadNpMin, double RadNpMax)
Integrate over r_np up to RadNp.