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

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

Removed Angstrom symbol (Mac code) and replaced it with simply a capital "A" so strange symbols won't appear anymore on Win.

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