source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/PolyHardSphereInten_v40.ipf @ 510

Last change on this file since 510 was 510, checked in by srkline, 14 years ago

Simple change in all of the model function files to include the name of the parameter wave in the Keyword=list that is generated when a model is plotted. This is becoming an issue where the proper parameter wave can't be deduced from just the suffix, then there is nothing to put in the table.

I should have added this when I initially wrote the wrapper...

File size: 8.9 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 (A^-1) for model: "
19        Prompt qmax "Enter maximum q-value (A^-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,1e-6,6.3e-6,0.1}
24        make/o/t parameters_phs = {"Radius (A)","polydispersity","volume fraction","SLD sphere (A^-2)","SLD solvent (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 (A\\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","parameters_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,1e-6,6.3e-6,0.1}
53        make/o/t smear_parameters_phs = {"Radius (A)","polydispersity","volume fraction","SLD sphere (A^-2)","SLD solvent (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 (A\\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","smear_parameters_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,sld,slds
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        sld = w[3]
120        slds = w[4]
121        cont = (sld - slds)*1.0e4               // contrast (odd units)
122        bkg = w[5]              // background (1/cm)
123
124        zz=1/(z2*z2)-1.0
125        bb = sigma/(zz+1)
126        cc = zz+1
127
128//*c   Compute the number density by <r-cubed>, not <r> cubed*/
129        r1 = sigma/2.0
130        r3 = r1*r1*r1
131        r3 *= (zz+2)*(zz+3)/((zz+1)*(zz+1))
132        rho=phi/(1.3333333333*pi*r3)
133        t1 = rho*bb*cc
134        t2 = rho*bb*bb*cc*(cc+1)
135        t3 = rho*bb*bb*bb*cc*(cc+1)*(cc+2)
136        capd = 1-pi*t3/6
137//************
138        v1=1/(1+bb*bb*k*k)
139        v2=1/(4+bb*bb*k*k)
140        pp=(v1^(cc/2))*sin(cc*atan(bb*k))
141        p1=bb*cc*(v1^((cc+1)/2))*sin((cc+1)*atan(bb*k))
142        p2=cc*(cc+1)*bb*bb*(v1^((cc+2)/2))*sin((cc+2)*atan(bb*k))
143        mu=(2^cc)*(v2^(cc/2))*sin(cc*atan(bb*k/2))
144        mu1=(2^(cc+1))*bb*cc*(v2^((cc+1)/2))*sin((cc+1)*atan(k*bb/2))
145        s1=bb*cc
146        s2=cc*(cc+1)*bb*bb
147        s3=cc*(cc+1)*(cc+2)*bb*bb*bb
148        chi=(v1^(cc/2))*cos(cc*atan(bb*k))
149        chi1=bb*cc*(v1^((cc+1)/2))*cos((cc+1)*atan(bb*k))
150        chi2=cc*(cc+1)*bb*bb*(v1^((cc+2)/2))*cos((cc+2)*atan(bb*k))
151        ll=(2^cc)*(v2^(cc/2))*cos(cc*atan(bb*k/2))
152        l1=(2^(cc+1))*bb*cc*(v2^((cc+1)/2))*cos((cc+1)*atan(k*bb/2))
153        d1=(pi/capd)*(2+(pi/capd)*(t3-(rho/k)*(k*s3-p2)))
154        d2=((pi/capd)^2)*(rho/k)*(k*s2-p1)
155        d3=(-1.0)*((pi/capd)^2)*(rho/k)*(k*s1-pp)
156        d4=(pi/capd)*(k-(pi/capd)*(rho/k)*(chi1-s1))
157        d5=((pi/capd)^2)*((rho/k)*(chi-1)+0.5*k*t2)
158        d6=((pi/capd)^2)*(rho/k)*(chi2-s2)
159 // e1,e,y1,y evaluated in one big ugly line instead - no continuation character in IGOR
160//            e1=pow((pi/capd),2)*pow((rho/k/k),2)*((chi-1)*(chi2-s2)
161//       -(chi1-s1)*(chi1-s1)-(k*s1-p)*(k*s3-p2)+pow((k*s2-p1),2));
162//            e=1-(2*pi/capd)*(1+0.5*pi*t3/capd)*(rho/k/k/k)*(k*s1-p)
163//       -(2*pi/capd)*rho/k/k*((chi1-s1)+(0.25*pi*t2/capd)*(chi2-s2))-e1;
164//            y1=pow((pi/capd),2)*pow((rho/k/k),2)*((k*s1-p)*(chi2-s2)
165//       -2*(k*s2-p1)*(chi1-s1)+(k*s3-p2)*(chi-1));
166//            y = (2*pi/capd)*(1+0.5*pi*t3/capd)*(rho/k/k/k)
167//       *(chi+0.5*k*k*s2-1)-(2*pi*rho/capd/k/k)*(k*s2-p1
168//       +(0.25*pi*t2/capd)*(k*s3-p2))-y1;
169       
170        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))
171        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
172        y1=((pi/capd)^2)*((rho/k/k)^2)*((k*s1-pp)*(chi2-s2)-2*(k*s2-p1)*(chi1-s1)+(k*s3-p2)*(chi-1))
173        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       
174
175        capl=2.0*pi*cont*rho/k/k/k*(pp-0.5*k*(s1+chi1))
176        capl1=2.0*pi*cont*rho/k/k/k*(p1-0.5*k*(s2+chi2))
177        capmu=2.0*pi*cont*rho/k/k/k*(1-chi-0.5*k*p1)
178        capmu1=2.0*pi*cont*rho/k/k/k*(s1-chi1-0.5*k*p2)
179 
180        h1=capl*(capl*(yy*d1-ee*d6)+capl1*(yy*d2-ee*d4)+capmu*(ee*d1+yy*d6)+capmu1*(ee*d2+yy*d4))
181        h2=capl1*(capl*(yy*d2-ee*d4)+capl1*(yy*d3-ee*d5)+capmu*(ee*d2+yy*d4)+capmu1*(ee*d3+yy*d5))
182        h3=capmu*(capl*(ee*d1+yy*d6)+capl1*(ee*d2+yy*d4)+capmu*(ee*d6-yy*d1)+capmu1*(ee*d4-yy*d2))
183        h4=capmu1*(capl*(ee*d2+yy*d4)+capl1*(ee*d3+yy*d5)+capmu*(ee*d4-yy*d2)+capmu1*(ee*d5-yy*d3))
184
185//*  This part computes the second integral in equation (1) of the paper.*/
186
187        hint1 = -2.0*(h1+h2+h3+h4)/(k*k*k*(ee*ee+yy*yy))
188 
189//*  This part computes the first integral in equation (1).  It also
190// generates the KC approximated effective structure factor.*/
191
192        pq=4*pi*cont*(sin(k*sigma/2)-0.5*k*sigma*cos(k*sigma/2))
193        hint2=8*pi*pi*rho*cont*cont/(k*k*k*k*k*k)*(1-chi-k*p1+0.25*k*k*(s2+chi2))
194 
195        ka=k*(sigma/2)
196//
197        hh=hint1+hint2          // this is the model intensity
198//
199        heff=1.0+hint1/hint2
200        ka2=ka*ka
201//*
202//  heff is PY analytical solution for intensity divided by the
203//   form factor.  happ is the KC approximated effective S(q)
204 
205//*******************
206//  add in the background then return the intensity value
207
208        return (hh+bkg)
209
210End   // end of fcngrif()
211
212// this is all there is to the smeared calculation!
213Function SmearedPolyHardSpheres(s) :FitFunc
214        Struct ResSmearAAOStruct &s
215
216////the name of your unsmeared model is the first argument
217        Smear_Model_20(PolyHardSpheres,s.coefW,s.xW,s.yW,s.resW)
218
219        return(0)
220End
221
222//wrapper to calculate the smeared model as an AAO-Struct
223// fills the struct and calls the ususal function with the STRUCT parameter
224//
225// used only for the dependency, not for fitting
226//
227Function fSmearedPolyHardSpheres(coefW,yW,xW)
228        Wave coefW,yW,xW
229       
230        String str = getWavesDataFolder(yW,0)
231        String DF="root:"+str+":"
232       
233        WAVE resW = $(DF+str+"_res")
234       
235        STRUCT ResSmearAAOStruct fs
236        WAVE fs.coefW = coefW   
237        WAVE fs.yW = yW
238        WAVE fs.xW = xW
239        WAVE fs.resW = resW
240       
241        Variable err
242        err = SmearedPolyHardSpheres(fs)
243       
244        return (0)
245End
Note: See TracBrowser for help on using the repository browser.