source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/HPMSA.ipf @ 146

Last change on this file since 146 was 130, checked in by srkline, 15 years ago

-Fixed RPA and HPMSA functions to be up-to-date and use data folders properly for the smeared calculations.

-Changed the coef/param naming scheme for RPA to be consistent with all other models.

-SANSModelPicker now finds Smeared Plot Procs with the correct number of parameters to populate the menus

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