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 |
---|
9 | Proc 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) |
---|
26 | End |
---|
27 | |
---|
28 | Proc 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) |
---|
53 | End |
---|
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 | // |
---|
60 | Function 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) |
---|
150 | End |
---|
151 | |
---|
152 | // this is all there is to the smeared calculation! |
---|
153 | Function 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) |
---|
167 | End |
---|
168 | |
---|
169 | // normalization is correct, using 3rd moment of lognormal distribution |
---|
170 | // |
---|
171 | Function 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) |
---|
178 | End |
---|
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 |
---|
183 | Macro 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 |
---|
209 | End |
---|
210 | |
---|
211 | // plots the lognormal distribution based on the coefficient values |
---|
212 | // a static calculation, so re-run each time |
---|
213 | // |
---|
214 | Macro 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 |
---|
232 | End |
---|