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

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

Modified all procedure files to add to the keyword=value strings to identify the function, coefficients, and suffix once the model is plotted (and the objects will exist)

a one-liner in the Plot and PlotSmeared? macros.

  • necessary for smoother functioning of the wrapper panel.
File size: 8.7 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
3
4////////////////////////////////////////////////
5// GaussUtils.proc and PlotUtils.proc MUST be included for the smearing calculation to compile
6// Adopting these into the experiment will insure that they are always present
7////////////////////////////////////////////////
8// this example is for the scattered intensity from a dense dispersion of polydisperse spheres
9// hard sphere interactions are included (exact, multicomponent Percus-Yevick)
10// the polydispersity in radius is a Schulz distribution
11//
12// 06 NOV 98 SRK
13////////////////////////////////////////////////
14
15Proc PlotPolyHardSpheres(num,qmin,qmax)
16        Variable num=128,qmin=0.001,qmax=0.7
17        Prompt num "Enter number of data points for model: "
18        Prompt qmin "Enter minimum q-value (^-1) for model: "
19        Prompt qmax "Enter maximum q-value (^-1) for model: "
20       
21        Make/O/D/n=(num) xwave_phs,ywave_phs
22        xwave_phs =alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
23        Make/O/D coef_phs = {100,0.12,0.1,2.0e-6,0.1}
24        make/o/t parameters_phs = {"Radius (A)","polydispersity","volume fraction","contrast (A^-2)","background (cm^-1)"}
25        Edit parameters_phs,coef_phs
26        Variable/G root:g_phs
27        g_phs := PolyHardSpheres(coef_phs,ywave_phs,xwave_phs)
28//      ywave_phs := PolyHardSpheres(coef_phs,xwave_phs)
29        Display ywave_phs vs xwave_phs
30        ModifyGraph log=1,marker=29,msize=2,mode=4
31        Label bottom "q (\\S-1\\M)"
32        Label left "Intensity (cm\\S-1\\M)"
33        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
34       
35        AddModelToStrings("PolyHardSpheres","coef_phs","phs")
36End
37
38///////////////////////////////////////////////////////////
39// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
40Proc PlotSmearedPolyHardSpheres(str)                                                           
41        String str
42        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
43       
44        // if any of the resolution waves are missing => abort
45        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
46                Abort
47        endif
48       
49        SetDataFolder $("root:"+str)
50       
51        // Setup parameter table for model function
52        Make/O/D smear_coef_phs = {100,0.12,0.1,2.0e-6,0.1}
53        make/o/t smear_parameters_phs = {"Radius (A)","polydispersity","volume fraction","contrast (A^-2)","background (cm^-1)"}
54        Edit smear_parameters_phs,smear_coef_phs
55       
56        // output smeared intensity wave, dimensions are identical to experimental QSIG values
57        // make extra copy of experimental q-values for easy plotting
58        Duplicate/O $(str+"_q") smeared_phs,smeared_qvals
59        SetScale d,0,0,"1/cm",smeared_phs                                       
60               
61        Variable/G gs_phs=0
62        gs_phs := fSmearedPolyHardSpheres(smear_coef_phs,smeared_phs,smeared_qvals)     //this wrapper fills the STRUCT
63
64        Display smeared_phs vs smeared_qvals
65        ModifyGraph log=1,marker=29,msize=2,mode=4
66        Label bottom "q (\\S-1\\M)"
67        Label left "Intensity (cm\\S-1\\M)"
68        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2) 
69       
70        SetDataFolder root:
71        AddModelToStrings("SmearedPolyHardSpheres","smear_coef_phs","phs")
72End
73
74//AAO version
75Function PolyHardSpheres(cw,yw,xw) : FitFunc
76        Wave cw,yw,xw
77
78#if exists("PolyHardSpheresX")
79        yw = PolyHardSpheresX(cw,xw)
80#else
81        yw = fPolyHardSpheres(cw,xw)
82#endif
83        return(0)
84End
85
86///////////////////////////////////////////////////////////////
87// unsmeared model calculation
88//   This program calculates the effective structure factor for a suspension
89//   of spheres whose size distribution is given by a Schulz distribution
90//   PY closure was used to solve.  Equations are analytical.
91//   Follows paper by W.L. Griffith, Phys. Rev. A 35 (5) p.2200 1987
92//  Original coding (F) by Jon Bender, U. Delaware
93//      converted to c 2/97 SRK
94// converted to IGOR 12/97 SRK
95//
96// replace single letter variables like "e" with "ee" (to be done MAY04)
97//
98///////////////////////////
99Function fPolyHardSpheres(w,k) : FitFunc
100        Wave w                  // the coefficient wave
101        Variable k              // the x values, as a variable (single k is OK)
102
103        // assign local variables
104        Variable mu,mu1,d1,d2,d3,d4,d5,d6,capd,rho,beta
105        Variable ll,l1,bb,cc,chi,chi1,chi2,ee,t1,t2,t3,pp
106        Variable ka,zz,v1,v2,p1,p2
107        Variable h1,h2,h3,h4,e1,yy,y1,ss,s1,s2,s3,hint1,hint2
108        Variable capl,capl1,capmu,capmu1,r3,pq,happ
109        Variable ka2,r1,heff
110        Variable hh
111         
112        //* reassign names to the variable set */
113        Variable rad,z2,phi,cont,bkg,sigma
114         
115        rad = w[0]              // radius (A)
116        sigma = 2*rad
117        z2 = w[1]               //polydispersity (0<z2<1)
118        phi = w[2]              // volume fraction (0<phi<1)
119        cont = w[3]*1.0e4               // contrast (odd units)
120        bkg = w[4]              // background (1/cm)
121
122        zz=1/(z2*z2)-1.0
123        bb = sigma/(zz+1)
124        cc = zz+1
125
126//*c   Compute the number density by <r-cubed>, not <r> cubed*/
127        r1 = sigma/2.0
128        r3 = r1*r1*r1
129        r3 *= (zz+2)*(zz+3)/((zz+1)*(zz+1))
130        rho=phi/(1.3333333333*pi*r3)
131        t1 = rho*bb*cc
132        t2 = rho*bb*bb*cc*(cc+1)
133        t3 = rho*bb*bb*bb*cc*(cc+1)*(cc+2)
134        capd = 1-pi*t3/6
135//************
136        v1=1/(1+bb*bb*k*k)
137        v2=1/(4+bb*bb*k*k)
138        pp=(v1^(cc/2))*sin(cc*atan(bb*k))
139        p1=bb*cc*(v1^((cc+1)/2))*sin((cc+1)*atan(bb*k))
140        p2=cc*(cc+1)*bb*bb*(v1^((cc+2)/2))*sin((cc+2)*atan(bb*k))
141        mu=(2^cc)*(v2^(cc/2))*sin(cc*atan(bb*k/2))
142        mu1=(2^(cc+1))*bb*cc*(v2^((cc+1)/2))*sin((cc+1)*atan(k*bb/2))
143        s1=bb*cc
144        s2=cc*(cc+1)*bb*bb
145        s3=cc*(cc+1)*(cc+2)*bb*bb*bb
146        chi=(v1^(cc/2))*cos(cc*atan(bb*k))
147        chi1=bb*cc*(v1^((cc+1)/2))*cos((cc+1)*atan(bb*k))
148        chi2=cc*(cc+1)*bb*bb*(v1^((cc+2)/2))*cos((cc+2)*atan(bb*k))
149        ll=(2^cc)*(v2^(cc/2))*cos(cc*atan(bb*k/2))
150        l1=(2^(cc+1))*bb*cc*(v2^((cc+1)/2))*cos((cc+1)*atan(k*bb/2))
151        d1=(pi/capd)*(2+(pi/capd)*(t3-(rho/k)*(k*s3-p2)))
152        d2=((pi/capd)^2)*(rho/k)*(k*s2-p1)
153        d3=(-1.0)*((pi/capd)^2)*(rho/k)*(k*s1-pp)
154        d4=(pi/capd)*(k-(pi/capd)*(rho/k)*(chi1-s1))
155        d5=((pi/capd)^2)*((rho/k)*(chi-1)+0.5*k*t2)
156        d6=((pi/capd)^2)*(rho/k)*(chi2-s2)
157 // e1,e,y1,y evaluated in one big ugly line instead - no continuation character in IGOR
158//            e1=pow((pi/capd),2)*pow((rho/k/k),2)*((chi-1)*(chi2-s2)
159//       -(chi1-s1)*(chi1-s1)-(k*s1-p)*(k*s3-p2)+pow((k*s2-p1),2));
160//            e=1-(2*pi/capd)*(1+0.5*pi*t3/capd)*(rho/k/k/k)*(k*s1-p)
161//       -(2*pi/capd)*rho/k/k*((chi1-s1)+(0.25*pi*t2/capd)*(chi2-s2))-e1;
162//            y1=pow((pi/capd),2)*pow((rho/k/k),2)*((k*s1-p)*(chi2-s2)
163//       -2*(k*s2-p1)*(chi1-s1)+(k*s3-p2)*(chi-1));
164//            y = (2*pi/capd)*(1+0.5*pi*t3/capd)*(rho/k/k/k)
165//       *(chi+0.5*k*k*s2-1)-(2*pi*rho/capd/k/k)*(k*s2-p1
166//       +(0.25*pi*t2/capd)*(k*s3-p2))-y1;
167       
168        e1=((pi/capd)^2)*((rho/k/k)^2)*((chi-1)*(chi2-s2)-(chi1-s1)*(chi1-s1)-(k*s1-pp)*(k*s3-p2)+((k*s2-p1)^2))
169        ee=1-(2*pi/capd)*(1+0.5*pi*t3/capd)*(rho/k/k/k)*(k*s1-pp)-(2*pi/capd)*rho/k/k*((chi1-s1)+(0.25*pi*t2/capd)*(chi2-s2))-e1
170        y1=((pi/capd)^2)*((rho/k/k)^2)*((k*s1-pp)*(chi2-s2)-2*(k*s2-p1)*(chi1-s1)+(k*s3-p2)*(chi-1))
171        yy = (2*pi/capd)*(1+0.5*pi*t3/capd)*(rho/k/k/k)*(chi+0.5*k*k*s2-1)-(2*pi*rho/capd/k/k)*(k*s2-p1+(0.25*pi*t2/capd)*(k*s3-p2))-y1       
172
173        capl=2.0*pi*cont*rho/k/k/k*(pp-0.5*k*(s1+chi1))
174        capl1=2.0*pi*cont*rho/k/k/k*(p1-0.5*k*(s2+chi2))
175        capmu=2.0*pi*cont*rho/k/k/k*(1-chi-0.5*k*p1)
176        capmu1=2.0*pi*cont*rho/k/k/k*(s1-chi1-0.5*k*p2)
177 
178        h1=capl*(capl*(yy*d1-ee*d6)+capl1*(yy*d2-ee*d4)+capmu*(ee*d1+yy*d6)+capmu1*(ee*d2+yy*d4))
179        h2=capl1*(capl*(yy*d2-ee*d4)+capl1*(yy*d3-ee*d5)+capmu*(ee*d2+yy*d4)+capmu1*(ee*d3+yy*d5))
180        h3=capmu*(capl*(ee*d1+yy*d6)+capl1*(ee*d2+yy*d4)+capmu*(ee*d6-yy*d1)+capmu1*(ee*d4-yy*d2))
181        h4=capmu1*(capl*(ee*d2+yy*d4)+capl1*(ee*d3+yy*d5)+capmu*(ee*d4-yy*d2)+capmu1*(ee*d5-yy*d3))
182
183//*  This part computes the second integral in equation (1) of the paper.*/
184
185        hint1 = -2.0*(h1+h2+h3+h4)/(k*k*k*(ee*ee+yy*yy))
186 
187//*  This part computes the first integral in equation (1).  It also
188// generates the KC approximated effective structure factor.*/
189
190        pq=4*pi*cont*(sin(k*sigma/2)-0.5*k*sigma*cos(k*sigma/2))
191        hint2=8*pi*pi*rho*cont*cont/(k*k*k*k*k*k)*(1-chi-k*p1+0.25*k*k*(s2+chi2))
192 
193        ka=k*(sigma/2)
194//
195        hh=hint1+hint2          // this is the model intensity
196//
197        heff=1.0+hint1/hint2
198        ka2=ka*ka
199//*
200//  heff is PY analytical solution for intensity divided by the
201//   form factor.  happ is the KC approximated effective S(q)
202 
203//*******************
204//  add in the background then return the intensity value
205
206        return (hh+bkg)
207
208End   // end of fcngrif()
209
210// this is all there is to the smeared calculation!
211Function SmearedPolyHardSpheres(s) :FitFunc
212        Struct ResSmearAAOStruct &s
213
214////the name of your unsmeared model is the first argument
215        Smear_Model_20(PolyHardSpheres,s.coefW,s.xW,s.yW,s.resW)
216
217        return(0)
218End
219
220//wrapper to calculate the smeared model as an AAO-Struct
221// fills the struct and calls the ususal function with the STRUCT parameter
222//
223// used only for the dependency, not for fitting
224//
225Function fSmearedPolyHardSpheres(coefW,yW,xW)
226        Wave coefW,yW,xW
227       
228        String str = getWavesDataFolder(yW,0)
229        String DF="root:"+str+":"
230       
231        WAVE resW = $(DF+str+"_res")
232       
233        STRUCT ResSmearAAOStruct fs
234        WAVE fs.coefW = coefW   
235        WAVE fs.yW = yW
236        WAVE fs.xW = xW
237        WAVE fs.resW = resW
238       
239        Variable err
240        err = SmearedPolyHardSpheres(fs)
241       
242        return (0)
243End
Note: See TracBrowser for help on using the repository browser.