source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/LogNormalSphere_v40.ipf @ 570

Last change on this file since 570 was 570, checked in by srkline, 13 years ago

Change (1):
In preparation for release, updated pragma IgorVersion?=6.1 in all procedures

Change (2):
As a side benefit of requiring 6.1, we can use the MultiThread? keyword to thread any model function we like. The speed benefit is only noticeable on functions that require at least one integration and at least 100 points (resolution smearing is NOT threaded, too many threadSafe issues, too little benefit). I have chosen to use the MultiThread? only on the XOP assignment. In the Igor code there are too many functions that are not explicitly declared threadsafe, making for a mess.

File size: 8.6 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.1
3
4#include "sphere_v40"
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 (A^-1) for model: "
14        Prompt qmax "Enter maximum q-value (A^-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 (A\\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","parameters_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 (A\\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","smear_parameters_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//      Variable t1=StopMSTimer(-2)
78
79#if exists("LogNormalSphereX")
80        MultiThread yw = LogNormalSphereX(cw,xw)
81#else
82        yw = fLogNormalSphere(cw,xw)
83#endif
84
85//      Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
86
87        return(0)
88End
89
90
91// calculates the model at each q-value by integrating over the normalized size distribution
92// integration is done by gauss quadrature of either 20 or 76 points (nord)
93// 76 points is slower, but reccommended to remove high-q oscillations
94//
95Function fLogNormalSphere(w,xx): FitFunc
96        wave w
97        variable xx
98       
99        Variable scale,rad,sig,rho,rhos,bkg,delrho,mu,r3
100       
101        //set up the coefficient values
102        scale=w[0]
103        rad=w[1]                //rad is the median radius
104        mu = ln(w[1])
105        sig=w[2]
106        rho=w[3]
107        rhos=w[4]
108        delrho=rho-rhos
109        bkg=w[5]
110       
111//temp set scale=1 and bkg=0 for quadrature calc
112        Make/O/D/N=4 sphere_temp
113        sphere_temp[0] = 1
114        sphere_temp[1] = rad            //changed in loop
115        sphere_temp[2] = rho
116        sphere_temp[3] = rhos
117        sphere_temp[4] = 0
118       
119        //currently using 20 pts...
120        Variable va,vb,ii,zi,nord,yy,summ,inten
121        String weightStr,zStr
122       
123        //select number of gauss points by setting nord=20 or76 points
124//      nord = 20
125        nord = 76
126       
127        weightStr = "gauss"+num2str(nord)+"wt"
128        zStr = "gauss"+num2str(nord)+"z"
129       
130        if (WaveExists($weightStr) == 0) // wave reference is not valid,
131                Make/D/N=(nord) $weightStr,$zStr
132                Wave gauWt = $weightStr
133                Wave gauZ = $zStr               // wave references to pass
134                if(nord==20)
135                        Make20GaussPoints(gauWt,gauZ)
136                else
137                        Make76GaussPoints(gauWt,gauZ)
138                endif   
139        else
140                if(exists(weightStr) > 1)
141                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
142                endif
143                Wave gauWt = $weightStr
144                Wave gauZ = $zStr               // create the wave references
145        endif
146       
147        // end points of integration
148        // limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero
149        // +/- 3 sigq catches 99.73% of distrubution
150        // change limits (and spacing of zi) at each evaluation based on R()
151        //integration from va to vb
152       
153//      va = -3*sig + rad
154        va = -3.5*sig +mu               //in ln(R) space
155        va = exp(va)                    //in R space
156        if (va<0)
157                va=0            //to avoid numerical error when  va<0 (-ve q-value)
158        endif
159//      vb = 3*exp(sig) +rad
160        vb = 3.5*sig*(1+sig)+ mu
161        vb = exp(vb)
162       
163        summ = 0.0              // initialize integral
164        Make/O/N=1 tmp_yw,tmp_xw
165        tmp_xw[0] = xx
166        for(ii=0;ii<nord;ii+=1)
167                // calculate Gauss points on integration interval (r-value for evaluation)
168                zi = ( gauZ[ii]*(vb-va) + vb + va )/2.0
169                sphere_temp[1] = zi
170                // calculate sphere scattering
171                SphereForm(sphere_temp,tmp_yw,tmp_xw)                   //AAO calculation, one point wave
172                yy = gauWt[ii] *  LogNormal_distr(sig,mu,zi) * tmp_yw[0]
173                yy *= 4*pi/3*zi*zi*zi           //un-normalize by current sphere volume
174               
175                summ += yy              //add to the running total of the quadrature
176        endfor
177// calculate value of integral to return
178        inten = (vb-va)/2.0*summ
179       
180        //re-normalize by polydisperse sphere volume
181        //third moment
182        r3 = exp(3*mu + 9/2*sig^2)              // <R^3> directly
183        inten /= (4*pi/3*r3)            //polydisperse volume
184       
185        inten *= scale
186        inten+=bkg
187       
188        Return(inten)
189End
190
191//wrapper to calculate the smeared model as an AAO-Struct
192// fills the struct and calls the ususal function with the STRUCT parameter
193//
194// used only for the dependency, not for fitting
195//
196Function fSmearedLogNormalSphere(coefW,yW,xW)
197        Wave coefW,yW,xW
198       
199        String str = getWavesDataFolder(yW,0)
200        String DF="root:"+str+":"
201       
202        WAVE resW = $(DF+str+"_res")
203       
204        STRUCT ResSmearAAOStruct fs
205        WAVE fs.coefW = coefW   
206        WAVE fs.yW = yW
207        WAVE fs.xW = xW
208        WAVE fs.resW = resW
209       
210        Variable err
211        err = SmearedLogNormalSphere(fs)
212       
213        return (0)
214End
215
216// this is all there is to the smeared calculation!
217Function SmearedLogNormalSphere(s) :FitFunc
218        Struct ResSmearAAOStruct &s
219
220//      the name of your unsmeared model (AAO) is the first argument
221        Smear_Model_20(LogNormalSphere,s.coefW,s.xW,s.yW,s.resW)
222
223        return(0)
224End
225       
226
227
228// normalization is correct, using 3rd moment of lognormal distribution
229//
230Function LogNormal_distr(sig,mu,pt)
231        Variable sig,mu,pt
232       
233        Variable retval
234       
235        retval = (1/ ( sig*pt*sqrt(2*pi)) )*exp( -0.5*((ln(pt) - mu)^2)/sig^2 )
236        return(retval)
237End
238
239//calculates number density given the coefficients of the lognormal distribution
240// the scale factor is the volume fraction
241// then nden = phi/<V> where <V> is calculated using the 3rd moment of the radius
242Macro NumberDensity_LogN()
243       
244        Variable nden,r3,rg,sv,i0,ravg,rpk
245       
246        if(Exists("coef_lns")!=1)
247                abort "You need to plot the unsmeared model first to create the coefficient table"
248        Endif
249       
250        Print "median radius (A) = ",coef_lns[1]
251        Print "sigma = ",coef_lns[2]
252        Print "volume fraction = ",coef_lns[0]
253       
254        r3 = exp(3*ln(coef_lns[1]) + 9/2*coef_lns[2]^2)         // <R^3> directly,[A^3]
255        nden = coef_lns[0]/(4*pi/3*r3)          //nden in 1/A^3
256        ravg = exp(ln(coef_lns[1]) + 0.5*coef_lns[2]^2)
257        rpk = exp(ln(coef_lns[1]) - coef_lns[2]^2)
258        rg = (3./5.)^0.5*exp(ln(coef_lns[1]) + 7.*coef_lns[2]^2)
259        sv = 1.0e8*3*coef_lns[0]*exp(-ln(coef_lns[1]) - 2.5*coef_lns[2]^2)
260        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)
261       
262        Print "number density (A^-3) = ",nden
263        Print "mean radius (A) = ",ravg
264        Print "peak dis. radius (A) = ",rpk
265        Print "Guinier radius (A) = ",rg
266        Print "Interfacial surface area / volume (cm-1) Sv = ",sv
267        Print "Forward cross section (cm-1 sr-1) I(0) = ",i0
268End
269
270// plots the lognormal distribution based on the coefficient values
271// a static calculation, so re-run each time
272//
273Macro PlotLogNormalDistribution()
274
275        variable sig,mu,maxr
276       
277        if(Exists("coef_lns")!=1)
278                abort "You need to plot the unsmeared model first to create the coefficient table"
279        Endif
280        sig=coef_lns[2]
281        mu = ln(coef_lns[1])
282       
283        Make/O/D/N=1000 lognormal_distribution
284        maxr = 5*sig*(1+sig)+ mu
285        maxr = exp(maxr)
286        SetScale/I x, 0, maxr, lognormal_distribution
287        lognormal_distribution = LogNormal_distr(sig,mu,x)
288        Display lognormal_distribution
289        modifygraph log(bottom)=1
290        legend
291End
Note: See TracBrowser for help on using the repository browser.