source: sans/Analysis/trunk/Put in User Procedures/SANS_Models_v3.00/HPMSA.ipf @ 42

Last change on this file since 42 was 42, checked in by srkline, 16 years ago

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