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