source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/NewModels_2006/LogNormalSphere.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.4 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
3
4#include "sphere"
5// plots the form factor of spherical particles with a log-normal radius distribution
6//
7// for the integration it may be better to use adaptive routine
8
9//Proc to setup data and coefficients to plot the model
10Proc PlotLogNormalSphere(num,qmin,qmax)
11        Variable num=128,qmin=0.001,qmax=0.7
12        Prompt num "Enter number of data points for model: "
13        Prompt qmin "Enter minimum q-value (^-1) for model: "
14        Prompt qmax "Enter maximum q-value (^-1) for model: "
15       
16        Make/O/D/N=(num) xwave_lns,ywave_lns
17        xwave_lns = alog( log(qmin) + x*((log(qmax)-log(qmin))/num) )
18        Make/O/D coef_lns = {0.01,60,0.2,1e-6,2e-6,0}
19        make/O/T parameters_lns = {"Volume Fraction (scale)","exp(mu)=median Radius (A)","sigma","SLD sphere (A-2)","SLD solvent (A-2)","bkg (cm-1 sr-1)"}
20        Edit parameters_lns,coef_lns
21       
22        Variable/G root:g_lns
23        g_lns := LogNormalSphere(coef_lns,ywave_lns,xwave_lns)
24        Display ywave_lns vs xwave_lns
25        ModifyGraph log=1,marker=29,msize=2,mode=4
26        Label bottom "q (\\S-1\\M)"
27        Label left "Intensity (cm\\S-1\\M)"
28        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
29       
30        AddModelToStrings("LogNormalSphere","coef_lns","lns")
31End
32
33// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
34Proc PlotSmearedLogNormalSphere(str)                                                           
35        String str
36        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
37       
38        // if any of the resolution waves are missing => abort
39        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
40                Abort
41        endif
42       
43        SetDataFolder $("root:"+str)
44       
45        // Setup parameter table for model function
46        Make/O/D smear_coef_lns = {0.01,60,0.2,1e-6,2e-6,0}                                     
47        make/o/t smear_parameters_lns = {"Volume Fraction (scale)","exp(mu)=median Radius (A)","sigma","SLD sphere (A-2)","SLD solvent (A-2)","bkg (cm-1 sr-1)"}               
48        Edit smear_parameters_lns,smear_coef_lns                                       
49       
50        // output smeared intensity wave, dimensions are identical to experimental QSIG values
51        // make extra copy of experimental q-values for easy plotting
52        Duplicate/O $(str+"_q") smeared_lns,smeared_qvals                               
53        SetScale d,0,0,"1/cm",smeared_lns                                                       
54                                       
55        Variable/G gs_lns=0
56        gs_lns := fSmearedLogNormalSphere(smear_coef_lns,smeared_lns,smeared_qvals)     //this wrapper fills the STRUCT
57       
58        Display smeared_lns vs smeared_qvals                                                                   
59        ModifyGraph log=1,marker=29,msize=2,mode=4
60        Label bottom "q (\\S-1\\M)"
61        Label left "Intensity (cm\\S-1\\M)"
62        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
63       
64        SetDataFolder root:
65        AddModelToStrings("SmearedLogNormalSphere","smear_coef_lns","lns")
66End
67
68
69
70
71//AAO version, uses XOP if available
72// simply calls the original single point calculation with
73// a wave assignment (this will behave nicely if given point ranges)
74Function LogNormalSphere(cw,yw,xw) : FitFunc
75        Wave cw,yw,xw
76       
77#if exists("LogNormalSphereX")
78        yw = LogNormalSphereX(cw,xw)
79#else
80        yw = fLogNormalSphere(cw,xw)
81#endif
82        return(0)
83End
84
85
86// calculates the model at each q-value by integrating over the normalized size distribution
87// integration is done by gauss quadrature of either 20 or 76 points (nord)
88// 76 points is slower, but reccommended to remove high-q oscillations
89//
90Function fLogNormalSphere(w,xx): FitFunc
91        wave w
92        variable xx
93       
94        Variable scale,rad,sig,rho,rhos,bkg,delrho,mu,r3
95       
96        //set up the coefficient values
97        scale=w[0]
98        rad=w[1]                //rad is the median radius
99        mu = ln(w[1])
100        sig=w[2]
101        rho=w[3]
102        rhos=w[4]
103        delrho=rho-rhos
104        bkg=w[5]
105       
106//temp set scale=1 and bkg=0 for quadrature calc
107        Make/O/D/N=4 sphere_temp
108        sphere_temp[0] = 1
109        sphere_temp[1] = rad            //changed in loop
110        sphere_temp[2] = delrho
111        sphere_temp[3] = 0
112       
113        //currently using 20 pts...
114        Variable va,vb,ii,zi,nord,yy,summ,inten
115        String weightStr,zStr
116       
117        //select number of gauss points by setting nord=20 or76 points
118//      nord = 20
119        nord = 76
120       
121        weightStr = "gauss"+num2str(nord)+"wt"
122        zStr = "gauss"+num2str(nord)+"z"
123       
124        if (WaveExists($weightStr) == 0) // wave reference is not valid,
125                Make/D/N=(nord) $weightStr,$zStr
126                Wave gauWt = $weightStr
127                Wave gauZ = $zStr               // wave references to pass
128                if(nord==20)
129                        Make20GaussPoints(gauWt,gauZ)
130                else
131                        Make76GaussPoints(gauWt,gauZ)
132                endif   
133        else
134                if(exists(weightStr) > 1)
135                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
136                endif
137                Wave gauWt = $weightStr
138                Wave gauZ = $zStr               // create the wave references
139        endif
140       
141        // end points of integration
142        // limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero
143        // +/- 3 sigq catches 99.73% of distrubution
144        // change limits (and spacing of zi) at each evaluation based on R()
145        //integration from va to vb
146       
147//      va = -3*sig + rad
148        va = -3.5*sig +mu               //in ln(R) space
149        va = exp(va)                    //in R space
150        if (va<0)
151                va=0            //to avoid numerical error when  va<0 (-ve q-value)
152        endif
153//      vb = 3*exp(sig) +rad
154        vb = 3.5*sig*(1+sig)+ mu
155        vb = exp(vb)
156       
157        summ = 0.0              // initialize integral
158        Make/O/N=1 tmp_yw,tmp_xw
159        tmp_xw[0] = xx
160        for(ii=0;ii<nord;ii+=1)
161                // calculate Gauss points on integration interval (r-value for evaluation)
162                zi = ( gauZ[ii]*(vb-va) + vb + va )/2.0
163                sphere_temp[1] = zi
164                // calculate sphere scattering
165                SphereForm(sphere_temp,tmp_yw,tmp_xw)                   //AAO calculation, one point wave
166                yy = gauWt[ii] *  LogNormal_distr(sig,mu,zi) * tmp_yw[0]
167                yy *= 4*pi/3*zi*zi*zi           //un-normalize by current sphere volume
168               
169                summ += yy              //add to the running total of the quadrature
170        endfor
171// calculate value of integral to return
172        inten = (vb-va)/2.0*summ
173       
174        //re-normalize by polydisperse sphere volume
175        //third moment
176        r3 = exp(3*mu + 9/2*sig^2)              // <R^3> directly
177        inten /= (4*pi/3*r3)            //polydisperse volume
178       
179        inten *= scale
180        inten+=bkg
181       
182        Return(inten)
183End
184
185//wrapper to calculate the smeared model as an AAO-Struct
186// fills the struct and calls the ususal function with the STRUCT parameter
187//
188// used only for the dependency, not for fitting
189//
190Function fSmearedLogNormalSphere(coefW,yW,xW)
191        Wave coefW,yW,xW
192       
193        String str = getWavesDataFolder(yW,0)
194        String DF="root:"+str+":"
195       
196        WAVE resW = $(DF+str+"_res")
197       
198        STRUCT ResSmearAAOStruct fs
199        WAVE fs.coefW = coefW   
200        WAVE fs.yW = yW
201        WAVE fs.xW = xW
202        WAVE fs.resW = resW
203       
204        Variable err
205        err = SmearedLogNormalSphere(fs)
206       
207        return (0)
208End
209
210// this is all there is to the smeared calculation!
211Function SmearedLogNormalSphere(s) :FitFunc
212        Struct ResSmearAAOStruct &s
213
214//      the name of your unsmeared model (AAO) is the first argument
215        Smear_Model_20(LogNormalSphere,s.coefW,s.xW,s.yW,s.resW)
216
217        return(0)
218End
219       
220
221
222// normalization is correct, using 3rd moment of lognormal distribution
223//
224Function LogNormal_distr(sig,mu,pt)
225        Variable sig,mu,pt
226       
227        Variable retval
228       
229        retval = (1/ ( sig*pt*sqrt(2*pi)) )*exp( -0.5*((ln(pt) - mu)^2)/sig^2 )
230        return(retval)
231End
232
233//calculates number density given the coefficients of the lognormal distribution
234// the scale factor is the volume fraction
235// then nden = phi/<V> where <V> is calculated using the 3rd moment of the radius
236Macro NumberDensity_LogN()
237       
238        Variable nden,r3,rg,sv,i0,ravg,rpk
239       
240        if(WaveExists(coef_lns)==0)
241                abort "You need to plot the model first to create the coefficient table"
242        Endif
243       
244        Print "median radius (A) = ",coef_lns[1]
245        Print "sigma = ",coef_lns[2]
246        Print "volume fraction = ",coef_lns[0]
247       
248        r3 = exp(3*ln(coef_lns[1]) + 9/2*coef_lns[2]^2)         // <R^3> directly,[A^3]
249        nden = coef_lns[0]/(4*pi/3*r3)          //nden in 1/A^3
250        ravg = exp(ln(coef_lns[1]) + 0.5*coef_lns[2]^2)
251        rpk = exp(ln(coef_lns[1]) - coef_lns[2]^2)
252        rg = (3./5.)^0.5*exp(ln(coef_lns[1]) + 7.*coef_lns[2]^2)
253        sv = 1.0e8*3*coef_lns[0]*exp(-ln(coef_lns[1]) - 2.5*coef_lns[2]^2)
254        i0 = 1.0e8*(4*pi/3)*coef_lns[0]*(coef_lns[3]-coef_lns[4])^2*exp(3*ln(coef_lns[1]) + 13.5*coef_lns[2]^2)
255       
256        Print "number density (A^-3) = ",nden
257        Print "mean radius (A) = ",ravg
258        Print "peak dis. radius (A) = ",rpk
259        Print "Guinier radius (A) = ",rg
260        Print "Interfacial surface area / volume (cm-1) Sv = ",sv
261        Print "Forward cross section (cm-1 sr-1) I(0) = ",i0
262End
263
264// plots the lognormal distribution based on the coefficient values
265// a static calculation, so re-run each time
266//
267Macro PlotLogNormalDistribution()
268
269        variable sig,mu,maxr
270       
271        if(WaveExists(coef_lns)==0)
272                abort "You need to plot the model first to create the coefficient table"
273        Endif
274        sig=coef_lns[2]
275        mu = ln(coef_lns[1])
276       
277        Make/O/D/N=1000 lognormal_distribution
278        maxr = 5*sig*(1+sig)+ mu
279        maxr = exp(maxr)
280        SetScale/I x, 0, maxr, lognormal_distribution
281        lognormal_distribution = LogNormal_distr(sig,mu,x)
282        Display lognormal_distribution
283        modifygraph log(bottom)=1
284        legend
285End
Note: See TracBrowser for help on using the repository browser.