source: sasview/sansmodels/src/libigor/libStructureFactor.c @ 196608d

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload
Last change on this file since 196608d was 572beba, checked in by Jae Cho <jhjcho@…>, 15 years ago

fixed mistake
  • Property mode set to 100644
File size: 18.6 KB
Line 
1/*      SimpleFit.c
2
3A simplified project designed to act as a template for your curve fitting function.
4The fitting function is a simple polynomial. It works but is of no practical use.
5*/
6
7#include "StandardHeaders.h"                    // Include ANSI headers, Mac headers, IgorXOP.h, XOP.h and XOPSupport.h
8#include "libStructureFactor.h"
9//Hard Sphere Structure Factor
10//
11double
12HardSphereStruct(double dp[], double q)
13{
14        double denom,dnum,alpha,beta,gamm,a,asq,ath,afor,rca,rsa;
15        double calp,cbeta,cgam,prefac,c,vstruc;
16        double r,phi,struc;
17
18        r = dp[0];
19        phi = dp[1];
20        //  compute constants
21        denom = pow((1.0-phi),4);
22        dnum = pow((1.0 + 2.0*phi),2);
23        alpha = dnum/denom;
24        beta = -6.0*phi*pow((1.0 + phi/2.0),2)/denom;
25        gamm = 0.50*phi*dnum/denom;
26        //
27        //  calculate the structure factor
28        //
29        a = 2.0*q*r;
30        asq = a*a;
31        ath = asq*a;
32        afor = ath*a;
33        rca = cos(a);
34        rsa = sin(a);
35        calp = alpha*(rsa/asq - rca/a);
36        cbeta = beta*(2.0*rsa/asq - (asq - 2.0)*rca/ath - 2.0/ath);
37        cgam = gamm*(-rca/a + (4.0/a)*((3.0*asq - 6.0)*rca/afor + (asq - 6.0)*rsa/ath + 6.0/afor));
38        prefac = -24.0*phi/a;
39        c = prefac*(calp + cbeta + cgam);
40        vstruc = 1.0/(1.0-c);
41        struc = vstruc;
42
43        return(struc);
44}
45
46//Sticky Hard Sphere Structure Factor
47//
48double
49StickyHS_Struct(double dp[], double q)
50{
51        double qv;
52        double rad,phi,eps,tau,eta;
53        double sig,aa,etam1,qa,qb,qc,radic;
54        double lam,lam2,test,mu,alpha,beta;
55        double kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,sq;
56
57        qv= q;
58        rad = dp[0];
59        phi = dp[1];
60        eps = dp[2];
61        tau = dp[3];
62
63        eta = phi/(1.0-eps)/(1.0-eps)/(1.0-eps);
64
65        sig = 2.0 * rad;
66        aa = sig/(1.0 - eps);
67        etam1 = 1.0 - eta;
68        //C
69        //C  SOLVE QUADRATIC FOR LAMBDA
70        //C
71        qa = eta/12.0;
72        qb = -1.0*(tau + eta/(etam1));
73        qc = (1.0 + eta/2.0)/(etam1*etam1);
74        radic = qb*qb - 4.0*qa*qc;
75        if(radic<0) {
76                //if(x>0.01 && x<0.015)
77                //      Print "Lambda unphysical - both roots imaginary"
78                //endif
79                return(-1.0);
80        }
81        //C   KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL
82        lam = (-1.0*qb-sqrt(radic))/(2.0*qa);
83        lam2 = (-1.0*qb+sqrt(radic))/(2.0*qa);
84        if(lam2<lam) {
85                lam = lam2;
86        }
87        test = 1.0 + 2.0*eta;
88        mu = lam*eta*etam1;
89        if(mu>test) {
90                //if(x>0.01 && x<0.015)
91                // Print "Lambda unphysical mu>test"
92                //endif
93                return(-1.0);
94        }
95        alpha = (1.0 + 2.0*eta - mu)/(etam1*etam1);
96        beta = (mu - 3.0*eta)/(2.0*etam1*etam1);
97        //C
98        //C   CALCULATE THE STRUCTURE FACTOR
99        //C
100        kk = qv*aa;
101        k2 = kk*kk;
102        k3 = kk*k2;
103        ds = sin(kk);
104        dc = cos(kk);
105        aq1 = ((ds - kk*dc)*alpha)/k3;
106        aq2 = (beta*(1.0-dc))/k2;
107        aq3 = (lam*ds)/(12.0*kk);
108        aq = 1.0 + 12.0*eta*(aq1+aq2-aq3);
109        //
110        bq1 = alpha*(0.5/kk - ds/k2 + (1.0 - dc)/k3);
111        bq2 = beta*(1.0/kk - ds/k2);
112        bq3 = (lam/12.0)*((1.0 - dc)/kk);
113        bq = 12.0*eta*(bq1+bq2-bq3);
114        //
115        sq = 1.0/(aq*aq +bq*bq);
116
117        return(sq);
118}
119
120
121
122//     SUBROUTINE SQWELL: CALCULATES THE STRUCTURE FACTOR FOR A
123//                        DISPERSION OF MONODISPERSE HARD SPHERES
124//     IN THE Mean Spherical APPROXIMATION ASSUMING THE SPHERES
125//     INTERACT THROUGH A SQUARE WELL POTENTIAL.
126//** not the best choice of closure ** see note below
127//     REFS:  SHARMA,SHARMA, PHYSICA 89A,(1977),212
128double
129SquareWellStruct(double dp[], double q)
130{
131        double req,phis,edibkb,lambda,struc;
132        double sigma,eta,eta2,eta3,eta4,etam1,etam14,alpha,beta,gamm;
133        double x,sk,sk2,sk3,sk4,t1,t2,t3,t4,ck;
134
135        x= q;
136
137        req = dp[0];
138        phis = dp[1];
139        edibkb = dp[2];
140        lambda = dp[3];
141
142        sigma = req*2.;
143        eta = phis;
144        eta2 = eta*eta;
145        eta3 = eta*eta2;
146        eta4 = eta*eta3;
147        etam1 = 1. - eta;
148        etam14 = etam1*etam1*etam1*etam1;
149        alpha = (  pow((1. + 2.*eta),2) + eta3*( eta-4.0 )  )/etam14;
150        beta = -(eta/3.0) * ( 18. + 20.*eta - 12.*eta2 + eta4 )/etam14;
151        gamm = 0.5*eta*( pow((1. + 2.*eta),2) + eta3*(eta-4.) )/etam14;
152        //
153        //  calculate the structure factor
154
155        sk = x*sigma;
156        sk2 = sk*sk;
157        sk3 = sk*sk2;
158        sk4 = sk3*sk;
159        t1 = alpha * sk3 * ( sin(sk) - sk * cos(sk) );
160        t2 = beta * sk2 * ( 2.*sk*sin(sk) - (sk2-2.)*cos(sk) - 2.0 );
161        t3 =   ( 4.0*sk3 - 24.*sk ) * sin(sk);
162        t3 = t3 - ( sk4 - 12.0*sk2 + 24.0 )*cos(sk) + 24.0;
163        t3 = gamm*t3;
164        t4 = -edibkb*sk3*(sin(lambda*sk) - lambda*sk*cos(lambda*sk)+ sk*cos(sk) - sin(sk) );
165        ck =  -24.0*eta*( t1 + t2 + t3 + t4 )/sk3/sk3;
166        struc  = 1./(1.-ck);
167
168        return(struc);
169}
170
171// Hayter-Penfold (rescaled) MSA structure factor for screened Coulomb interactions
172//
173double
174HayterPenfoldMSA(double dp[], double q)
175{
176        double Elcharge=1.602189e-19;           // electron charge in Coulombs (C)
177        double kB=1.380662e-23;                         // Boltzman constant in J/K
178        double FrSpPerm=8.85418782E-12; //Permittivity of free space in C^2/(N m^2)
179        double SofQ, QQ, Qdiam, Vp, csalt, ss;
180        double VolFrac, SIdiam, diam, Kappa, cs, IonSt;
181        double dialec, Perm, Beta, Temp, zz, charge;
182        double pi;
183        int ierr;
184
185        pi = 4.0*atan(1.);
186        QQ= q;
187
188        diam=dp[0];             //in dp[0] coming from python is radius : cahnged on Mar. 30, 2009
189        zz = dp[1];             //# of charges
190        VolFrac=dp[2];
191        Temp=dp[3];             //in degrees Kelvin
192        csalt=dp[4];            //in molarity
193        dialec=dp[5];           // unitless
194                                                ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
195                                                //////////////////////////// convert to USEFUL inputs in SI units                                                //
196                                                ////////////////////////////    NOTE: easiest to do EVERYTHING in SI units                               //
197                                                ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
198        Beta=1.0/(kB*Temp);             // in Joules^-1
199        Perm=dialec*FrSpPerm;   //in C^2/(N  m^2)
200        charge=zz*Elcharge;             //in Coulomb (C)
201        SIdiam = diam*1E-10;            //in m
202        Vp=4.0*pi/3.0*(SIdiam/2.0)*(SIdiam/2.0)*(SIdiam/2.0);   //in m^3
203        cs=csalt*6.022E23*1E3;  //# salt molecules/m^3
204
205        //         Compute the derived values of :
206        //                       Ionic strength IonSt (in C^2/m^3)
207        //                      Kappa (Debye-Huckel screening length in m)
208        //      and             gamma Exp(-k)
209        IonSt=0.5 * Elcharge*Elcharge*(zz*VolFrac/Vp+2.0*cs);
210        Kappa=sqrt(2*Beta*IonSt/Perm);     //Kappa calc from Ionic strength
211                                                                           //   Kappa=2/SIdiam                                  // Use to compare with HP paper
212        gMSAWave[5]=Beta*charge*charge/(pi*Perm*SIdiam*pow((2.0+Kappa*SIdiam),2));
213
214        //         Finally set up dimensionless parameters
215        Qdiam=QQ*diam;
216        gMSAWave[6] = Kappa*SIdiam;
217        gMSAWave[4] = VolFrac;
218
219        //Function sqhpa(qq)  {this is where Hayter-Penfold program began}
220
221        //       FIRST CALCULATE COUPLING
222
223        ss=pow(gMSAWave[4],(1.0/3.0));
224        gMSAWave[9] = 2.0*ss*gMSAWave[5]*exp(gMSAWave[6]-gMSAWave[6]/ss);
225
226        //        CALCULATE COEFFICIENTS, CHECK ALL IS WELL
227        //        AND IF SO CALCULATE S(Q*SIG)
228
229        ierr=0;
230        ierr=sqcoef(ierr);
231        if (ierr>=0) {
232                SofQ=sqhcal(Qdiam);
233        }else{
234        //SofQ=NaN;
235                SofQ=-1.0;
236                //      print "Error Level = ",ierr
237                //      print "Please report HPMSA problem with above error code"
238        }
239
240        return(SofQ);
241}
242
243
244
245/////////////////////////////////////////////////////////////
246/////////////////////////////////////////////////////////////
247//
248//
249//      CALCULATES RESCALED VOLUME FRACTION AND CORRESPONDING
250//      COEFFICIENTS FOR "SQHPA"
251//
252//      JOHN B. HAYTER   (I.L.L.)    14-SEP-81
253//
254//      ON EXIT:
255//
256//      SETA IS THE RESCALED VOLUME FRACTION
257//      SGEK IS THE RESCALED CONTACT POTENTIAL
258//      SAK IS THE RESCALED SCREENING CONSTANT
259//      A,B,C,F,U,V ARE THE MSA COEFFICIENTS
260//      G1= G(1+) IS THE CONTACT VALUE OF G(R/SIG):
261//      FOR THE GILLAN CONDITION, THE DIFFERENCE FROM
262//      ZERO INDICATES THE COMPUTATIONAL ACCURACY.
263//
264//      IR > 0:    NORMAL EXIT,  IR IS THE NUMBER OF ITERATIONS.
265//         < 0:    FAILED TO CONVERGE
266//
267int
268sqcoef(int ir)
269{
270        int itm=40,ix,ig,ii;
271        double acc=5.0E-6,del,e1,e2,f1,f2;
272
273        //      WAVE gMSAWave = $"root:HayPenMSA:gMSAWave"
274        f1=0;           //these were never properly initialized...
275        f2=0;
276
277        ig=1;
278        if (gMSAWave[6]>=(1.0+8.0*gMSAWave[4])) {
279                ig=0;
280                gMSAWave[15]=gMSAWave[14];
281                gMSAWave[16]=gMSAWave[4];
282                ix=1;
283                ir = sqfun(ix,ir);
284                gMSAWave[14]=gMSAWave[15];
285                gMSAWave[4]=gMSAWave[16];
286                if((ir<0.0) || (gMSAWave[14]>=0.0)) {
287                        return ir;
288                }
289        }
290        gMSAWave[10]=fmin(gMSAWave[4],0.20);
291        if ((ig!=1) || ( gMSAWave[9]>=0.15)) {
292                ii=0;
293                do {
294                        ii=ii+1;
295                        if(ii>itm) {
296                                ir=-1;
297                                return ir;
298                        }
299                        if (gMSAWave[10]<=0.0) {
300                            gMSAWave[10]=gMSAWave[4]/ii;
301                        }
302                        if(gMSAWave[10]>0.6) {
303                            gMSAWave[10] = 0.35/ii;
304                        }
305                        e1=gMSAWave[10];
306                        gMSAWave[15]=f1;
307                        gMSAWave[16]=e1;
308                        ix=2;
309                        ir = sqfun(ix,ir);
310                        f1=gMSAWave[15];
311                        e1=gMSAWave[16];
312                        e2=gMSAWave[10]*1.01;
313                        gMSAWave[15]=f2;
314                        gMSAWave[16]=e2;
315                        ix=2;
316                        ir = sqfun(ix,ir);
317                        f2=gMSAWave[15];
318                        e2=gMSAWave[16];
319                        e2=e1-(e2-e1)*f1/(f2-f1);
320                        gMSAWave[10] = e2;
321                        del = fabs((e2-e1)/e1);
322                } while (del>acc);
323                gMSAWave[15]=gMSAWave[14];
324                gMSAWave[16]=e2;
325                ix=4;
326                ir = sqfun(ix,ir);
327                gMSAWave[14]=gMSAWave[15];
328                e2=gMSAWave[16];
329                ir=ii;
330                if ((ig!=1) || (gMSAWave[10]>=gMSAWave[4])) {
331                    return ir;
332                }
333        }
334        gMSAWave[15]=gMSAWave[14];
335        gMSAWave[16]=gMSAWave[4];
336        ix=3;
337        ir = sqfun(ix,ir);
338        gMSAWave[14]=gMSAWave[15];
339        gMSAWave[4]=gMSAWave[16];
340        if ((ir>=0) && (gMSAWave[14]<0.0)) {
341                ir=-3;
342        }
343        return ir;
344}
345
346
347int
348sqfun(int ix, int ir)
349{
350        double acc=1.0e-6;
351        double reta,eta2,eta21,eta22,eta3,eta32,eta2d,eta2d2,eta3d,eta6d,e12,e24,rgek;
352        double rak,ak1,ak2,dak,dak2,dak4,d,d2,dd2,dd4,dd45,ex1,ex2,sk,ck,ckma,skma;
353        double al1,al2,al3,al4,al5,al6,be1,be2,be3,vu1,vu2,vu3,vu4,vu5,ph1,ph2,ta1,ta2,ta3,ta4,ta5;
354        double a1,a2,a3,b1,b2,b3,v1,v2,v3,p1,p2,p3,pp,pp1,pp2,p1p2,t1,t2,t3,um1,um2,um3,um4,um5,um6;
355        double w0,w1,w2,w3,w4,w12,w13,w14,w15,w16,w24,w25,w26,w32,w34,w3425,w35,w3526,w36,w46,w56;
356        double fa,fap,ca,e24g,pwk,qpw,pg,del,fun,fund,g24;
357        int ii,ibig,itm=40;
358        //      WAVE gMSAWave = $"root:HayPenMSA:gMSAWave"
359        a2=0;
360        a3=0;
361        b2=0;
362        b3=0;
363        v2=0;
364        v3=0;
365        p2=0;
366        p3=0;
367
368        //     CALCULATE CONSTANTS; NOTATION IS HAYTER PENFOLD (1981)
369
370        reta = gMSAWave[16];
371        eta2 = reta*reta;
372        eta3 = eta2*reta;
373        e12 = 12.0*reta;
374        e24 = e12+e12;
375        gMSAWave[13] = pow( (gMSAWave[4]/gMSAWave[16]),(1.0/3.0));
376        gMSAWave[12]=gMSAWave[6]/gMSAWave[13];
377        ibig=0;
378        if (( gMSAWave[12]>15.0) && (ix==1)) {
379                ibig=1;
380        }
381
382        gMSAWave[11] = gMSAWave[5]*gMSAWave[13]*exp(gMSAWave[6]- gMSAWave[12]);
383        rgek =  gMSAWave[11];
384        rak =  gMSAWave[12];
385        ak2 = rak*rak;
386        ak1 = 1.0+rak;
387        dak2 = 1.0/ak2;
388        dak4 = dak2*dak2;
389        d = 1.0-reta;
390        d2 = d*d;
391        dak = d/rak;
392        dd2 = 1.0/d2;
393        dd4 = dd2*dd2;
394        dd45 = dd4*2.0e-1;
395        eta3d=3.0*reta;
396        eta6d = eta3d+eta3d;
397        eta32 = eta3+ eta3;
398        eta2d = reta+2.0;
399        eta2d2 = eta2d*eta2d;
400        eta21 = 2.0*reta+1.0;
401        eta22 = eta21*eta21;
402
403        //     ALPHA(I)
404
405        al1 = -eta21*dak;
406        al2 = (14.0*eta2-4.0*reta-1.0)*dak2;
407        al3 = 36.0*eta2*dak4;
408
409        //      BETA(I)
410
411        be1 = -(eta2+7.0*reta+1.0)*dak;
412        be2 = 9.0*reta*(eta2+4.0*reta-2.0)*dak2;
413        be3 = 12.0*reta*(2.0*eta2+8.0*reta-1.0)*dak4;
414
415        //      NU(I)
416
417        vu1 = -(eta3+3.0*eta2+45.0*reta+5.0)*dak;
418        vu2 = (eta32+3.0*eta2+42.0*reta-2.0e1)*dak2;
419        vu3 = (eta32+3.0e1*reta-5.0)*dak4;
420        vu4 = vu1+e24*rak*vu3;
421        vu5 = eta6d*(vu2+4.0*vu3);
422
423        //      PHI(I)
424
425        ph1 = eta6d/rak;
426        ph2 = d-e12*dak2;
427
428        //      TAU(I)
429
430        ta1 = (reta+5.0)/(5.0*rak);
431        ta2 = eta2d*dak2;
432        ta3 = -e12*rgek*(ta1+ta2);
433        ta4 = eta3d*ak2*(ta1*ta1-ta2*ta2);
434        ta5 = eta3d*(reta+8.0)*1.0e-1-2.0*eta22*dak2;
435
436        //     double PRECISION SINH(K), COSH(K)
437
438        ex1 = exp(rak);
439        ex2 = 0.0;
440        if ( gMSAWave[12]<20.0) {
441                ex2=exp(-rak);
442        }
443        sk=0.5*(ex1-ex2);
444        ck = 0.5*(ex1+ex2);
445        ckma = ck-1.0-rak*sk;
446        skma = sk-rak*ck;
447
448        //      a(I)
449
450        a1 = (e24*rgek*(al1+al2+ak1*al3)-eta22)*dd4;
451        if (ibig==0) {
452                a2 = e24*(al3*skma+al2*sk-al1*ck)*dd4;
453                a3 = e24*(eta22*dak2-0.5*d2+al3*ckma-al1*sk+al2*ck)*dd4;
454        }
455
456        //      b(I)
457
458        b1 = (1.5*reta*eta2d2-e12*rgek*(be1+be2+ak1*be3))*dd4;
459        if (ibig==0) {
460                b2 = e12*(-be3*skma-be2*sk+be1*ck)*dd4;
461                b3 = e12*(0.5*d2*eta2d-eta3d*eta2d2*dak2-be3*ckma+be1*sk-be2*ck)*dd4;
462        }
463
464        //      V(I)
465
466        v1 = (eta21*(eta2-2.0*reta+1.0e1)*2.5e-1-rgek*(vu4+vu5))*dd45;
467        if (ibig==0) {
468                v2 = (vu4*ck-vu5*sk)*dd45;
469                v3 = ((eta3-6.0*eta2+5.0)*d-eta6d*(2.0*eta3-3.0*eta2+18.0*reta+1.0e1)*dak2+e24*vu3+vu4*sk-vu5*ck)*dd45;
470        }
471
472
473        //       P(I)
474
475        pp1 = ph1*ph1;
476        pp2 = ph2*ph2;
477        pp = pp1+pp2;
478        p1p2 = ph1*ph2*2.0;
479        p1 = (rgek*(pp1+pp2-p1p2)-0.5*eta2d)*dd2;
480        if (ibig==0) {
481                p2 = (pp*sk+p1p2*ck)*dd2;
482                p3 = (pp*ck+p1p2*sk+pp1-pp2)*dd2;
483        }
484
485        //       T(I)
486
487        t1 = ta3+ta4*a1+ta5*b1;
488        if (ibig!=0) {
489
490                //              VERY LARGE SCREENING:  ASYMPTOTIC SOLUTION
491
492                v3 = ((eta3-6.0*eta2+5.0)*d-eta6d*(2.0*eta3-3.0*eta2+18.0*reta+1.0e1)*dak2+e24*vu3)*dd45;
493                t3 = ta4*a3+ta5*b3+e12*ta2 - 4.0e-1*reta*(reta+1.0e1)-1.0;
494                p3 = (pp1-pp2)*dd2;
495                b3 = e12*(0.5*d2*eta2d-eta3d*eta2d2*dak2+be3)*dd4;
496                a3 = e24*(eta22*dak2-0.5*d2-al3)*dd4;
497                um6 = t3*a3-e12*v3*v3;
498                um5 = t1*a3+a1*t3-e24*v1*v3;
499                um4 = t1*a1-e12*v1*v1;
500                al6 = e12*p3*p3;
501                al5 = e24*p1*p3-b3-b3-ak2;
502                al4 = e12*p1*p1-b1-b1;
503                w56 = um5*al6-al5*um6;
504                w46 = um4*al6-al4*um6;
505                fa = -w46/w56;
506                ca = -fa;
507                gMSAWave[3] = fa;
508                gMSAWave[2] = ca;
509                gMSAWave[1] = b1+b3*fa;
510                gMSAWave[0] = a1+a3*fa;
511                gMSAWave[8] = v1+v3*fa;
512                gMSAWave[14] = -(p1+p3*fa);
513                gMSAWave[15] = gMSAWave[14];
514                if (fabs(gMSAWave[15])<1.0e-3) {
515                        gMSAWave[15] = 0.0;
516                }
517                gMSAWave[10] = gMSAWave[16];
518
519        } else {
520
521                t2 = ta4*a2+ta5*b2+e12*(ta1*ck-ta2*sk);
522                t3 = ta4*a3+ta5*b3+e12*(ta1*sk-ta2*(ck-1.0))-4.0e-1*reta*(reta+1.0e1)-1.0;
523
524                //              MU(i)
525
526                um1 = t2*a2-e12*v2*v2;
527                um2 = t1*a2+t2*a1-e24*v1*v2;
528                um3 = t2*a3+t3*a2-e24*v2*v3;
529                um4 = t1*a1-e12*v1*v1;
530                um5 = t1*a3+t3*a1-e24*v1*v3;
531                um6 = t3*a3-e12*v3*v3;
532
533                //                      GILLAN CONDITION ?
534                //
535                //                      YES - G(X=1+) = 0
536                //
537                //                      COEFFICIENTS AND FUNCTION VALUE
538                //
539                if ((ix==1) || (ix==3)) {
540
541                        //                      NO - CALCULATE REMAINING COEFFICIENTS.
542
543                        //                      LAMBDA(I)
544
545                        al1 = e12*p2*p2;
546                        al2 = e24*p1*p2-b2-b2;
547                        al3 = e24*p2*p3;
548                        al4 = e12*p1*p1-b1-b1;
549                        al5 = e24*p1*p3-b3-b3-ak2;
550                        al6 = e12*p3*p3;
551
552                        //                      OMEGA(I)
553
554                        w16 = um1*al6-al1*um6;
555                        w15 = um1*al5-al1*um5;
556                        w14 = um1*al4-al1*um4;
557                        w13 = um1*al3-al1*um3;
558                        w12 = um1*al2-al1*um2;
559
560                        w26 = um2*al6-al2*um6;
561                        w25 = um2*al5-al2*um5;
562                        w24 = um2*al4-al2*um4;
563
564                        w36 = um3*al6-al3*um6;
565                        w35 = um3*al5-al3*um5;
566                        w34 = um3*al4-al3*um4;
567                        w32 = um3*al2-al3*um2;
568                        //
569                        w46 = um4*al6-al4*um6;
570                        w56 = um5*al6-al5*um6;
571                        w3526 = w35+w26;
572                        w3425 = w34+w25;
573
574                        //                      QUARTIC COEFFICIENTS
575
576                        w4 = w16*w16-w13*w36;
577                        w3 = 2.0*w16*w15-w13*w3526-w12*w36;
578                        w2 = w15*w15+2.0*w16*w14-w13*w3425-w12*w3526;
579                        w1 = 2.0*w15*w14-w13*w24-w12*w3425;
580                        w0 = w14*w14-w12*w24;
581
582                        //                      ESTIMATE THE STARTING VALUE OF f
583
584                        if (ix==1) {
585                                //                              LARGE K
586                                fap = (w14-w34-w46)/(w12-w15+w35-w26+w56-w32);
587                        } else {
588                                //                              ASSUME NOT TOO FAR FROM GILLAN CONDITION.
589                                //                              IF BOTH RGEK AND RAK ARE SMALL, USE P-W ESTIMATE.
590                                gMSAWave[14]=0.5*eta2d*dd2*exp(-rgek);
591                                if (( gMSAWave[11]<=2.0) && ( gMSAWave[11]>=0.0) && ( gMSAWave[12]<=1.0)) {
592                                        e24g = e24*rgek*exp(rak);
593                                        pwk = sqrt(e24g);
594                                        qpw = (1.0-sqrt(1.0+2.0*d2*d*pwk/eta22))*eta21/d;
595                                        gMSAWave[14] = -qpw*qpw/e24+0.5*eta2d*dd2;
596                                }
597                                pg = p1+gMSAWave[14];
598                                ca = ak2*pg+2.0*(b3*pg-b1*p3)+e12*gMSAWave[14]*gMSAWave[14]*p3;
599                                ca = -ca/(ak2*p2+2.0*(b3*p2-b2*p3));
600                                fap = -(pg+p2*ca)/p3;
601                        }
602
603                        //                      AND REFINE IT ACCORDING TO NEWTON
604                        ii=0;
605                        do {
606                                ii = ii+1;
607                                if (ii>itm) {
608                                        //                                      FAILED TO CONVERGE IN ITM ITERATIONS
609                                        ir=-2;
610                                        return (ir);
611                                }
612                                fa = fap;
613                                fun = w0+(w1+(w2+(w3+w4*fa)*fa)*fa)*fa;
614                                fund = w1+(2.0*w2+(3.0*w3+4.0*w4*fa)*fa)*fa;
615                                fap = fa-fun/fund;
616                                del=fabs((fap-fa)/fa);
617                        } while (del>acc);
618
619                        ir = ir+ii;
620                        fa = fap;
621                        ca = -(w16*fa*fa+w15*fa+w14)/(w13*fa+w12);
622                        gMSAWave[14] = -(p1+p2*ca+p3*fa);
623                        gMSAWave[15] = gMSAWave[14];
624                        if (fabs(gMSAWave[15])<1.0e-3) {
625                                gMSAWave[15] = 0.0;
626                        }
627                        gMSAWave[10] = gMSAWave[16];
628                } else {
629                        ca = ak2*p1+2.0*(b3*p1-b1*p3);
630                        ca = -ca/(ak2*p2+2.0*(b3*p2-b2*p3));
631                        fa = -(p1+p2*ca)/p3;
632                        if (ix==2) {
633                                gMSAWave[15] = um1*ca*ca+(um2+um3*fa)*ca+um4+um5*fa+um6*fa*fa;
634                        }
635                        if (ix==4) {
636                                gMSAWave[15] = -(p1+p2*ca+p3*fa);
637                        }
638                }
639                gMSAWave[3] = fa;
640                gMSAWave[2] = ca;
641                gMSAWave[1] = b1+b2*ca+b3*fa;
642                gMSAWave[0] = a1+a2*ca+a3*fa;
643                gMSAWave[8] = (v1+v2*ca+v3*fa)/gMSAWave[0];
644        }
645        g24 = e24*rgek*ex1;
646        gMSAWave[7] = (rak*ak2*ca-g24)/(ak2*g24);
647        return (ir);
648}
649
650// called as DiamCyl(hcyl,rcyl)
651//modified from Igor NIST package XOP
652double
653DiamCyl(double hcyl, double rcyl)
654{
655        double diam,a,b,t1,t2,ddd;
656        double pi;
657
658        pi = 4.0*atan(1.0);
659        if (rcyl == 0 || hcyl == 0) {
660                return 0.0;
661        }
662        a = rcyl;
663        b = hcyl/2.0;
664        t1 = a*a*2.0*b/2.0;
665        t2 = 1.0 + (b/a)*(1.0+a/b/2)*(1.0+pi*a/b/2.0);
666        ddd = 3.0*t1*t2;
667        diam = pow(ddd,(1.0/3.0));
668
669        return(diam);
670}
671
672//prolate OR oblate ellipsoids
673//aa is the axis of rotation
674//if aa>bb, then PROLATE
675//if aa<bb, then OBLATE
676// A. Isihara, J. Chem. Phys. 18, 1446 (1950)
677//returns DIAMETER
678// called as DiamEllip(aa,bb)
679
680//modified from Igor NIST package XOP
681double
682DiamEllip(double aa, double bb)
683{
684        double ee,e1,bd,b1,bL,b2,del,ddd,diam;
685
686        if (aa == 0 || bb == 0) {
687                return 0.0;
688        }
689        if (aa == bb) {
690                return 2*aa;
691        }
692        if(aa>bb) {
693                ee = (aa*aa - bb*bb)/(aa*aa);
694        }else{
695                ee = (bb*bb - aa*aa)/(bb*bb);
696        }
697
698        bd = 1.0-ee;
699        e1 = sqrt(ee);
700        b1 = 1.0 + asin(e1)/(e1*sqrt(bd));
701        bL = (1.0+e1)/(1.0-e1);
702        b2 = 1.0 + bd/2/e1*log(bL);
703        del = 0.75*b1*b2;
704
705        ddd = 2.0*(del+1.0)*aa*bb*bb;           //volume is always calculated correctly
706        diam = pow(ddd,(1.0/3.0));
707
708        return(diam);
709}
710
711double
712sqhcal(double qq)
713{
714    double SofQ,etaz,akz,gekz,e24,x1,x2,ck,sk,ak2,qk,q2k,qk2,qk3,qqk,sink,cosk,asink,qcosk,aqk,inter;
715        //      WAVE gMSAWave = $"root:HayPenMSA:gMSAWave"
716
717        etaz = gMSAWave[10];
718        akz =  gMSAWave[12];
719        gekz =  gMSAWave[11];
720        e24 = 24.0*etaz;
721        x1 = exp(akz);
722        x2 = 0.0;
723        if ( gMSAWave[12]<20.0) {
724                x2 = exp(-akz);
725        }
726        ck = 0.5*(x1+x2);
727        sk = 0.5*(x1-x2);
728        ak2 = akz*akz;
729
730        if (qq<=0.0) {
731                SofQ = -1.0/gMSAWave[0];
732        } else {
733                qk = qq/gMSAWave[13];
734                q2k = qk*qk;
735                qk2 = 1.0/q2k;
736                qk3 = qk2/qk;
737                qqk = 1.0/(qk*(q2k+ak2));
738                sink = sin(qk);
739                cosk = cos(qk);
740                asink = akz*sink;
741                qcosk = qk*cosk;
742                aqk = gMSAWave[0]*(sink-qcosk);
743                aqk=aqk+gMSAWave[1]*((2.0*qk2-1.0)*qcosk+2.0*sink-2.0/qk);
744                inter=24.0*qk3+4.0*(1.0-6.0*qk2)*sink;
745                aqk=(aqk+0.5*etaz*gMSAWave[0]*(inter-(1.0-12.0*qk2+24.0*qk2*qk2)*qcosk))*qk3;
746                aqk=aqk +gMSAWave[2]*(ck*asink-sk*qcosk)*qqk;
747                aqk=aqk +gMSAWave[3]*(sk*asink-qk*(ck*cosk-1.0))*qqk;
748                aqk=aqk +gMSAWave[3]*(cosk-1.0)*qk2;
749                aqk=aqk -gekz*(asink+qcosk)*qqk;
750                SofQ = 1.0/(1.0-e24*aqk);
751        }
752        return (SofQ);
753}
754
755///////////end of XOP
756
757
Note: See TracBrowser for help on using the repository browser.