# source:sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/HPMSA_v40.ipf@629

Last change on this file since 629 was 629, checked in by srkline, 13 years ago

1 - added function in HPMSA to calculate U(r). See the function for details of usage.

2 - fixed bug in NSORT to trim the _res wave at the point when it is loaded.

File size: 24.4 KB
Line
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.1
3
4// converted from John Hayter's Fortran to Igor code by Paul Butler
5//
6Proc PlotHayterPenfoldMSA(num,qmin,qmax)
7        Variable num=128, qmin=.001, qmax=.3
8        Prompt num "Enter number of data points for model: "
9        Prompt qmin "Enter minimum q-value (A^-1) for model: "
10        Prompt qmin "Enter minimum q-value (A^-1) for model: "
11
12        //      Set up data folder for global variables
13        if(!DataFolderExists(":HayPenMSA"))
14                NewDataFolder :HayPenMSA
15        endif
16        Make/O/D/N=17 :HayPenMSA:gMSAWave
17        //
18        make/o/d/n=(num) xwave_hpmsa, ywave_hpmsa, xdiamwave_hpmsa
19        xwave_hpmsa = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
20        make/o/d coef_hpmsa = {41.5,19,0.0192,298,0.0,78}                    //**** numerical vals, # of variables
21        make/o/t parameters_hpmsa = {"Diameter (A)","Charge","Volume Fraction","Temperature(K)","monovalent salt conc. (M)","dielectric constant of solvent"}
22        Edit parameters_hpmsa,coef_hpmsa
23        Variable/G root:g_hpmsa
24        g_hpmsa := HayterPenfoldMSA(coef_hpmsa,ywave_hpmsa,xwave_hpmsa)
25//      ywave_hpmsa := HayterPenfoldMSA(coef_hpmsa,xwave_hpmsa)
26        xdiamwave_hpmsa:=xwave_hpmsa*coef_hpmsa[0]
27        //Display ywave_hpmsa vs xdiamwave_hpmsa
28        Display ywave_hpmsa vs xwave_hpmsa
29        ModifyGraph log=0,marker=29,msize=2,mode=4,grid=1                       //**** log=0 if linear scale desired
30        Label bottom "q (A\\S-1\\M)"
31        Label left "Structure Factor"
32        AutoPositionWindow/M=1/R=\$(WinName(0,1)) \$WinName(0,2)
33
35End
36
37//AAO version
38Function HayterPenfoldMSA(cw,yw,xw) : FitFunc
39        Wave cw,yw,xw
40
41#if exists("HayterPenfoldMSAX")
42        yw = HayterPenfoldMSAX(cw,xw)
43#else
44        yw = fHayterPenfoldMSA(cw,xw)
45#endif
46        return(0)
47End
48
49//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
50//
51//   Wrapper for Hayter Penfold MSA routines:
52//              SETS UP THE PARAMETERS FOR THE
53//              CALCULATION OF THE STRUCTURE FACTOR ,S(Q)
54//              GIVEN THE THREE REQUIRED PARAMETERS VALK, GAMMA, ETA.
55//
56//      *** NOTE ****  THIS CALCULATION REQUIRES THAT THE NUMBER OF
57//                     Q-VALUES AT WHICH THE S(Q) IS CALCULATED BE
58//                     A POWER OF 2
59//
60Function fHayterPenfoldMSA(w,x) : FitFunc
61        wave w
62        variable x
63
64//      variable timer=StartMSTimer
65
66        variable Elcharge=1.602189e-19          // electron charge in Coulombs (C)
67        variable kB=1.380662e-23                                // Boltzman constant in J/K
68        variable FrSpPerm=8.85418782E-12        //Permittivity of free space in C^2/(N m^2)
69
70        variable SofQ, QQ, Qdiam, gammaek, Vp, csalt, ss
71        variable VolFrac, SIdiam, diam, Kappa, cs, IonSt
72        variable dialec, Perm, Beta, Temp, zz, charge, ierr
73
74        diam=w[0]               //in A  (not SI .. should force people to think in nm!!!)
75        zz = w[1]               //# of charges
76        VolFrac=w[2]
77        QQ=x                    //in A^-1 (not SI .. should force people to think in nm^-1!!!)
78        Temp=w[3]               //in degrees Kelvin
79        csalt=w[4]              //in molarity
80        dialec=w[5]             // unitless
81
82//                       Set up data folder for global variables in the calling macro
83//      string SaveDF=GetDataFolder(1)
84 //     if (DataFolderExists("root:HayPenMSA"))
85//              SetDataFolder root:HayPenMSA
86 //     else
87//              NewDataFolder/S root:HayPenMSA
88//      endif
89//      variable/G a,b,c,f,eta,gek,ak,u,v,gamk,seta,sgek,sak,scal,g1, fval, evar
90//      SetDataFolder SaveDF
91
92//need NVAR references to ALL global variables in HayPenMSA subfolder
93//      NVAR a=root:HayPenMSA:a
94//      NVAR b=root:HayPenMSA:b
95//      NVAR c=root:HayPenMSA:c
96//      NVAR f=root:HayPenMSA:f
97//      NVAR eta=root:HayPenMSA:eta
98//      NVAR gek=root:HayPenMSA:gek
99//      NVAR ak=root:HayPenMSA:ak
100//      NVAR u=root:HayPenMSA:u
101//      NVAR v=root:HayPenMSA:v
102//      NVAR gamk=root:HayPenMSA:gamk
103//      NVAR seta=root:HayPenMSA:seta
104//      NVAR sgek=root:HayPenMSA:sgek
105//      NVAR sak=root:HayPenMSA:sak
106//      NVAR scal=root:HayPenMSA:scal
107//      NVAR g1=root:HayPenMSA:g1
108//      NVAR fval=root:HayPenMSA:fval
109//      NVAR evar=root:HayPenMSA:evar
110
111//      NVAR a=root:HayPenMSA:a = gMSAWave[0]
112//      NVAR b=root:HayPenMSA:b = gMSAWave[1]
113//      NVAR c=root:HayPenMSA:c = gMSAWave[2]
114//      NVAR f=root:HayPenMSA:f = gMSAWave[3]
115//      NVAR eta=root:HayPenMSA:eta = gMSAWave[4]
116//      NVAR gek=root:HayPenMSA:gek = gMSAWave[5]
117//      NVAR ak=root:HayPenMSA:ak = gMSAWave[6]
118//      NVAR u=root:HayPenMSA:u = gMSAWave[7]
119//      NVAR v=root:HayPenMSA:v = gMSAWave[8]
120//      NVAR gamk=root:HayPenMSA:gamk = gMSAWave[9]
121//      NVAR seta=root:HayPenMSA:seta = gMSAWave[10]
122//      NVAR sgek=root:HayPenMSA:sgek = gMSAWave[11]
123//      NVAR sak=root:HayPenMSA:sak = gMSAWave[12]
124//      NVAR scal=root:HayPenMSA:scal = gMSAWave[13]
125//      NVAR g1=root:HayPenMSA:g1 = gMSAWave[14]
126//      NVAR fval=root:HayPenMSA:fval = gMSAWave[15]
127//      NVAR evar=root:HayPenMSA:evar = gMSAWave[16]
129        WAVE gMSAWave = \$":HayPenMSA:gMSAWave"
130
131
132////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
133//////////////////////////// convert to USEFUL inputs in SI units                                                //
134////////////////////////////    NOTE: easiest to do EVERYTHING in SI units                               //
135////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
136        Beta=1/(kB*Temp)                // in Joules^-1
137        Perm=dialec*FrSpPerm    //in C^2/(N  m^2)
138        charge=zz*Elcharge              //in Coulomb (C)
139        SIdiam = diam*1E-10             //in m
140        Vp=4*pi/3*(SIdiam/2)^3  //in m^3
141        cs=csalt*6.022E23*1E3   //# salt molecules/m^3
142
143//         Compute the derived values of :
144//                       Ionic strength IonSt (in C^2/m^3)
145//                      Kappa (Debye-Huckel screening length in m)
146//      and             gamma Exp(-k)
147        IonSt=0.5 * Elcharge^2*(zz*VolFrac/Vp+2*cs)
148        Kappa=sqrt(2*Beta*IonSt/Perm)     //Kappa calc from Ionic strength
149//      Kappa=2/SIdiam                                  // Use to compare with HP paper
150        gMSAWave[5]=Beta*charge^2/(pi*Perm*SIdiam*(2+Kappa*SIdiam)^2)
151
152//         Finally set up dimensionless parameters
153        Qdiam=QQ*diam
154      gMSAWave[6] = Kappa*SIdiam
155      gMSAWave[4] = VolFrac
156
157
158
159//  ***************  now go to John Hayter and Jeff Penfold setup routine************
160
161
162//    *** ALL FURTHER PROGRAMS COMMENTS ARE FROM J. HAYTER
163//        EXCEPT WHERE INDICATED  ^*
164//
165//
166//       ROUTINE TO CALCULATE S(Q*SIG) FOR A SCREENED COULOMB
167//       POTENTIAL BETWEEN FINITE PARTICLES OF DIAMETER 'SIG'
168//       AT ANY VOLUME FRACTION.  THIS ROUTINE IS MUCH MORE POWER-
169//       FUL THAN "SQHP" AND SHOULD BE USED TO REPLACE THE LATTER
170//       IN EXISTING PROGRAMS.  NOTE THAT THE COMMON AREA IS
171//       CHANGED;  IN PARTICULAR THE POTENTIAL IS PASSED
172//       DIRECTLY AS 'GEK' = GAMMA*EXP(-K) IN THE PRESENT ROUTINE.
173//
174//     JOHN B. HAYTER
175//
176//  ***** THIS VERSION ENTERED ON 5/30/85 BY JOHN F. BILLMAN
177//
178//       CALLING SEQUENCE:
179//
180//            CALL SQHPA(QQ,SQ,NPT,IERR)
181//
182//      QQ:   ARRAY OF DIMENSION NPT CONTAINING THE VALUES
183//            OF Q*SIG AT WHICH S(Q*SIG) WILL BE CALCULATED.
184//      SQ:   ARRAY OF DIMENSION NPT INTO WHICH THE VALUES OF
185//            S(Q*SIG) WILL BE RETURNED
186//      NPT:  NUMBER OF VALUES OF Q*SIG
187//
188//      IERR  > 0:   NORMAL EXIT; IERR=NUMBER OF ITERATIONS
189//             -1:   NEWTON ITERATION NON-CONVERGENT IN "SQCOEF"
190//             -2:   NEWTON ITERATION NON-CONVERGENT IN "SQFUN"
191//             -3:   CANNOT RESCALE TO G(1+) > 0.
192//
193//        ALL OTHER PARAMETERS ARE TRANSMITTED THROUGH A SINGLE
194//        NAMED COMMON AREA:
195//
196//     REAL*8 a,b,//,f
197//     COMMON /SQHPB/ ETA,GEK,AK,A,B,C,F,U,V,GAMK,SETA,SGEK,SAK,SCAL,G1
198//
199//     ON ENTRY:
200//
201//       ETA:    VOLUME FRACTION
202//       GEK:    THE CONTACT POTENTIAL GAMMA*EXP(-K)
203//        AK:    THE DIMENSIONLESS SCREENING CONSTANT
204//               K=KAPPA*SIG  WHERE KAPPA IS THE INVERSE SCREENING
205//               LENGTH AND SIG IS THE PARTICLE DIAMETER
206//
207//     ON EXIT:
208//
209//       GAMK IS THE COUPLING:  2*GAMMA*SS*EXP(-K/SS), SS=ETA^(1/3).
210//       SETA,SGEK AND SAK ARE THE RESCALED INPUT PARAMETERS.
211//       SCAL IS THE RESCALING FACTOR:  (ETA/SETA)^(1/3).
212//       G1=G(1+), THE CONTACT VALUE OF G(R/SIG).
213//       A.B,C,F,U,V ARE THE CONSTANTS APPEARING IN THE ANALYTIC
214//       SOLUTION OF THE MSA [HAYTER-PENFOLD; MOL. PHYS. 42: 109 (1981)
215//
216//     NOTES:
217//
218//       (A)  AFTER THE FIRST CALL TO SQHPA, S(Q*SIG) MAY BE EVALUATED
219//            AT OTHER Q*SIG VALUES BY REDEFINING THE ARRAY QQ AND CALLING
220//            "SQHCAL" DIRECTLY FROM THE MAIN PROGRAM.
221//
222//       (B)  THE RESULTING S(Q*SIG) MAY BE TRANSFORMED TO G(SS/SIG)
223//            BY USING THE ROUTINE "TROGS"
224//
225//       (C)  NO ERROR CHECKING OF INPUT PARAMETERS IS PERFORMED;
226//            IT IS THE RESPONSIBILITY OF THE CALLING PROGRAM TO
227//            VERIFY VALIDITY.
228//
229//      SUBROUTINES CALLED BY SQHPA:
230//
231//         (1) SQCOEF:    RESCALES THE PROBLEM AND CALCULATES THE
232//                        APPROPRIATE COEFFICIENTS FOR "SQHCAL"
233//
234//         (2) SQFUN:     CALCULATES VARIOUS VALUES FOR "SQCOEF"
235//
236//         (3) SQHCAL:    CALCULATES H-P  S(Q*SIG)  GIVEN IN A,B,C,F
237//
238
239
240
241
242//Function sqhpa(qq)  {this is where Hayter-Penfold program began}
243
244
245//       FIRST CALCULATE COUPLING
246
247      ss=gMSAWave[4]^(1.0/3.0)
248       gMSAWave[9] = 2.0*ss*gMSAWave[5]*exp(gMSAWave[6]-gMSAWave[6]/ss)
249
250//        CALCULATE COEFFICIENTS, CHECK ALL IS WELL
251//        AND IF SO CALCULATE S(Q*SIG)
252
253      ierr=0
254      ierr=sqcoef(ierr)
255      if (ierr>=0)
256            SofQ=sqhcal(Qdiam)
257       else
258        SofQ=NaN
259        print "Error Level = ",ierr
260             print "Please report HPMSA problem with above error code"
261      endif
262      //KillDataFolder root:HayPenMSA
263//      variable elapsed=StopMSTimer(timer)
264//      Print "elapsed time = ",elapsed
265
266      return SofQ
267end
268
269/////////////////////////////////////////////////////////////
270/////////////////////////////////////////////////////////////
271//
272//
273//      CALCULATES RESCALED VOLUME FRACTION AND CORRESPONDING
274//      COEFFICIENTS FOR "SQHPA"
275//
276//      JOHN B. HAYTER   (I.L.L.)    14-SEP-81
277//
278//      ON EXIT:
279//
280//      SETA IS THE RESCALED VOLUME FRACTION
281//      SGEK IS THE RESCALED CONTACT POTENTIAL
282//      SAK IS THE RESCALED SCREENING CONSTANT
283//      A,B,C,F,U,V ARE THE MSA COEFFICIENTS
284//      G1= G(1+) IS THE CONTACT VALUE OF G(R/SIG):
285//      FOR THE GILLAN CONDITION, THE DIFFERENCE FROM
286//      ZERO INDICATES THE COMPUTATIONAL ACCURACY.
287//
288//      IR > 0:    NORMAL EXIT,  IR IS THE NUMBER OF ITERATIONS.
289//         < 0:    FAILED TO CONVERGE
290//
291
292
293Function sqcoef(ir)
294        variable ir
295
296      variable itm=40.0, acc=5.0E-6, ix,ig,ii,del,e1,e2,f1,f2
297      WAVE gMSAWave = \$":HayPenMSA:gMSAWave"
298
299//      NVAR a=root:HayPenMSA:a
300//      NVAR b=root:HayPenMSA:b
301//      NVAR c=root:HayPenMSA:c
302//      NVAR f=root:HayPenMSA:f
303//      NVAR eta=root:HayPenMSA:eta
304//      NVAR gek=root:HayPenMSA:gek
305//      NVAR ak=root:HayPenMSA:ak
306//      NVAR u=root:HayPenMSA:u
307//      NVAR v=root:HayPenMSA:v
308//      NVAR gamk=root:HayPenMSA:gamk
309//      NVAR seta=root:HayPenMSA:seta
310//      NVAR sgek=root:HayPenMSA:sgek
311//      NVAR sak=root:HayPenMSA:sak
312//      NVAR scal=root:HayPenMSA:scal
313//      NVAR g1=root:HayPenMSA:g1
314//      NVAR fval=root:HayPenMSA:fval
315//      NVAR evar=root:HayPenMSA:evar
316
317
318      ig=1
319      if (gMSAWave[6]>=(1.0+8.0*gMSAWave[4]))
320                ig=0
321                gMSAWave[15]=gMSAWave[14]
322                 gMSAWave[16]=gMSAWave[4]
323                  ix=1
324                ir = sqfun(ix,ir)
325                gMSAWave[14]=gMSAWave[15]
326                 gMSAWave[4]=gMSAWave[16]
327                if((ir<0.0) %| (gMSAWave[14]>=0.0))
328                   return ir
329                endif
330        endif
331      gMSAWave[10]=min(gMSAWave[4],0.20)
332      if ((ig!=1) %| ( gMSAWave[9]>=0.15))
333                ii=0
334                do
335                        ii=ii+1
336                        if(ii>itm)
337                                ir=-1
338                                return ir
339                        endif
340                        if (gMSAWave[10]<=0.0)
341                            gMSAWave[10]=gMSAWave[4]/ii
342                        endif
343                        if(gMSAWave[10]>0.6)
344                            gMSAWave[10] = 0.35/ii
345                         endif
346                        e1=gMSAWave[10]
347                        gMSAWave[15]=f1
348                        gMSAWave[16]=e1
349                        ix=2
350                        ir = sqfun(ix,ir)
351                        f1=gMSAWave[15]
352                         e1=gMSAWave[16]
353                        e2=gMSAWave[10]*1.01
354                        gMSAWave[15]=f2
355                         gMSAWave[16]=e2
356                         ix=2
357                        ir = sqfun(ix,ir)
358                        f2=gMSAWave[15]
359                         e2=gMSAWave[16]
360                        e2=e1-(e2-e1)*f1/(f2-f1)
361                        gMSAWave[10] = e2
362                        del = abs((e2-e1)/e1)
363                while (del>acc)
364                gMSAWave[15]=gMSAWave[14]
365                 gMSAWave[16]=e2
366                 ix=4
367                ir = sqfun(ix,ir)
368                gMSAWave[14]=gMSAWave[15]
369                 e2=gMSAWave[16]
370                ir=ii
371                if ((ig!=1) %| (gMSAWave[10]>=gMSAWave[4]))
372                    return ir
373                endif
374        endif
375        gMSAWave[15]=gMSAWave[14]
376        gMSAWave[16]=gMSAWave[4]
377         ix=3
378        ir = sqfun(ix,ir)
379        gMSAWave[14]=gMSAWave[15]
380         gMSAWave[4]=gMSAWave[16]
381        if ((ir>=0) %& (gMSAWave[14]<0.0))
382              ir=-3
383        endif
384        return ir
385end
386
387///////////////////////////////////////////////////////
388///////////////////////////////////////////////////////
389//
390//
391//       CALCULATES VARIOUS COEFFICIENTS AND FUNCTION
392//       VALUES FOR "SQCOEF"  (USED BY "SQHPA").
393//
394//   **  THIS ROUTINE WORKS LOCALLY IN DOUBLE PRECISION  **
395//
396//       JOHN B. HAYTER    (I.L.L.)  23-MAR-82
397//
398//       IX=1:   SOLVE FOR LARGE K, RETURN G(1+).
399//          2:   RETURN FUNCTION TO SOLVE FOR ETA(GILLAN).
400//          3:   ASSUME NEAR GILLAN, SOLVE, RETURN G(1+).
401//          4:   RETURN G(1+) FOR ETA=ETA(GILLAN).
402//
403//
404
405
406
407Function sqfun(ix,ir)
408        variable ix, ir
409
410      variable acc=1.0e-6, itm=40.0
411      variable reta,eta2,eta21,eta22,eta3,eta32,eta2d,eta2d2,eta3d,eta6d,e12,e24,ibig,rgek
412      variable rak,ak1,ak2,dak,dak2,dak4,d,d2,dd2,dd4,dd45,ex1,ex2,sk,ck,ckma,skma
413      variable al1,al2,al3,al4,al5,al6,be1,be2,be3,vu1,vu2,vu3,vu4,vu5,ph1,ph2,ta1,ta2,ta3,ta4,ta5
414      variable 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
415      variable w0,w1,w2,w3,w4,w12,w13,w14,w15,w16,w24,w25,w26,w32,w34,w3425,w35,w3526,w36,w46,w56
416      variable fa,fap,ca,e24g,pwk,qpw,pg,ii,del,fun,fund,g24
417      WAVE gMSAWave = \$":HayPenMSA:gMSAWave"
418
419//      NVAR a=root:HayPenMSA:a
420//      NVAR b=root:HayPenMSA:b
421//      NVAR c=root:HayPenMSA:c
422//      NVAR f=root:HayPenMSA:f
423//      NVAR eta=root:HayPenMSA:eta
424//      NVAR gek=root:HayPenMSA:gek
425//      NVAR ak=root:HayPenMSA:ak
426//      NVAR u=root:HayPenMSA:u
427//      NVAR v=root:HayPenMSA:v
428//      NVAR gamk=root:HayPenMSA:gamk
429//      NVAR seta=root:HayPenMSA:seta
430//      NVAR sgek=root:HayPenMSA:sgek
431//      NVAR sak=root:HayPenMSA:sak
432//      NVAR scal=root:HayPenMSA:scal
433//      NVAR g1=root:HayPenMSA:g1
434//      NVAR fval=root:HayPenMSA:fval
435//      NVAR evar=root:HayPenMSA:evar
436
437
438//     CALCULATE CONSTANTS; NOTATION IS HAYTER PENFOLD (1981)
439
440      reta = gMSAWave[16]
441      eta2 = reta*reta
442      eta3 = eta2*reta
443      e12 = 12.0*reta
444      e24 = e12+e12
445      gMSAWave[13] = (gMSAWave[4]/gMSAWave[16])^(1.0/3.0)
446       gMSAWave[12]=gMSAWave[6]/gMSAWave[13]
447      ibig=0
448      if (( gMSAWave[12]>15.0) %& (ix==1))
449            ibig=1
450      endif
451       gMSAWave[11] = gMSAWave[5]*gMSAWave[13]*exp(gMSAWave[6]- gMSAWave[12])
452      rgek =  gMSAWave[11]
453      rak =  gMSAWave[12]
454      ak2 = rak*rak
455      ak1 = 1.0+rak
456      dak2 = 1.0/ak2
457      dak4 = dak2*dak2
458      d = 1.0-reta
459      d2 = d*d
460      dak = d/rak
461      dd2 = 1.0/d2
462      dd4 = dd2*dd2
463      dd45 = dd4*2.0e-1
464      eta3d=3.0*reta
465      eta6d = eta3d+eta3d
466      eta32 = eta3+ eta3
467      eta2d = reta+2.0
468      eta2d2 = eta2d*eta2d
469      eta21 = 2.0*reta+1.0
470      eta22 = eta21*eta21
471
472//     ALPHA(I)
473
474      al1 = -eta21*dak
475      al2 = (14.0*eta2-4.0*reta-1.0)*dak2
476      al3 = 36.0*eta2*dak4
477
478//      BETA(I)
479
480      be1 = -(eta2+7.0*reta+1.0)*dak
481      be2 = 9.0*reta*(eta2+4.0*reta-2.0)*dak2
482      be3 = 12.0*reta*(2.0*eta2+8.0*reta-1.0)*dak4
483
484//      NU(I)
485
486      vu1 = -(eta3+3.0*eta2+45.0*reta+5.0)*dak
487      vu2 = (eta32+3.0*eta2+42.0*reta-2.0e1)*dak2
488      vu3 = (eta32+3.0e1*reta-5.0)*dak4
489      vu4 = vu1+e24*rak*vu3
490      vu5 = eta6d*(vu2+4.0*vu3)
491
492//      PHI(I)
493
494      ph1 = eta6d/rak
495      ph2 = d-e12*dak2
496
497//      TAU(I)
498
499      ta1 = (reta+5.0)/(5.0*rak)
500      ta2 = eta2d*dak2
501      ta3 = -e12*rgek*(ta1+ta2)
502      ta4 = eta3d*ak2*(ta1*ta1-ta2*ta2)
503      ta5 = eta3d*(reta+8.0)*1.0e-1-2.0*eta22*dak2
504
505//     DOUBLE PRECISION SINH(K), COSH(K)
506
507      ex1 = exp(rak)
508      ex2 = 0.0
509      if ( gMSAWave[12]<20.0)
510           ex2=exp(-rak)
511      endif
512      sk=0.5*(ex1-ex2)
513      ck = 0.5*(ex1+ex2)
514      ckma = ck-1.0-rak*sk
515      skma = sk-rak*ck
516
517//      a(I)
518
519      a1 = (e24*rgek*(al1+al2+ak1*al3)-eta22)*dd4
520      if (ibig==0)
521                a2 = e24*(al3*skma+al2*sk-al1*ck)*dd4
522                a3 = e24*(eta22*dak2-0.5*d2+al3*ckma-al1*sk+al2*ck)*dd4
523        endif
524
525//      b(I)
526
527      b1 = (1.5*reta*eta2d2-e12*rgek*(be1+be2+ak1*be3))*dd4
528      if (ibig==0)
529                b2 = e12*(-be3*skma-be2*sk+be1*ck)*dd4
530                b3 = e12*(0.5*d2*eta2d-eta3d*eta2d2*dak2-be3*ckma+be1*sk-be2*ck)*dd4
531        endif
532
533//      V(I)
534
535      v1 = (eta21*(eta2-2.0*reta+1.0e1)*2.5e-1-rgek*(vu4+vu5))*dd45
536      if (ibig==0)
537                v2 = (vu4*ck-vu5*sk)*dd45
538                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
539        endif
540
541
542//       P(I)
543
544      pp1 = ph1*ph1
545      pp2 = ph2*ph2
546      pp = pp1+pp2
547      p1p2 = ph1*ph2*2.0
548      p1 = (rgek*(pp1+pp2-p1p2)-0.5*eta2d)*dd2
549      if (ibig==0)
550                p2 = (pp*sk+p1p2*ck)*dd2
551                p3 = (pp*ck+p1p2*sk+pp1-pp2)*dd2
552        endif
553
554//       T(I)
555
556      t1 = ta3+ta4*a1+ta5*b1
557      if (ibig!=0)
558
559//              VERY LARGE SCREENING:  ASYMPTOTIC SOLUTION
560
561                v3 = ((eta3-6.0*eta2+5.0)*d-eta6d*(2.0*eta3-3.0*eta2+18.0*reta+1.0e1)*dak2+e24*vu3)*dd45
562                t3 = ta4*a3+ta5*b3+e12*ta2 - 4.0e-1*reta*(reta+1.0e1)-1.0
563                p3 = (pp1-pp2)*dd2
564                b3 = e12*(0.5*d2*eta2d-eta3d*eta2d2*dak2+be3)*dd4
565                a3 = e24*(eta22*dak2-0.5*d2-al3)*dd4
566                um6 = t3*a3-e12*v3*v3
567                um5 = t1*a3+a1*t3-e24*v1*v3
568                um4 = t1*a1-e12*v1*v1
569                al6 = e12*p3*p3
570                al5 = e24*p1*p3-b3-b3-ak2
571                al4 = e12*p1*p1-b1-b1
572                w56 = um5*al6-al5*um6
573                w46 = um4*al6-al4*um6
574                fa = -w46/w56
575                ca = -fa
576                gMSAWave[3] = fa
577                gMSAWave[2] = ca
578                gMSAWave[1] = b1+b3*fa
579                gMSAWave[0] = a1+a3*fa
580                gMSAWave[8] = v1+v3*fa
581                gMSAWave[14] = -(p1+p3*fa)
582                gMSAWave[15] = gMSAWave[14]
583                if (abs(gMSAWave[15])<1.0e-3)
584                     gMSAWave[15] = 0.0
585                endif
586                gMSAWave[10] = gMSAWave[16]
587        else
588                t2 = ta4*a2+ta5*b2+e12*(ta1*ck-ta2*sk)
589                t3 = ta4*a3+ta5*b3+e12*(ta1*sk-ta2*(ck-1.0))-4.0e-1*reta*(reta+1.0e1)-1.0
590
591//              MU(i)
592
593                um1 = t2*a2-e12*v2*v2
594                um2 = t1*a2+t2*a1-e24*v1*v2
595                um3 = t2*a3+t3*a2-e24*v2*v3
596                um4 = t1*a1-e12*v1*v1
597                um5 = t1*a3+t3*a1-e24*v1*v3
598                um6 = t3*a3-e12*v3*v3
599
600//                      GILLAN CONDITION ?
601//
602//                      YES - G(X=1+) = 0
603//
604//                      COEFFICIENTS AND FUNCTION VALUE
605//
606                IF ((IX==1) %| (IX==3))
607
608//                      NO - CALCULATE REMAINING COEFFICIENTS.
609
610//                      LAMBDA(I)
611
612                        al1 = e12*p2*p2
613                        al2 = e24*p1*p2-b2-b2
614                        al3 = e24*p2*p3
615                        al4 = e12*p1*p1-b1-b1
616                        al5 = e24*p1*p3-b3-b3-ak2
617                        al6 = e12*p3*p3
618
619//                      OMEGA(I)
620
621                        w16 = um1*al6-al1*um6
622                        w15 = um1*al5-al1*um5
623                        w14 = um1*al4-al1*um4
624                        w13 = um1*al3-al1*um3
625                        w12 = um1*al2-al1*um2
626
627                        w26 = um2*al6-al2*um6
628                        w25 = um2*al5-al2*um5
629                        w24 = um2*al4-al2*um4
630
631                        w36 = um3*al6-al3*um6
632                        w35 = um3*al5-al3*um5
633                        w34 = um3*al4-al3*um4
634                        w32 = um3*al2-al3*um2
635//
636                        w46 = um4*al6-al4*um6
637                        w56 = um5*al6-al5*um6
638                        w3526 = w35+w26
639                        w3425 = w34+w25
640
641//                      QUARTIC COEFFICIENTS
642
643                        w4 = w16*w16-w13*w36
644                        w3 = 2.0*w16*w15-w13*w3526-w12*w36
645                        w2 = w15*w15+2.0*w16*w14-w13*w3425-w12*w3526
646                        w1 = 2.0*w15*w14-w13*w24-w12*w3425
647                        w0 = w14*w14-w12*w24
648
649//                      ESTIMATE THE STARTING VALUE OF f
650
651                        if (ix==1)
652
653//                              LARGE K
654
655                                fap = (w14-w34-w46)/(w12-w15+w35-w26+w56-w32)
656                        else
657
658
659//                              ASSUME NOT TOO FAR FROM GILLAN CONDITION.
660//                              IF BOTH RGEK AND RAK ARE SMALL, USE P-W ESTIMATE.
661
662                                gMSAWave[14]=0.5*eta2d*dd2*exp(-rgek)
663                                if (( gMSAWave[11]<=2.0) %& ( gMSAWave[11]>=0.0) %& ( gMSAWave[12]<=1.0))
664                                        e24g = e24*rgek*exp(rak)
665                                        pwk = sqrt(e24g)
666                                        qpw = (1.0-sqrt(1.0+2.0*d2*d*pwk/eta22))*eta21/d
667                                        gMSAWave[14] = -qpw*qpw/e24+0.5*eta2d*dd2
668                                endif
669                                pg = p1+gMSAWave[14]
670                                ca = ak2*pg+2.0*(b3*pg-b1*p3)+e12*gMSAWave[14]*gMSAWave[14]*p3
671                                ca = -ca/(ak2*p2+2.0*(b3*p2-b2*p3))
672                                fap = -(pg+p2*ca)/p3
673                        endif
674
675//                      AND REFINE IT ACCORDING TO NEWTON
676
677                        ii=0
678                        do
679                                ii = ii+1
680                                if (ii>itm)
681
682//                                      FAILED TO CONVERGE IN ITM ITERATIONS
683
684                                        ir=-2
685                                        return ir
686                                endif
687                                fa = fap
688                                fun = w0+(w1+(w2+(w3+w4*fa)*fa)*fa)*fa
689                                fund = w1+(2.0*w2+(3.0*w3+4.0*w4*fa)*fa)*fa
690                                fap = fa-fun/fund
691                                del=abs((fap-fa)/fa)
692                        while (del>acc)
693                        ir = ir+ii
694                        fa = fap
695                        ca = -(w16*fa*fa+w15*fa+w14)/(w13*fa+w12)
696                        gMSAWave[14] = -(p1+p2*ca+p3*fa)
697                        gMSAWave[15] = gMSAWave[14]
698                        if (abs(gMSAWave[15])<1.0e-3)
699                              gMSAWave[15] = 0.0
700                        endif
701                        gMSAWave[10] = gMSAWave[16]
702                else
703                        ca = ak2*p1+2.0*(b3*p1-b1*p3)
704                        ca = -ca/(ak2*p2+2.0*(b3*p2-b2*p3))
705                        fa = -(p1+p2*ca)/p3
706                        if (ix==2)
707                                gMSAWave[15] = um1*ca*ca+(um2+um3*fa)*ca+um4+um5*fa+um6*fa*fa
708                        endif
709                        if (ix==4)
710                               gMSAWave[15] = -(p1+p2*ca+p3*fa)
711                        endif
712                endif
713                gMSAWave[3] = fa
714                gMSAWave[2] = ca
715                gMSAWave[1] = b1+b2*ca+b3*fa
716                gMSAWave[0] = a1+a2*ca+a3*fa
717                gMSAWave[8] = (v1+v2*ca+v3*fa)/gMSAWave[0]
718        endif
719        g24 = e24*rgek*ex1
720         gMSAWave[7] = (rak*ak2*ca-g24)/(ak2*g24)
721      return ir
722      end
723
724//////////////////////////////////////////////////////////
725//////////////////////////////////////////////////////////
726//
727//      CALCULATES S(Q) FOR "SQHPA"
728//
729//  **  THIS ROUTINE WORKS LOCALLY IN DOUBLE PRECISION **
730//
731//    JOHN B. HAYTER  (I.L.L.)    19-AUG-81
732//
733
734
735Function sqhcal(qq)
736        variable qq
737
738        variable SofQ,etaz,akz,gekz,e24,x1,x2,ck,sk,ak2,qk,q2k,qk2,qk3,qqk,sink,cosk,asink,qcosk,aqk,inter
739        WAVE gMSAWave = \$":HayPenMSA:gMSAWave"
740
741
742        //NVAR a=root:HayPenMSA:a
743        //NVAR b=root:HayPenMSA:b
744        //NVAR c=root:HayPenMSA:c
745        //NVAR f=root:HayPenMSA:f
746        //NVAR eta=root:HayPenMSA:eta
747        //NVAR gek=root:HayPenMSA:gek
748        //NVAR ak=root:HayPenMSA:ak
749        //NVAR u=root:HayPenMSA:u
750        //NVAR v=root:HayPenMSA:v
751        //NVAR gamk=root:HayPenMSA:gamk
752        //NVAR seta=root:HayPenMSA:seta
753        //NVAR sgek=root:HayPenMSA:sgek
754        //NVAR sak=root:HayPenMSA:sak
755        //NVAR scal=root:HayPenMSA:scal
756        //NVAR g1=root:HayPenMSA:g1
757
758      etaz = gMSAWave[10]
759      akz =  gMSAWave[12]
760      gekz =  gMSAWave[11]
761      e24 = 24.0*etaz
762      x1 = exp(akz)
763      x2 = 0.0
764      if ( gMSAWave[12]<20.0)
765             x2 = exp(-akz)
766      endif
767      ck = 0.5*(x1+x2)
768      sk = 0.5*(x1-x2)
769      ak2 = akz*akz
770
771      if (qq<=0.0)
772                SofQ = -1.0/gMSAWave[0]
773        else
774                qk = qq/gMSAWave[13]
775                q2k = qk*qk
776                qk2 = 1.0/q2k
777                qk3 = qk2/qk
778                qqk = 1.0/(qk*(q2k+ak2))
779                sink = sin(qk)
780                cosk = cos(qk)
781                asink = akz*sink
782                qcosk = qk*cosk
783                aqk = gMSAWave[0]*(sink-qcosk)
784        aqk=aqk+gMSAWave[1]*((2.0*qk2-1.0)*qcosk+2.0*sink-2.0/qk)
785                inter=24.0*qk3+4.0*(1.0-6.0*qk2)*sink
786                aqk=(aqk+0.5*etaz*gMSAWave[0]*(inter-(1.0-12.0*qk2+24.0*qk2*qk2)*qcosk))*qk3
787                aqk=aqk +gMSAWave[2]*(ck*asink-sk*qcosk)*qqk
788                aqk=aqk +gMSAWave[3]*(sk*asink-qk*(ck*cosk-1.0))*qqk
789                aqk=aqk +gMSAWave[3]*(cosk-1.0)*qk2
790                aqk=aqk -gekz*(asink+qcosk)*qqk
791                SofQ = 1.0/(1.0-e24*aqk)
792      endif
793      return SofQ
794end
795
796
797//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
798//
799// calculate U(r) given the HPMSA parameters
800//
801// pass in the x wave as scaled r/diam values
802//
803Function fHPMSA_Ur(w,x)
804        wave w
805        variable x
806
807//      variable timer=StartMSTimer
808
809        variable Elcharge=1.602189e-19          // electron charge in Coulombs (C)
810        variable kB=1.380662e-23                                // Boltzman constant in J/K
811        variable FrSpPerm=8.85418782E-12        //Permittivity of free space in C^2/(N m^2)
812
813        variable SofQ, QQ, Qdiam, gammaek, Vp, csalt, ss
814        variable VolFrac, SIdiam, diam, Kappa, cs, IonSt
815        variable dialec, Perm, beta1, Temp, zz, charge, ierr
816        Variable ur,gam,kapSig
817
818        diam=w[0]               //in A  (not SI .. should force people to think in nm!!!)
819        zz = w[1]               //# of charges
820        VolFrac=w[2]
821        QQ=x                    //in A^-1 (not SI .. should force people to think in nm^-1!!!)
822        Temp=w[3]               //in degrees Kelvin
823        csalt=w[4]              //in molarity
824        dialec=w[5]             // unitless
825
826
827////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
828//////////////////////////// convert to USEFUL inputs in SI units                                                //
829////////////////////////////    NOTE: easiest to do EVERYTHING in SI units                               //
830////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
831        beta1=1/(kB*Temp)               // in Joules^-1
832        Perm=dialec*FrSpPerm    //in C^2/(N  m^2)
833        charge=zz*Elcharge              //in Coulomb (C)
834        SIdiam = diam*1E-10             //in m
835        Vp=4*pi/3*(SIdiam/2)^3  //in m^3
836        cs=csalt*6.022E23*1E3   //# salt molecules/m^3
837
838//         Compute the derived values of :
839//                       Ionic strength IonSt (in C^2/m^3)
840//                      Kappa (Debye-Huckel screening length in m)
841//      and             gamma Exp(-k)
842        IonSt=0.5 * Elcharge^2*(zz*VolFrac/Vp+2*cs)
843        Kappa=sqrt(2*beta1*IonSt/Perm)     //Kappa calc from Ionic strength
844//      Kappa=2/SIdiam                                  // Use to compare with HP paper
845
846        kapSig = kappa*SIDiam
847
848        gam=beta1*charge^2/(pi*Perm*SIdiam*(2+Kappa*SIdiam)^2)*exp(kapSig)
849        ur = gam*exp(-kapSig*x)/x
850
851
852        Return(ur)
853End
Note: See TracBrowser for help on using the repository browser.