source: sans/Analysis/branches/ajj_23APR07/XOPs/SANSAnalysis/StructureFactor.c @ 93

Last change on this file since 93 was 93, checked in by ajj, 15 years ago

Add combined XOP code. Currently only contains XCode project file to build Universal binary suitable for Igor 6.

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