source: sans/Release/trunk/NCNR_User_Procedures/SANS/Analysis/Models/HPMSA_v40.ipf @ 381

Last change on this file since 381 was 381, checked in by srkline, 14 years ago

Merging Dev/trunk revision 374+ into Release/trunk for version 6.004

File size: 22.3 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
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("root:HayPenMSA"))
14                NewDataFolder root:HayPenMSA
15        endif
16        Make/O/D/N=17 root: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       
34        AddModelToStrings("HayterPenfoldMSA","coef_hpmsa","hpmsa")
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]
128//use wave instead
129        WAVE gMSAWave = $"root: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 = $"root: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 = $"root: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 = $"root: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     
Note: See TracBrowser for help on using the repository browser.