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

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

Changed Plot* and PlotSmeared?* naming schemes to be all consistent prefixes for the actual function name, so that the macros can be constructed from the function name, or vice versa.

also some tweaks to the wrapper to make sure that plot and append really work

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