1 | #pragma rtGlobals=1 // Use modern global access method. |
---|
2 | #pragma IgorVersion = 6.0 |
---|
3 | |
---|
4 | //#include "Sphere" |
---|
5 | //plots the form factor for spheres with a Sculz distribution of radius |
---|
6 | // at low polydispersity (< 0.2), it is very similar to the Gaussian distribution |
---|
7 | // at larger polydispersities, it is more skewed and similar to log-normal |
---|
8 | // |
---|
9 | |
---|
10 | // |
---|
11 | Proc PlotSchulzSpheres(num,qmin,qmax) |
---|
12 | Variable num=128,qmin=0.001,qmax=0.7 |
---|
13 | Prompt num "Enter number of data points for model: " |
---|
14 | Prompt qmin "Enter minimum q-value (^-1) for model: " |
---|
15 | Prompt qmax "Enter maximum q-value (^-1) for model: " |
---|
16 | |
---|
17 | Make/O/D/N=(num) xwave_sch,ywave_sch |
---|
18 | xwave_sch = alog( log(qmin) + x*((log(qmax)-log(qmin))/num) ) |
---|
19 | Make/O/D coef_sch = {0.01,60,0.2,1e-6,3e-6,0.001} |
---|
20 | make/O/T parameters_sch = {"Volume Fraction (scale)","mean radius (A)","polydisp (sig/avg)","SLD sphere (A-2)","SLD solvent (A-2)","bkg (cm-1 sr-1)"} |
---|
21 | Edit parameters_sch,coef_sch |
---|
22 | |
---|
23 | Variable/G root:g_sch |
---|
24 | g_sch := SchulzSpheres(coef_sch,ywave_sch,xwave_sch) |
---|
25 | Display ywave_sch vs xwave_sch |
---|
26 | ModifyGraph log=1,marker=29,msize=2,mode=4 |
---|
27 | Label bottom "q (\\S-1\\M)" |
---|
28 | Label left "Intensity (cm\\S-1\\M)" |
---|
29 | AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2) |
---|
30 | End |
---|
31 | |
---|
32 | |
---|
33 | // - sets up a dependency to a wrapper, not the actual SmearedModelFunction |
---|
34 | Proc PlotSmearedSchulzSpheres(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_sch = {0.01,60,0.2,1e-6,3e-6,0.001} |
---|
47 | make/o/t smear_parameters_sch = {"Volume Fraction (scale)","mean radius (A)","polydisp (sig/avg)","SLD sphere (A-2)","SLD solvent (A-2)","bkg (cm-1 sr-1)"} |
---|
48 | Edit smear_parameters_sch,smear_coef_sch |
---|
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_sch,smeared_qvals |
---|
53 | SetScale d,0,0,"1/cm",smeared_sch |
---|
54 | |
---|
55 | Variable/G gs_sch=0 |
---|
56 | gs_sch := fSmearedSchulzSpheres(smear_coef_sch,smeared_sch,smeared_qvals) //this wrapper fills the STRUCT |
---|
57 | |
---|
58 | Display smeared_sch vs smeared_qvals |
---|
59 | ModifyGraph log=1,marker=29,msize=2,mode=4 |
---|
60 | Label bottom "q (\\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 | End |
---|
66 | |
---|
67 | |
---|
68 | |
---|
69 | Function Schulz_Point(x,avg,zz) |
---|
70 | Variable x,avg,zz |
---|
71 | |
---|
72 | Variable dr |
---|
73 | |
---|
74 | dr = zz*ln(x) - gammln(zz+1)+(zz+1)*ln((zz+1)/avg)-(x/avg*(zz+1)) |
---|
75 | return (exp(dr)) |
---|
76 | End |
---|
77 | |
---|
78 | |
---|
79 | //AAO version, uses XOP if available |
---|
80 | // simply calls the original single point calculation with |
---|
81 | // a wave assignment (this will behave nicely if given point ranges) |
---|
82 | Function SchulzSpheres(cw,yw,xw) : FitFunc |
---|
83 | Wave cw,yw,xw |
---|
84 | |
---|
85 | #if exists("SchulzSpheresX") |
---|
86 | yw = SchulzSpheresX(cw,xw) |
---|
87 | #else |
---|
88 | yw = fSchulzSpheres(cw,xw) |
---|
89 | #endif |
---|
90 | return(0) |
---|
91 | End |
---|
92 | |
---|
93 | //use the analytic formula from Kotlarchyk & Chen, JCP 79 (1983) 2461 |
---|
94 | //equations 23-30 |
---|
95 | // |
---|
96 | // need to calculate in terms of logarithms to avoid numerical errors |
---|
97 | // |
---|
98 | Function fSchulzSpheres(w,x) : FitFunc |
---|
99 | Wave w |
---|
100 | Variable x |
---|
101 | |
---|
102 | Variable scale,ravg,pd,delrho,bkg,zz,rho,rhos,vpoly |
---|
103 | scale = w[0] |
---|
104 | ravg = w[1] |
---|
105 | pd = w[2] |
---|
106 | rho = w[3] |
---|
107 | rhos = w[4] |
---|
108 | bkg = w[5] |
---|
109 | |
---|
110 | delrho=rho-rhos |
---|
111 | zz = (1/pd)^2-1 |
---|
112 | |
---|
113 | Variable zp1,zp2,zp3,zp4,zp5,zp6,zp7 |
---|
114 | Variable aa,b1,b2,b3,at1,at2,rt1,rt2,rt3,t1,t2,t3 |
---|
115 | Variable v1,v2,v3,g1,g11,gd,pq,g2,g22 |
---|
116 | |
---|
117 | ZP1 = zz + 1 |
---|
118 | ZP2 = zz + 2 |
---|
119 | ZP3 = zz + 3 |
---|
120 | ZP4 = zz + 4 |
---|
121 | ZP5 = zz + 5 |
---|
122 | ZP6 = zz + 6 |
---|
123 | ZP7 = zz + 7 |
---|
124 | |
---|
125 | //small QR limit - use Guinier approx |
---|
126 | Variable i_zero,Rg2,zp8 |
---|
127 | zp8 = zz+8 |
---|
128 | if(x*ravg < 0.1) |
---|
129 | i_zero = scale*delrho*delrho*1e8*4*Pi/3*ravg^3 |
---|
130 | i_zero *= zp6*zp5*zp4/zp1/zp1/zp1 //6th moment / 3rd moment |
---|
131 | Rg2 = 3*zp8*zp7/5/(zp1^2)*ravg*ravg |
---|
132 | pq = i_zero*exp(-x*x*Rg2/3) |
---|
133 | pq += bkg |
---|
134 | return(pq) |
---|
135 | endif |
---|
136 | // |
---|
137 | aa = (zz+1)/x/Ravg |
---|
138 | |
---|
139 | AT1 = atan(1/aa) |
---|
140 | AT2 = atan(2/aa) |
---|
141 | // |
---|
142 | // CALCULATIONS ARE PERFORMED TO AVOID LARGE # ERRORS |
---|
143 | // - trick is to propogate the a^(z+7) term through the G1 |
---|
144 | // |
---|
145 | T1 = ZP7*log(aa) - zp1/2*log(aa*aa+4) |
---|
146 | T2 = ZP7*log(aa) - zp3/2*log(aa*aa+4) |
---|
147 | T3 = ZP7*log(aa) - zp2/2*log(aa*aa+4) |
---|
148 | // Print T1,T2,T3 |
---|
149 | RT1 = alog(T1) |
---|
150 | RT2 = alog(T2) |
---|
151 | RT3 = alog(T3) |
---|
152 | V1 = aa^6 - RT1*cos(zp1*at2) |
---|
153 | V2 = ZP1*ZP2*( aa^4 + RT2*cos(zp3*at2) ) |
---|
154 | V3 = -2*ZP1*RT3*SIN(zp2*at2) |
---|
155 | G1 = (V1+V2+V3) |
---|
156 | |
---|
157 | Pq = log(G1) - 6*log(ZP1) + 6*log(Ravg) |
---|
158 | Pq = alog(Pq)*8*PI*PI*delrho*delrho |
---|
159 | |
---|
160 | // |
---|
161 | // beta factor is not used here, but could be for the |
---|
162 | // decoupling approximation |
---|
163 | // |
---|
164 | // G11 = G1 |
---|
165 | // GD = -ZP7*log(aa) |
---|
166 | // G1 = log(G11) + GD |
---|
167 | // |
---|
168 | // T1 = ZP1*at1 |
---|
169 | // T2 = ZP2*at1 |
---|
170 | // G2 = SIN( T1 ) - ZP1/SQRT(aa*aa+1)*COS( T2 ) |
---|
171 | // G22 = G2*G2 |
---|
172 | // BETA = ZP1*log(aa) - ZP1*log(aa*aa+1) - G1 + log(G22) |
---|
173 | // BETA = 2*alog(BETA) |
---|
174 | |
---|
175 | //re-normalize by the average volume |
---|
176 | vpoly = 4*Pi/3*zp3*zp2/zp1/zp1*(ravg)^3 |
---|
177 | Pq /= vpoly |
---|
178 | //scale, convert to cm^-1 |
---|
179 | Pq *= scale * 1e8 |
---|
180 | // add in the background |
---|
181 | Pq += bkg |
---|
182 | |
---|
183 | //return (g1) |
---|
184 | Return (Pq) |
---|
185 | End |
---|
186 | |
---|
187 | //wrapper to calculate the smeared model as an AAO-Struct |
---|
188 | // fills the struct and calls the ususal function with the STRUCT parameter |
---|
189 | // |
---|
190 | // used only for the dependency, not for fitting |
---|
191 | // |
---|
192 | Function fSmearedSchulzSpheres(coefW,yW,xW) |
---|
193 | Wave coefW,yW,xW |
---|
194 | |
---|
195 | String str = getWavesDataFolder(yW,0) |
---|
196 | String DF="root:"+str+":" |
---|
197 | |
---|
198 | WAVE resW = $(DF+str+"_res") |
---|
199 | |
---|
200 | STRUCT ResSmearAAOStruct fs |
---|
201 | WAVE fs.coefW = coefW |
---|
202 | WAVE fs.yW = yW |
---|
203 | WAVE fs.xW = xW |
---|
204 | WAVE fs.resW = resW |
---|
205 | |
---|
206 | Variable err |
---|
207 | err = SmearedSchulzSpheres(fs) |
---|
208 | |
---|
209 | return (0) |
---|
210 | End |
---|
211 | |
---|
212 | // this is all there is to the smeared calculation! |
---|
213 | Function SmearedSchulzSpheres(s) :FitFunc |
---|
214 | Struct ResSmearAAOStruct &s |
---|
215 | |
---|
216 | // the name of your unsmeared model (AAO) is the first argument |
---|
217 | Smear_Model_20(SchulzSpheres,s.coefW,s.xW,s.yW,s.resW) |
---|
218 | |
---|
219 | return(0) |
---|
220 | End |
---|
221 | |
---|
222 | |
---|
223 | //calculates number density given the coefficients of the Schulz distribution |
---|
224 | // the scale factor is the volume fraction |
---|
225 | // then nden = phi/<V> where <V> is calculated using the 3rd moment of the radius |
---|
226 | Macro NumberDensity_Schulz() |
---|
227 | |
---|
228 | Variable nden,zz,zp1,zp2,zp3,zp4,zp5,zp6,zp7,zp8,vpoly,v2poly,rg,sv,i0 |
---|
229 | |
---|
230 | if(WaveExists(coef_sch)==0) |
---|
231 | abort "You need to plot the model first to create the coefficient table" |
---|
232 | Endif |
---|
233 | |
---|
234 | Print "mean radius (A) = ",coef_sch[1] |
---|
235 | Print "polydispersity (sig/avg) = ",coef_sch[2] |
---|
236 | Print "volume fraction = ",coef_sch[0] |
---|
237 | |
---|
238 | zz = (1/coef_sch[2])^2-1 |
---|
239 | zp1 = zz + 1. |
---|
240 | zp2 = zz + 2. |
---|
241 | zp3 = zz + 3. |
---|
242 | zp4 = zz + 4. |
---|
243 | zp5 = zz + 5. |
---|
244 | zp6 = zz + 6. |
---|
245 | zp7 = zz + 7. |
---|
246 | zp8 = zz + 8. |
---|
247 | // average particle volume |
---|
248 | vpoly = 4*Pi/3*zp3*zp2/zp1/zp1*(coef_sch[1])^3 |
---|
249 | // average particle volume |
---|
250 | v2poly = (4*Pi/3)^2*zp6*zp5*zp4*zp3*zp2/(zp1^5)*(coef_sch[1])^6 |
---|
251 | nden = coef_sch[0]/vpoly //nden in 1/A^3 |
---|
252 | rg = coef_sch[1]*((3*zp8*zp7)/5/zp1/zp1)^0.5 // in A |
---|
253 | sv = 1.0e8*3*coef_sch[0]*zp1/coef_sch[1]/zp3 // in 1/cm |
---|
254 | i0 = 1.0e8*nden*v2poly*(coef_sch[3]-coef_sch[4])^2 // 1/cm/sr |
---|
255 | |
---|
256 | Print "number density (A^-3) = ",nden |
---|
257 | Print "Guinier radius (A) = ",rg |
---|
258 | Print "Interfacial surface area / volume (cm-1) Sv = ",sv |
---|
259 | Print "Forward cross section (cm-1 sr-1) I(0) = ",i0 |
---|
260 | End |
---|
261 | |
---|
262 | // plots the Schulz distribution based on the coefficient values |
---|
263 | // a static calculation, so re-run each time |
---|
264 | // |
---|
265 | Macro PlotSchulzDistribution() |
---|
266 | |
---|
267 | variable pd,avg,zz,maxr |
---|
268 | |
---|
269 | if(WaveExists(coef_sch)==0) |
---|
270 | abort "You need to plot the model first to create the coefficient table" |
---|
271 | Endif |
---|
272 | pd=coef_sch[2] |
---|
273 | avg = coef_sch[1] |
---|
274 | zz = (1/pd)^2-1 |
---|
275 | |
---|
276 | Make/O/D/N=1000 Schulz_distribution |
---|
277 | maxr = avg*(1+10*pd) |
---|
278 | |
---|
279 | SetScale/I x, 0, maxr, Schulz_distribution |
---|
280 | Schulz_distribution = Schulz_Point(x,avg,zz) |
---|
281 | Display Schulz_distribution |
---|
282 | legend |
---|
283 | End |
---|
284 | |
---|
285 | // don't use the integral technique here since there's an analytic solution |
---|
286 | // available |
---|
287 | ////////////////////// |
---|
288 | //requires that sphere.ipf is included |
---|
289 | // |
---|
290 | // |
---|
291 | //Function SchulzSpheres_Integrated(w,x) |
---|
292 | // Wave w |
---|
293 | // Variable x |
---|
294 | // |
---|
295 | // Variable scale,radius,pd,delrho,bkg,zz,rho,rhos |
---|
296 | // scale = w[0] |
---|
297 | // radius = w[1] |
---|
298 | // pd = w[2] |
---|
299 | // rho = w[3] |
---|
300 | // rhos = w[4] |
---|
301 | // bkg = w[5] |
---|
302 | // |
---|
303 | // delrho=rho-rhos |
---|
304 | // zz = (1/pd)^2-1 |
---|
305 | //// |
---|
306 | //// local variables |
---|
307 | // Variable nord,ii,a,b,va,vb,contr,vcyl,nden,summ,yyy,zi,qq |
---|
308 | // Variable answer,zp1,zp2,zp3,vpoly |
---|
309 | // String weightStr,zStr |
---|
310 | // |
---|
311 | // //select number of gauss points by setting nord=20 or76 points |
---|
312 | //// nord = 20 |
---|
313 | // nord = 76 |
---|
314 | // |
---|
315 | // weightStr = "gauss"+num2str(nord)+"wt" |
---|
316 | // zStr = "gauss"+num2str(nord)+"z" |
---|
317 | // |
---|
318 | // if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
319 | // Make/D/N=(nord) $weightStr,$zStr |
---|
320 | // Wave gauWt = $weightStr |
---|
321 | // Wave gauZ = $zStr // wave references to pass |
---|
322 | // if(nord==20) |
---|
323 | // Make20GaussPoints(gauWt,gauZ) |
---|
324 | // else |
---|
325 | // Make76GaussPoints(gauWt,gauZ) |
---|
326 | // endif |
---|
327 | // else |
---|
328 | // if(exists(weightStr) > 1) |
---|
329 | // Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
330 | // endif |
---|
331 | // Wave gauWt = $weightStr |
---|
332 | // Wave gauZ = $zStr // create the wave references |
---|
333 | // endif |
---|
334 | // |
---|
335 | //// set up the integration |
---|
336 | //// end points and weights |
---|
337 | //// limits are technically 0-inf, but wisely choose non-zero region of distribution |
---|
338 | // Variable range=8 //multiples of the std. dev. from the mean |
---|
339 | // a = radius*(1-range*pd) |
---|
340 | // if (a<0) |
---|
341 | // a=0 //otherwise numerical error when pd >= 0.3, making a<0 |
---|
342 | // endif |
---|
343 | // If(pd>0.3) |
---|
344 | // range = 3.4 + (pd-0.3)*18 //to account for skewed tail |
---|
345 | // Endif |
---|
346 | // b = radius*(1+range*pd) // is this far enough past avg radius? |
---|
347 | // va =a |
---|
348 | // vb =b |
---|
349 | // |
---|
350 | ////temp set scale=1 and bkg=0 for quadrature calc |
---|
351 | // Make/O/D/N=4 sphere_temp |
---|
352 | // sphere_temp[0] = 1 |
---|
353 | // sphere_temp[1] = radius //changed in loop |
---|
354 | // sphere_temp[2] = delrho |
---|
355 | // sphere_temp[3] = 0 |
---|
356 | // |
---|
357 | //// evaluate at Gauss points |
---|
358 | // summ = 0.0 // initialize integral |
---|
359 | // for(ii=0;ii<nord;ii+=1) |
---|
360 | // zi = ( gauZ[ii]*(vb-va) + vb + va )/2.0 |
---|
361 | // sphere_temp[1] = zi |
---|
362 | // yyy = gauWt[ii] * Schulz_Point(zi,radius,zz) * SphereForm(sphere_temp,x) |
---|
363 | // //un-normalize by volume |
---|
364 | // yyy *= 4*pi/3*zi^3 |
---|
365 | // summ += yyy |
---|
366 | // endfor |
---|
367 | //// calculate value of integral to return |
---|
368 | // answer = (vb-va)/2.0*summ |
---|
369 | // |
---|
370 | // //re-normalize by the average volume |
---|
371 | // zp1 = zz + 1. |
---|
372 | // zp2 = zz + 2. |
---|
373 | // zp3 = zz + 3. |
---|
374 | // vpoly = 4*Pi/3*zp3*zp2/zp1/zp1*(radius)^3 |
---|
375 | // answer /= vpoly |
---|
376 | ////scale |
---|
377 | // answer *= scale |
---|
378 | //// add in the background |
---|
379 | // answer += bkg |
---|
380 | // |
---|
381 | // Return (answer) |
---|
382 | //End |
---|
383 | // |
---|