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

Last change on this file since 1028 was 1028, checked in by srkline, 6 years ago

adding marquee operations to do a simple box sum, for manual calculation of transmission

Error in LamellarParacrystalline? model has been fixed.

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/D/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.