source: sasview/sansmodels/src/libigor/libStructureFactor.c @ 560c4d8

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 560c4d8 was 67424cd, checked in by Mathieu Doucet <doucetm@…>, 13 years ago

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