1 | #pragma rtGlobals=1 // Use modern global access method. |
---|
2 | #pragma IgorVersion = 6.0 |
---|
3 | |
---|
4 | #include "sphere" |
---|
5 | // plots the form factor of spheres with a Gaussian radius distribution |
---|
6 | // |
---|
7 | // also can plot the distribution itself, based on the current model parameters |
---|
8 | // |
---|
9 | // integration is currently done using 20-pt quadrature, but may benefit from |
---|
10 | //switching to an adaptive integration. |
---|
11 | // |
---|
12 | |
---|
13 | Proc PlotGaussPolySphere(num,qmin,qmax) |
---|
14 | Variable num=128,qmin=0.001,qmax=0.7 |
---|
15 | Prompt num "Enter number of data points for model: " |
---|
16 | Prompt qmin "Enter minimum q-value (^-1) for model: " |
---|
17 | Prompt qmax "Enter maximum q-value (^-1) for model: " |
---|
18 | |
---|
19 | Make/O/D/N=(num) xwave_pgs,ywave_pgs |
---|
20 | xwave_pgs = alog( log(qmin) + x*((log(qmax)-log(qmin))/num) ) |
---|
21 | Make/O/D coef_pgs = {0.01,60,0.2,1e-6,3e-6,0.001} |
---|
22 | make/O/T parameters_pgs = {"Volume Fraction (scale)","mean radius (A)","polydisp (sig/avg)","SLD sphere (A-2)","SLD solvent (A-2)","bkg (cm-1 sr-1)"} |
---|
23 | Edit parameters_pgs,coef_pgs |
---|
24 | |
---|
25 | Variable/G root:g_pgs |
---|
26 | g_pgs := GaussPolySphere(coef_pgs,ywave_pgs,xwave_pgs) |
---|
27 | Display ywave_pgs vs xwave_pgs |
---|
28 | ModifyGraph log=1,marker=29,msize=2,mode=4 |
---|
29 | Label bottom "q (\\S-1\\M)" |
---|
30 | Label left "Intensity (cm\\S-1\\M)" |
---|
31 | AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2) |
---|
32 | End |
---|
33 | |
---|
34 | // - sets up a dependency to a wrapper, not the actual SmearedModelFunction |
---|
35 | Proc PlotSmearedGaussPolySpheres(str) |
---|
36 | String str |
---|
37 | Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4) |
---|
38 | |
---|
39 | // if any of the resolution waves are missing => abort |
---|
40 | if(ResolutionWavesMissingDF(str)) //updated to NOT use global strings (in GaussUtils) |
---|
41 | Abort |
---|
42 | endif |
---|
43 | |
---|
44 | SetDataFolder $("root:"+str) |
---|
45 | |
---|
46 | // Setup parameter table for model function |
---|
47 | Make/O/D smear_coef_pgs = {0.01,60,0.2,1e-6,3e-6,0.001} |
---|
48 | make/o/t smear_parameters_pgs = {"Volume Fraction (scale)","mean radius (A)","polydisp (sig/avg)","SLD sphere (A-2)","SLD solvent (A-2)","bkg (cm-1 sr-1)"} |
---|
49 | Edit smear_parameters_pgs,smear_coef_pgs |
---|
50 | |
---|
51 | // output smeared intensity wave, dimensions are identical to experimental QSIG values |
---|
52 | // make extra copy of experimental q-values for easy plotting |
---|
53 | Duplicate/O $(str+"_q") smeared_pgs,smeared_qvals |
---|
54 | SetScale d,0,0,"1/cm",smeared_pgs |
---|
55 | |
---|
56 | Variable/G gs_pgs=0 |
---|
57 | gs_pgs := fSmearedGaussPolySphere(smear_coef_pgs,smeared_pgs,smeared_qvals) //this wrapper fills the STRUCT |
---|
58 | |
---|
59 | Display smeared_pgs vs smeared_qvals |
---|
60 | ModifyGraph log=1,marker=29,msize=2,mode=4 |
---|
61 | Label bottom "q (\\S-1\\M)" |
---|
62 | Label left "Intensity (cm\\S-1\\M)" |
---|
63 | AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2) |
---|
64 | |
---|
65 | SetDataFolder root: |
---|
66 | End |
---|
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) |
---|
74 | Function GaussPolySphere(cw,yw,xw) : FitFunc |
---|
75 | Wave cw,yw,xw |
---|
76 | |
---|
77 | #if exists("GaussPolySphereX") |
---|
78 | yw = GaussPolySphereX(cw,xw) |
---|
79 | #else |
---|
80 | yw = fGaussPolySphere(cw,xw) |
---|
81 | #endif |
---|
82 | return(0) |
---|
83 | End |
---|
84 | |
---|
85 | Function fGaussPolySphere(w,xx) : FitFunc |
---|
86 | wave w |
---|
87 | variable xx |
---|
88 | |
---|
89 | Variable scale,rad,pd,sig,rho,rhos,bkg,delrho |
---|
90 | |
---|
91 | //set up the coefficient values |
---|
92 | scale=w[0] |
---|
93 | rad=w[1] |
---|
94 | pd=w[2] |
---|
95 | sig=pd*rad |
---|
96 | rho=w[3] |
---|
97 | rhos=w[4] |
---|
98 | delrho=rho-rhos |
---|
99 | bkg=w[5] |
---|
100 | |
---|
101 | //temp set scale=1 and bkg=0 for quadrature calc |
---|
102 | Make/O/D/N=4 sphere_temp |
---|
103 | sphere_temp[0] = 1 |
---|
104 | sphere_temp[1] = rad //changed in loop |
---|
105 | sphere_temp[2] = delrho |
---|
106 | sphere_temp[3] = 0 |
---|
107 | |
---|
108 | //could use 5 pt quadrature to integrate over the size distribution, since it's a gaussian |
---|
109 | //currently using 20 pts... |
---|
110 | Variable va,vb,ii,zi,nord,yy,summ,inten |
---|
111 | String weightStr,zStr |
---|
112 | |
---|
113 | //select number of gauss points by setting nord=20 or76 points |
---|
114 | nord = 20 |
---|
115 | // nord = 76 |
---|
116 | |
---|
117 | weightStr = "gauss"+num2str(nord)+"wt" |
---|
118 | zStr = "gauss"+num2str(nord)+"z" |
---|
119 | |
---|
120 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
121 | Make/D/N=(nord) $weightStr,$zStr |
---|
122 | Wave gauWt = $weightStr |
---|
123 | Wave gauZ = $zStr // wave references to pass |
---|
124 | if(nord==20) |
---|
125 | Make20GaussPoints(gauWt,gauZ) |
---|
126 | else |
---|
127 | Make76GaussPoints(gauWt,gauZ) |
---|
128 | endif |
---|
129 | else |
---|
130 | if(exists(weightStr) > 1) |
---|
131 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
132 | endif |
---|
133 | Wave gauWt = $weightStr |
---|
134 | Wave gauZ = $zStr // create the wave references |
---|
135 | endif |
---|
136 | |
---|
137 | // end points of integration |
---|
138 | // limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero |
---|
139 | // +/- 3 sigq catches 99.73% of distrubution |
---|
140 | // change limits (and spacing of zi) at each evaluation based on R() |
---|
141 | //integration from va to vb |
---|
142 | va = -4*sig + rad |
---|
143 | if (va<0) |
---|
144 | va=0 //to avoid numerical error when va<0 (-ve q-value) |
---|
145 | endif |
---|
146 | vb = 4*sig +rad |
---|
147 | |
---|
148 | summ = 0.0 // initialize integral |
---|
149 | Make/O/N=1 tmp_yw,tmp_xw |
---|
150 | tmp_xw[0] = xx |
---|
151 | for(ii=0;ii<nord;ii+=1) |
---|
152 | // calculate Gauss points on integration interval (r-value for evaluation) |
---|
153 | zi = ( gauZ[ii]*(vb-va) + vb + va )/2.0 |
---|
154 | sphere_temp[1] = zi |
---|
155 | // calculate sphere scattering |
---|
156 | SphereForm(sphere_temp,tmp_yw,tmp_xw) // AAO calculation, one point |
---|
157 | yy = gauWt[ii] * Gauss_distr(sig,rad,zi) * tmp_yw[0] |
---|
158 | yy *= 4*pi/3*zi*zi*zi //un-normalize by current sphere volume |
---|
159 | |
---|
160 | summ += yy //add to the running total of the quadrature |
---|
161 | endfor |
---|
162 | // calculate value of integral to return |
---|
163 | inten = (vb-va)/2.0*summ |
---|
164 | |
---|
165 | //re-normalize by polydisperse sphere volume |
---|
166 | inten /= (4*pi/3*rad^3)*(1+3*pd^2) |
---|
167 | |
---|
168 | inten *= scale |
---|
169 | inten+=bkg |
---|
170 | |
---|
171 | Return(inten) |
---|
172 | End |
---|
173 | |
---|
174 | //wrapper to calculate the smeared model as an AAO-Struct |
---|
175 | // fills the struct and calls the ususal function with the STRUCT parameter |
---|
176 | // |
---|
177 | // used only for the dependency, not for fitting |
---|
178 | // |
---|
179 | Function fSmearedGaussPolySphere(coefW,yW,xW) |
---|
180 | Wave coefW,yW,xW |
---|
181 | |
---|
182 | String str = getWavesDataFolder(yW,0) |
---|
183 | String DF="root:"+str+":" |
---|
184 | |
---|
185 | WAVE resW = $(DF+str+"_res") |
---|
186 | |
---|
187 | STRUCT ResSmearAAOStruct fs |
---|
188 | WAVE fs.coefW = coefW |
---|
189 | WAVE fs.yW = yW |
---|
190 | WAVE fs.xW = xW |
---|
191 | WAVE fs.resW = resW |
---|
192 | |
---|
193 | Variable err |
---|
194 | err = SmearedGaussPolySphere(fs) |
---|
195 | |
---|
196 | return (0) |
---|
197 | End |
---|
198 | |
---|
199 | // this is all there is to the smeared calculation! |
---|
200 | Function SmearedGaussPolySphere(s) :FitFunc |
---|
201 | Struct ResSmearAAOStruct &s |
---|
202 | |
---|
203 | // the name of your unsmeared model (AAO) is the first argument |
---|
204 | s.yW = Smear_Model_20(GaussPolySphere,s.coefW,s.xW,s.resW) |
---|
205 | |
---|
206 | return(0) |
---|
207 | End |
---|
208 | |
---|
209 | |
---|
210 | |
---|
211 | //calculates number density given the coefficients of the Gaussian distribution |
---|
212 | // the scale factor is the volume fraction |
---|
213 | // then nden = phi/<V> where <V> is calculated using the 3rd moment of the radius |
---|
214 | Macro NumberDensity_Gauss() |
---|
215 | |
---|
216 | Variable nden,vpoly,Ravg,p,Rg,I0,Sv,phi |
---|
217 | |
---|
218 | if(WaveExists(coef_pgs)==0) |
---|
219 | abort "You need to plot the model first to create the coefficient table" |
---|
220 | Endif |
---|
221 | Ravg = coef_pgs[1] // mean radius (A) |
---|
222 | p = coef_pgs[2] // polydispersity |
---|
223 | phi = coef_pgs[0] // volume fraction |
---|
224 | Print "mean radius (A) = ",Ravg |
---|
225 | Print "polydispersity (sig/avg) = ",p |
---|
226 | Print "volume fraction = ",phi |
---|
227 | |
---|
228 | //re-normalize by polydisperse sphere volume |
---|
229 | vpoly = (4*pi/3*Ravg^3)*(1+3*p^2) |
---|
230 | nden = phi/vpoly //nden in 1/A^3 |
---|
231 | |
---|
232 | Print "number density (A^-3) = ",nden |
---|
233 | |
---|
234 | Rg = Ravg*(0.6)^0.5*(1+28*p^2+210*p^4+420*p^6+105*p^8) / (1+15*p^2+45*p^4+15*p^6) |
---|
235 | Print "Guinier Radius (A) = ",Rg |
---|
236 | I0 = 1.0e8*(4./3.)*pi*phi*(coef_pgs[3]- coef_pgs[4])^2*Ravg^3*(1+15*p^2+45*p^4+15*p^6)/(1+3*p^2) |
---|
237 | Print "Forward scattering cross-section (cm-1 sr-1) I(0)= ",I0 |
---|
238 | Sv=1.0e8*(3*phi/Ravg)*(1+p^2)/(1+3*p^2) |
---|
239 | Print "Interfacial surface area per unit sample volume (cm-1) Sv= ",Sv |
---|
240 | End |
---|
241 | |
---|
242 | // plots the Gauss distribution based on the coefficient values |
---|
243 | // a static calculation, so re-run each time |
---|
244 | // |
---|
245 | Macro PlotGaussDistribution() |
---|
246 | |
---|
247 | variable pd,avg,zz,maxr |
---|
248 | |
---|
249 | if(WaveExists(coef_pgs)==0) |
---|
250 | abort "You need to plot the model first to create the coefficient table" |
---|
251 | Endif |
---|
252 | pd=coef_pgs[2] |
---|
253 | avg = coef_pgs[1] |
---|
254 | |
---|
255 | Make/O/D/N=1000 Gauss_distribution |
---|
256 | maxr = avg*(1+10*pd) |
---|
257 | |
---|
258 | SetScale/I x, 0, maxr, Gauss_distribution |
---|
259 | Gauss_distribution = Gauss_distr(pd*avg,avg,x) |
---|
260 | Display Gauss_distribution |
---|
261 | legend |
---|
262 | End |
---|
263 | |
---|
264 | |
---|
265 | Function Gauss_distr(sig,avg,pt) |
---|
266 | Variable sig,avg,pt |
---|
267 | |
---|
268 | Variable retval |
---|
269 | |
---|
270 | retval = (1/ ( sig*sqrt(2*Pi)) )*exp(-(avg-pt)^2/sig^2/2) |
---|
271 | return(retval) |
---|
272 | End |
---|