source: sasview/sansmodels/src/libigor/libStructureFactor.c @ 384647dc

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 384647dc was 6e93a02, checked in by Jae Cho <jhjcho@…>, 15 years ago

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