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

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

LOTS of changes to the analysis ipf files:

-- see sphere.ipf for the simplest example of the changes --

  • #pragma Igor 6
  • #if directive to look for XOP
  • AAO unsmeared functions
  • STRUCT functions for smearing (also AAO)
  • new wrappers for dependencies to struct functions

(2006 models have NOT been completed yet, only the old models)

  • loading data files into data folders (PlotUtils?) + some streamlining of the loaders
  • Smear_Model_N is now AAO + some streamlining of the quadrature code

-- SHS and SW structure factor XOPs are crashing (need DP wave, I may have old XOP)
-- this breaks fitting of the smeared models until wrappers can be devised
-- all packages will be broken due to the new data folder structure
-- multiple instances of functions will now cause problems (MSA)
-- RPA model is a problem with its odd functional form (extra wave)

-- lots of other carnage to follow as the bugs and typos are shaken out

24 JUL 2007 SRK

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