source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/NewModels_2006/LogNormalSphere.ipf @ 92

Last change on this file since 92 was 42, checked in by srkline, 16 years ago

initial checkin of Analysis v.3.00 files

File size: 7.1 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2
3#include "sphere"
4// plots the form factor of spherical particles with a log-normal radius distribution
5//
6// for the integration it may be better to use adaptive routine
7
8//Proc to setup data and coefficients to plot the model
9Proc PlotLogNormalPolySphere(num,qmin,qmax)
10        Variable num=128,qmin=0.001,qmax=0.7
11        Prompt num "Enter number of data points for model: "
12        Prompt qmin "Enter minimum q-value (^-1) for model: "
13        Prompt qmax "Enter maximum q-value (^-1) for model: "
14       
15        Make/O/D/N=(num) xwave_lns,ywave_lns
16        xwave_lns = alog( log(qmin) + x*((log(qmax)-log(qmin))/num) )
17        Make/O/D coef_lns = {0.01,60,0.2,1e-6,2e-6,0}
18        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)"}
19        Edit parameters_lns,coef_lns
20        ywave_lns := LogNormalPolySphere(coef_lns,xwave_lns)
21        Display ywave_lns vs xwave_lns
22        ModifyGraph log=1,marker=29,msize=2,mode=4
23        Label bottom "q (\\S-1\\M)"
24        Label left "Intensity (cm\\S-1\\M)"
25        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
26End
27
28Proc PlotSmearLogNormPolySphere()                                                               
29        //no input parameters necessary, it MUST use the experimental q-values
30        // from the experimental data read in from an AVE/QSIG data file
31       
32        // if no gQvals wave, data must not have been loaded => abort
33        if(ResolutionWavesMissing())
34                Abort
35        endif
36       
37        // Setup parameter table for model function
38        Make/O/D smear_coef_lns = {0.01,60,0.2,1e-6,2e-6,0}                                     
39        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)"}               
40        Edit smear_parameters_lns,smear_coef_lns                                       
41       
42        // output smeared intensity wave, dimensions are identical to experimental QSIG values
43        // make extra copy of experimental q-values for easy plotting
44        Duplicate/O $gQvals smeared_lns,smeared_qvals                           
45        SetScale d,0,0,"1/cm",smeared_lns                                                       
46
47        smeared_lns := SmearLogNormSphere(smear_coef_lns,$gQvals)               
48        Display smeared_lns vs smeared_qvals                                                                   
49        ModifyGraph log=1,marker=29,msize=2,mode=4
50        Label bottom "q (\\S-1\\M)"
51        Label left "Intensity (cm\\S-1\\M)"
52        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
53End
54
55
56// calculates the model at each q-value by integrating over the normalized size distribution
57// integration is done by gauss quadrature of either 20 or 76 points (nord)
58// 76 points is slower, but reccommended to remove high-q oscillations
59//
60Function LogNormalPolySphere(w,x): FitFunc
61        wave w
62        variable x
63       
64        Variable scale,rad,sig,rho,rhos,bkg,delrho,mu,r3
65       
66        //set up the coefficient values
67        scale=w[0]
68        rad=w[1]                //rad is the median radius
69        mu = ln(w[1])
70        sig=w[2]
71        rho=w[3]
72        rhos=w[4]
73        delrho=rho-rhos
74        bkg=w[5]
75       
76//temp set scale=1 and bkg=0 for quadrature calc
77        Make/O/D/N=4 sphere_temp
78        sphere_temp[0] = 1
79        sphere_temp[1] = rad            //changed in loop
80        sphere_temp[2] = delrho
81        sphere_temp[3] = 0
82       
83        //currently using 20 pts...
84        Variable va,vb,ii,zi,nord,yy,summ,inten
85        String weightStr,zStr
86       
87        //select number of gauss points by setting nord=20 or76 points
88//      nord = 20
89        nord = 76
90       
91        weightStr = "gauss"+num2str(nord)+"wt"
92        zStr = "gauss"+num2str(nord)+"z"
93       
94        if (WaveExists($weightStr) == 0) // wave reference is not valid,
95                Make/D/N=(nord) $weightStr,$zStr
96                Wave gauWt = $weightStr
97                Wave gauZ = $zStr               // wave references to pass
98                if(nord==20)
99                        Make20GaussPoints(gauWt,gauZ)
100                else
101                        Make76GaussPoints(gauWt,gauZ)
102                endif   
103        else
104                if(exists(weightStr) > 1)
105                         Abort "wave name is already in use"            //executed only if name is in use elsewhere
106                endif
107                Wave gauWt = $weightStr
108                Wave gauZ = $zStr               // create the wave references
109        endif
110       
111        // end points of integration
112        // limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero
113        // +/- 3 sigq catches 99.73% of distrubution
114        // change limits (and spacing of zi) at each evaluation based on R()
115        //integration from va to vb
116       
117//      va = -3*sig + rad
118        va = -3.5*sig +mu               //in ln(R) space
119        va = exp(va)                    //in R space
120        if (va<0)
121                va=0            //to avoid numerical error when  va<0 (-ve q-value)
122        endif
123//      vb = 3*exp(sig) +rad
124        vb = 3.5*sig*(1+sig)+ mu
125        vb = exp(vb)
126       
127        summ = 0.0              // initialize integral
128        for(ii=0;ii<nord;ii+=1)
129                // calculate Gauss points on integration interval (r-value for evaluation)
130                zi = ( gauZ[ii]*(vb-va) + vb + va )/2.0
131                sphere_temp[1] = zi
132                // calculate sphere scattering
133                yy = gauWt[ii] *  LogNormal_distr(sig,mu,zi) * SphereForm(sphere_temp,x)
134                yy *= 4*pi/3*zi*zi*zi           //un-normalize by current sphere volume
135               
136                summ += yy              //add to the running total of the quadrature
137        endfor
138// calculate value of integral to return
139        inten = (vb-va)/2.0*summ
140       
141        //re-normalize by polydisperse sphere volume
142        //third moment
143        r3 = exp(3*mu + 9/2*sig^2)              // <R^3> directly
144        inten /= (4*pi/3*r3)            //polydisperse volume
145       
146        inten *= scale
147        inten+=bkg
148       
149        Return(inten)
150End
151
152// this is all there is to the smeared calculation!
153Function SmearLogNormSphere(w,x) :FitFunc
154        Wave w
155        Variable x
156       
157        Variable ans
158        SVAR sq = gSig_Q
159        SVAR qb = gQ_bar
160        SVAR sh = gShadow
161        SVAR gQ = gQVals
162       
163        //the name of your unsmeared model is the first argument
164        ans = Smear_Model_20(LogNormalPolySphere,$sq,$qb,$sh,$gQ,w,x)
165
166        return(ans)
167End
168
169// normalization is correct, using 3rd moment of lognormal distribution
170//
171Function LogNormal_distr(sig,mu,pt)
172        Variable sig,mu,pt
173       
174        Variable retval
175       
176        retval = (1/ ( sig*pt*sqrt(2*pi)) )*exp( -0.5*((ln(pt) - mu)^2)/sig^2 )
177        return(retval)
178End
179
180//calculates number density given the coefficients of the lognormal distribution
181// the scale factor is the volume fraction
182// then nden = phi/<V> where <V> is calculated using the 3rd moment of the radius
183Macro NumberDensity_LogN()
184       
185        Variable nden,r3,rg,sv,i0,ravg,rpk
186       
187        if(WaveExists(coef_lns)==0)
188                abort "You need to plot the model first to create the coefficient table"
189        Endif
190       
191        Print "median radius (A) = ",coef_lns[1]
192        Print "sigma = ",coef_lns[2]
193        Print "volume fraction = ",coef_lns[0]
194       
195        r3 = exp(3*ln(coef_lns[1]) + 9/2*coef_lns[2]^2)         // <R^3> directly,[A^3]
196        nden = coef_lns[0]/(4*pi/3*r3)          //nden in 1/A^3
197        ravg = exp(ln(coef_lns[1]) + 0.5*coef_lns[2]^2)
198        rpk = exp(ln(coef_lns[1]) - coef_lns[2]^2)
199        rg = (3./5.)^0.5*exp(ln(coef_lns[1]) + 7.*coef_lns[2]^2)
200        sv = 1.0e8*3*coef_lns[0]*exp(-ln(coef_lns[1]) - 2.5*coef_lns[2]^2)
201        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)
202       
203        Print "number density (A^-3) = ",nden
204        Print "mean radius (A) = ",ravg
205        Print "peak dis. radius (A) = ",rpk
206        Print "Guinier radius (A) = ",rg
207        Print "Interfacial surface area / volume (cm-1) Sv = ",sv
208        Print "Forward cross section (cm-1 sr-1) I(0) = ",i0
209End
210
211// plots the lognormal distribution based on the coefficient values
212// a static calculation, so re-run each time
213//
214Macro PlotLogNormalDistribution()
215
216        variable sig,mu,maxr
217       
218        if(WaveExists(coef_lns)==0)
219                abort "You need to plot the model first to create the coefficient table"
220        Endif
221        sig=coef_lns[2]
222        mu = ln(coef_lns[1])
223       
224        Make/O/D/N=1000 lognormal_distribution
225        maxr = 5*sig*(1+sig)+ mu
226        maxr = exp(maxr)
227        SetScale/I x, 0, maxr, lognormal_distribution
228        lognormal_distribution = LogNormal_distr(sig,mu,x)
229        Display lognormal_distribution
230        modifygraph log(bottom)=1
231        legend
232End
Note: See TracBrowser for help on using the repository browser.