1 | #pragma rtGlobals=1 // Use modern global access method. |
---|
2 | #pragma IgorVersion=6.1 |
---|
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 | // JAN 2012 |
---|
11 | //-- added the functions to calculate the Beta factor for the decoupling approximation to |
---|
12 | // modify the S(q) calculation. A separate Sq_Beta.ipf file was created to hold all of the new P*S' combinations |
---|
13 | // |
---|
14 | |
---|
15 | // |
---|
16 | Proc PlotSchulzSpheres(num,qmin,qmax) |
---|
17 | Variable num=128,qmin=0.001,qmax=0.7 |
---|
18 | Prompt num "Enter number of data points for model: " |
---|
19 | Prompt qmin "Enter minimum q-value (A^-1) for model: " |
---|
20 | Prompt qmax "Enter maximum q-value (A^-1) for model: " |
---|
21 | |
---|
22 | Make/O/D/N=(num) xwave_sch,ywave_sch |
---|
23 | xwave_sch = alog( log(qmin) + x*((log(qmax)-log(qmin))/num) ) |
---|
24 | Make/O/D coef_sch = {0.01,60,0.2,1e-6,3e-6,0.001} |
---|
25 | 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)"} |
---|
26 | Edit parameters_sch,coef_sch |
---|
27 | |
---|
28 | Variable/G root:g_sch |
---|
29 | g_sch := SchulzSpheres(coef_sch,ywave_sch,xwave_sch) |
---|
30 | Display ywave_sch vs xwave_sch |
---|
31 | ModifyGraph log=1,marker=29,msize=2,mode=4 |
---|
32 | Label bottom "q (A\\S-1\\M)" |
---|
33 | Label left "Intensity (cm\\S-1\\M)" |
---|
34 | AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2) |
---|
35 | |
---|
36 | AddModelToStrings("SchulzSpheres","coef_sch","parameters_sch","sch") |
---|
37 | End |
---|
38 | |
---|
39 | |
---|
40 | // - sets up a dependency to a wrapper, not the actual SmearedModelFunction |
---|
41 | Proc PlotSmearedSchulzSpheres(str) |
---|
42 | String str |
---|
43 | Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4) |
---|
44 | |
---|
45 | // if any of the resolution waves are missing => abort |
---|
46 | if(ResolutionWavesMissingDF(str)) //updated to NOT use global strings (in GaussUtils) |
---|
47 | Abort |
---|
48 | endif |
---|
49 | |
---|
50 | SetDataFolder $("root:"+str) |
---|
51 | |
---|
52 | // Setup parameter table for model function |
---|
53 | Make/O/D smear_coef_sch = {0.01,60,0.2,1e-6,3e-6,0.001} |
---|
54 | 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)"} |
---|
55 | Edit smear_parameters_sch,smear_coef_sch |
---|
56 | |
---|
57 | // output smeared intensity wave, dimensions are identical to experimental QSIG values |
---|
58 | // make extra copy of experimental q-values for easy plotting |
---|
59 | Duplicate/O $(str+"_q") smeared_sch,smeared_qvals |
---|
60 | SetScale d,0,0,"1/cm",smeared_sch |
---|
61 | |
---|
62 | Variable/G gs_sch=0 |
---|
63 | gs_sch := fSmearedSchulzSpheres(smear_coef_sch,smeared_sch,smeared_qvals) //this wrapper fills the STRUCT |
---|
64 | |
---|
65 | Display smeared_sch vs smeared_qvals |
---|
66 | ModifyGraph log=1,marker=29,msize=2,mode=4 |
---|
67 | Label bottom "q (A\\S-1\\M)" |
---|
68 | Label left "Intensity (cm\\S-1\\M)" |
---|
69 | AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2) |
---|
70 | |
---|
71 | AddModelToStrings("SmearedSchulzSpheres","smear_coef_sch","smear_parameters_sch","sch") |
---|
72 | End |
---|
73 | |
---|
74 | // |
---|
75 | Proc PlotSchulzSpheresBeta(num,qmin,qmax) |
---|
76 | Variable num=128,qmin=0.001,qmax=0.7 |
---|
77 | Prompt num "Enter number of data points for model: " |
---|
78 | Prompt qmin "Enter minimum q-value (A^-1) for model: " |
---|
79 | Prompt qmax "Enter maximum q-value (A^-1) for model: " |
---|
80 | |
---|
81 | Make/O/D/N=(num) xwave_sch_beta,ywave_sch_beta |
---|
82 | xwave_sch_beta = alog( log(qmin) + x*((log(qmax)-log(qmin))/num) ) |
---|
83 | Make/O/D coef_sch_beta = {0.01,60,0.2,1e-6,3e-6,0.001} |
---|
84 | make/O/T parameters_sch_beta = {"Volume Fraction (scale)","mean radius (A)","polydisp (sig/avg)","SLD sphere (A-2)","SLD solvent (A-2)","bkg (cm-1 sr-1)"} |
---|
85 | Edit parameters_sch_beta,coef_sch_beta |
---|
86 | |
---|
87 | Variable/G root:g_sch_beta |
---|
88 | g_sch_beta := SchulzSpheresBeta(coef_sch_beta,ywave_sch_beta,xwave_sch_beta) |
---|
89 | Display ywave_sch_beta vs xwave_sch_beta |
---|
90 | ModifyGraph log=1,marker=29,msize=2,mode=4 |
---|
91 | Label bottom "q (A\\S-1\\M)" |
---|
92 | Label left "Intensity (cm\\S-1\\M)" |
---|
93 | AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2) |
---|
94 | |
---|
95 | AddModelToStrings("SchulzSpheresBeta","coef_sch_beta","parameters_sch_beta","sch_beta") |
---|
96 | End |
---|
97 | |
---|
98 | Function Schulz_Point(x,avg,zz) |
---|
99 | Variable x,avg,zz |
---|
100 | |
---|
101 | Variable dr |
---|
102 | |
---|
103 | dr = zz*ln(x) - gammln(zz+1)+(zz+1)*ln((zz+1)/avg)-(x/avg*(zz+1)) |
---|
104 | return (exp(dr)) |
---|
105 | End |
---|
106 | |
---|
107 | |
---|
108 | //AAO version, uses XOP if available |
---|
109 | // simply calls the original single point calculation with |
---|
110 | // a wave assignment (this will behave nicely if given point ranges) |
---|
111 | Function SchulzSpheres(cw,yw,xw) : FitFunc |
---|
112 | Wave cw,yw,xw |
---|
113 | |
---|
114 | #if exists("SchulzSpheresX") |
---|
115 | yw = SchulzSpheresX(cw,xw) |
---|
116 | #else |
---|
117 | yw = fSchulzSpheres(cw,xw) |
---|
118 | #endif |
---|
119 | return(0) |
---|
120 | End |
---|
121 | |
---|
122 | //use the analytic formula from Kotlarchyk & Chen, JCP 79 (1983) 2461 |
---|
123 | //equations 23-30 |
---|
124 | // |
---|
125 | // need to calculate in terms of logarithms to avoid numerical errors |
---|
126 | // |
---|
127 | Function fSchulzSpheres(w,x) : FitFunc |
---|
128 | Wave w |
---|
129 | Variable x |
---|
130 | |
---|
131 | Variable scale,ravg,pd,delrho,bkg,zz,rho,rhos,vpoly |
---|
132 | scale = w[0] |
---|
133 | ravg = w[1] |
---|
134 | pd = w[2] |
---|
135 | rho = w[3] |
---|
136 | rhos = w[4] |
---|
137 | bkg = w[5] |
---|
138 | |
---|
139 | delrho=rho-rhos |
---|
140 | zz = (1/pd)^2-1 |
---|
141 | |
---|
142 | Variable zp1,zp2,zp3,zp4,zp5,zp6,zp7 |
---|
143 | Variable aa,b1,b2,b3,at1,at2,rt1,rt2,rt3,t1,t2,t3 |
---|
144 | Variable v1,v2,v3,g1,g11,gd,pq,g2,g22 |
---|
145 | |
---|
146 | ZP1 = zz + 1 |
---|
147 | ZP2 = zz + 2 |
---|
148 | ZP3 = zz + 3 |
---|
149 | ZP4 = zz + 4 |
---|
150 | ZP5 = zz + 5 |
---|
151 | ZP6 = zz + 6 |
---|
152 | ZP7 = zz + 7 |
---|
153 | |
---|
154 | //small QR limit - use Guinier approx |
---|
155 | Variable i_zero,Rg2,zp8 |
---|
156 | zp8 = zz+8 |
---|
157 | if(x*ravg < 0.1) |
---|
158 | i_zero = scale*delrho*delrho*1e8*4*Pi/3*ravg^3 |
---|
159 | i_zero *= zp6*zp5*zp4/zp1/zp1/zp1 //6th moment / 3rd moment |
---|
160 | Rg2 = 3*zp8*zp7/5/(zp1^2)*ravg*ravg |
---|
161 | pq = i_zero*exp(-x*x*Rg2/3) |
---|
162 | pq += bkg |
---|
163 | return(pq) |
---|
164 | endif |
---|
165 | // |
---|
166 | aa = (zz+1)/x/Ravg |
---|
167 | |
---|
168 | AT1 = atan(1/aa) |
---|
169 | AT2 = atan(2/aa) |
---|
170 | // |
---|
171 | // CALCULATIONS ARE PERFORMED TO AVOID LARGE # ERRORS |
---|
172 | // - trick is to propogate the a^(z+7) term through the G1 |
---|
173 | // |
---|
174 | T1 = ZP7*log(aa) - zp1/2*log(aa*aa+4) |
---|
175 | T2 = ZP7*log(aa) - zp3/2*log(aa*aa+4) |
---|
176 | T3 = ZP7*log(aa) - zp2/2*log(aa*aa+4) |
---|
177 | // Print T1,T2,T3 |
---|
178 | RT1 = alog(T1) |
---|
179 | RT2 = alog(T2) |
---|
180 | RT3 = alog(T3) |
---|
181 | V1 = aa^6 - RT1*cos(zp1*at2) |
---|
182 | V2 = ZP1*ZP2*( aa^4 + RT2*cos(zp3*at2) ) |
---|
183 | V3 = -2*ZP1*RT3*SIN(zp2*at2) |
---|
184 | G1 = (V1+V2+V3) |
---|
185 | |
---|
186 | Pq = log(G1) - 6*log(ZP1) + 6*log(Ravg) |
---|
187 | Pq = alog(Pq)*8*PI*PI*delrho*delrho |
---|
188 | |
---|
189 | // |
---|
190 | // beta factor is not used here, but could be for the |
---|
191 | // decoupling approximation |
---|
192 | // |
---|
193 | // G11 = G1 |
---|
194 | // GD = -ZP7*log(aa) |
---|
195 | // G1 = log(G11) + GD |
---|
196 | // |
---|
197 | // T1 = ZP1*at1 |
---|
198 | // T2 = ZP2*at1 |
---|
199 | // G2 = SIN( T1 ) - ZP1/SQRT(aa*aa+1)*COS( T2 ) |
---|
200 | // G22 = G2*G2 |
---|
201 | // BETA = ZP1*log(aa) - ZP1*log(aa*aa+1) - G1 + log(G22) |
---|
202 | // BETA = 2*alog(BETA) |
---|
203 | |
---|
204 | //re-normalize by the average volume |
---|
205 | vpoly = 4*Pi/3*zp3*zp2/zp1/zp1*(ravg)^3 |
---|
206 | Pq /= vpoly |
---|
207 | //scale, convert to cm^-1 |
---|
208 | Pq *= scale * 1e8 |
---|
209 | // add in the background |
---|
210 | Pq += bkg |
---|
211 | |
---|
212 | //return (g1) |
---|
213 | Return (Pq) |
---|
214 | End |
---|
215 | |
---|
216 | //wrapper to calculate the smeared model as an AAO-Struct |
---|
217 | // fills the struct and calls the ususal function with the STRUCT parameter |
---|
218 | // |
---|
219 | // used only for the dependency, not for fitting |
---|
220 | // |
---|
221 | Function fSmearedSchulzSpheres(coefW,yW,xW) |
---|
222 | Wave coefW,yW,xW |
---|
223 | |
---|
224 | String str = getWavesDataFolder(yW,0) |
---|
225 | String DF="root:"+str+":" |
---|
226 | |
---|
227 | WAVE resW = $(DF+str+"_res") |
---|
228 | |
---|
229 | STRUCT ResSmearAAOStruct fs |
---|
230 | WAVE fs.coefW = coefW |
---|
231 | WAVE fs.yW = yW |
---|
232 | WAVE fs.xW = xW |
---|
233 | WAVE fs.resW = resW |
---|
234 | |
---|
235 | Variable err |
---|
236 | err = SmearedSchulzSpheres(fs) |
---|
237 | |
---|
238 | return (0) |
---|
239 | End |
---|
240 | |
---|
241 | // this is all there is to the smeared calculation! |
---|
242 | Function SmearedSchulzSpheres(s) :FitFunc |
---|
243 | Struct ResSmearAAOStruct &s |
---|
244 | |
---|
245 | // the name of your unsmeared model (AAO) is the first argument |
---|
246 | Smear_Model_20(SchulzSpheres,s.coefW,s.xW,s.yW,s.resW) |
---|
247 | |
---|
248 | return(0) |
---|
249 | End |
---|
250 | |
---|
251 | |
---|
252 | //calculates number density given the coefficients of the Schulz distribution |
---|
253 | // the scale factor is the volume fraction |
---|
254 | // then nden = phi/<V> where <V> is calculated using the 3rd moment of the radius |
---|
255 | Macro NumberDensity_Schulz() |
---|
256 | |
---|
257 | Variable nden,zz,zp1,zp2,zp3,zp4,zp5,zp6,zp7,zp8,vpoly,v2poly,rg,sv,i0 |
---|
258 | |
---|
259 | if(Exists("coef_sch")!=1) |
---|
260 | abort "You need to plot the unsmeared model first to create the coefficient table" |
---|
261 | Endif |
---|
262 | |
---|
263 | Print "mean radius (A) = ",coef_sch[1] |
---|
264 | Print "polydispersity (sig/avg) = ",coef_sch[2] |
---|
265 | Print "volume fraction = ",coef_sch[0] |
---|
266 | |
---|
267 | zz = (1/coef_sch[2])^2-1 |
---|
268 | zp1 = zz + 1. |
---|
269 | zp2 = zz + 2. |
---|
270 | zp3 = zz + 3. |
---|
271 | zp4 = zz + 4. |
---|
272 | zp5 = zz + 5. |
---|
273 | zp6 = zz + 6. |
---|
274 | zp7 = zz + 7. |
---|
275 | zp8 = zz + 8. |
---|
276 | // average particle volume |
---|
277 | vpoly = 4*Pi/3*zp3*zp2/zp1/zp1*(coef_sch[1])^3 |
---|
278 | // average particle volume |
---|
279 | v2poly = (4*Pi/3)^2*zp6*zp5*zp4*zp3*zp2/(zp1^5)*(coef_sch[1])^6 |
---|
280 | nden = coef_sch[0]/vpoly //nden in 1/A^3 |
---|
281 | rg = coef_sch[1]*((3*zp8*zp7)/5/zp1/zp1)^0.5 // in A |
---|
282 | sv = 1.0e8*3*coef_sch[0]*zp1/coef_sch[1]/zp3 // in 1/cm |
---|
283 | i0 = 1.0e8*nden*v2poly*(coef_sch[3]-coef_sch[4])^2 // 1/cm/sr |
---|
284 | |
---|
285 | Print "number density (A^-3) = ",nden |
---|
286 | Print "Guinier radius (A) = ",rg |
---|
287 | Print "Interfacial surface area / volume (cm-1) Sv = ",sv |
---|
288 | Print "Forward cross section (cm-1 sr-1) I(0) = ",i0 |
---|
289 | End |
---|
290 | |
---|
291 | // plots the Schulz distribution based on the coefficient values |
---|
292 | // a static calculation, so re-run each time |
---|
293 | // |
---|
294 | Macro PlotSchulzDistribution() |
---|
295 | |
---|
296 | variable pd,avg,zz,maxr |
---|
297 | |
---|
298 | if(Exists("coef_sch")!=1) |
---|
299 | abort "You need to plot the unsmeared model first to create the coefficient table" |
---|
300 | Endif |
---|
301 | pd=coef_sch[2] |
---|
302 | avg = coef_sch[1] |
---|
303 | zz = (1/pd)^2-1 |
---|
304 | |
---|
305 | Make/O/D/N=1000 Schulz_distribution |
---|
306 | maxr = avg*(1+10*pd) |
---|
307 | |
---|
308 | SetScale/I x, 0, maxr, Schulz_distribution |
---|
309 | Schulz_distribution = Schulz_Point(x,avg,zz) |
---|
310 | Display Schulz_distribution |
---|
311 | legend |
---|
312 | End |
---|
313 | |
---|
314 | // don't use the integral technique here since there's an analytic solution |
---|
315 | // available |
---|
316 | ////////////////////// |
---|
317 | //requires that sphere.ipf is included |
---|
318 | // |
---|
319 | // |
---|
320 | //Function SchulzSpheres_Integrated(w,x) |
---|
321 | // Wave w |
---|
322 | // Variable x |
---|
323 | // |
---|
324 | // Variable scale,radius,pd,delrho,bkg,zz,rho,rhos |
---|
325 | // scale = w[0] |
---|
326 | // radius = w[1] |
---|
327 | // pd = w[2] |
---|
328 | // rho = w[3] |
---|
329 | // rhos = w[4] |
---|
330 | // bkg = w[5] |
---|
331 | // |
---|
332 | // delrho=rho-rhos |
---|
333 | // zz = (1/pd)^2-1 |
---|
334 | //// |
---|
335 | //// local variables |
---|
336 | // Variable nord,ii,a,b,va,vb,contr,vcyl,nden,summ,yyy,zi,qq |
---|
337 | // Variable answer,zp1,zp2,zp3,vpoly |
---|
338 | // String weightStr,zStr |
---|
339 | // |
---|
340 | // //select number of gauss points by setting nord=20 or76 points |
---|
341 | //// nord = 20 |
---|
342 | // nord = 76 |
---|
343 | // |
---|
344 | // weightStr = "gauss"+num2str(nord)+"wt" |
---|
345 | // zStr = "gauss"+num2str(nord)+"z" |
---|
346 | // |
---|
347 | // if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
348 | // Make/D/N=(nord) $weightStr,$zStr |
---|
349 | // Wave gauWt = $weightStr |
---|
350 | // Wave gauZ = $zStr // wave references to pass |
---|
351 | // if(nord==20) |
---|
352 | // Make20GaussPoints(gauWt,gauZ) |
---|
353 | // else |
---|
354 | // Make76GaussPoints(gauWt,gauZ) |
---|
355 | // endif |
---|
356 | // else |
---|
357 | // if(exists(weightStr) > 1) |
---|
358 | // Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
359 | // endif |
---|
360 | // Wave gauWt = $weightStr |
---|
361 | // Wave gauZ = $zStr // create the wave references |
---|
362 | // endif |
---|
363 | // |
---|
364 | //// set up the integration |
---|
365 | //// end points and weights |
---|
366 | //// limits are technically 0-inf, but wisely choose non-zero region of distribution |
---|
367 | // Variable range=8 //multiples of the std. dev. from the mean |
---|
368 | // a = radius*(1-range*pd) |
---|
369 | // if (a<0) |
---|
370 | // a=0 //otherwise numerical error when pd >= 0.3, making a<0 |
---|
371 | // endif |
---|
372 | // If(pd>0.3) |
---|
373 | // range = 3.4 + (pd-0.3)*18 //to account for skewed tail |
---|
374 | // Endif |
---|
375 | // b = radius*(1+range*pd) // is this far enough past avg radius? |
---|
376 | // va =a |
---|
377 | // vb =b |
---|
378 | // |
---|
379 | ////temp set scale=1 and bkg=0 for quadrature calc |
---|
380 | // Make/O/D/N=4 sphere_temp |
---|
381 | // sphere_temp[0] = 1 |
---|
382 | // sphere_temp[1] = radius //changed in loop |
---|
383 | // sphere_temp[2] = delrho |
---|
384 | // sphere_temp[3] = 0 |
---|
385 | // |
---|
386 | //// evaluate at Gauss points |
---|
387 | // summ = 0.0 // initialize integral |
---|
388 | // for(ii=0;ii<nord;ii+=1) |
---|
389 | // zi = ( gauZ[ii]*(vb-va) + vb + va )/2.0 |
---|
390 | // sphere_temp[1] = zi |
---|
391 | // yyy = gauWt[ii] * Schulz_Point(zi,radius,zz) * SphereForm(sphere_temp,x) |
---|
392 | // //un-normalize by volume |
---|
393 | // yyy *= 4*pi/3*zi^3 |
---|
394 | // summ += yyy |
---|
395 | // endfor |
---|
396 | //// calculate value of integral to return |
---|
397 | // answer = (vb-va)/2.0*summ |
---|
398 | // |
---|
399 | // //re-normalize by the average volume |
---|
400 | // zp1 = zz + 1. |
---|
401 | // zp2 = zz + 2. |
---|
402 | // zp3 = zz + 3. |
---|
403 | // vpoly = 4*Pi/3*zp3*zp2/zp1/zp1*(radius)^3 |
---|
404 | // answer /= vpoly |
---|
405 | ////scale |
---|
406 | // answer *= scale |
---|
407 | //// add in the background |
---|
408 | // answer += bkg |
---|
409 | // |
---|
410 | // Return (answer) |
---|
411 | //End |
---|
412 | // |
---|
413 | |
---|
414 | //AAO version, uses XOP if available |
---|
415 | // simply calls the original single point calculation with |
---|
416 | // a wave assignment (this will behave nicely if given point ranges) |
---|
417 | Function SchulzSpheresBeta(cw,yw,xw) : FitFunc |
---|
418 | Wave cw,yw,xw |
---|
419 | |
---|
420 | #if exists("SchulzSpheresBetaX") |
---|
421 | yw = SchulzSpheresBetaX(cw,xw) |
---|
422 | #else |
---|
423 | yw = fSchulzSpheresBeta(cw,xw) |
---|
424 | #endif |
---|
425 | return(0) |
---|
426 | End |
---|
427 | |
---|
428 | //use the analytic formula from Kotlarchyk & Chen, JCP 79 (1983) 2461 |
---|
429 | //equations 26 + 23-30 |
---|
430 | // |
---|
431 | // need to calculate in terms of logarithms to avoid numerical errors |
---|
432 | // this returns the Beta factor |
---|
433 | // |
---|
434 | Function fSchulzSpheresBeta(w,x) : FitFunc |
---|
435 | Wave w |
---|
436 | Variable x |
---|
437 | |
---|
438 | Variable scale,ravg,pd,delrho,bkg,zz,rho,rhos,vpoly |
---|
439 | scale = w[0] |
---|
440 | ravg = w[1] |
---|
441 | pd = w[2] |
---|
442 | rho = w[3] |
---|
443 | rhos = w[4] |
---|
444 | bkg = w[5] |
---|
445 | |
---|
446 | delrho=rho-rhos |
---|
447 | zz = (1/pd)^2-1 |
---|
448 | |
---|
449 | Variable zp1,zp2,zp3,zp4,zp5,zp6,zp7 |
---|
450 | Variable aa,b1,b2,b3,at1,at2,rt1,rt2,rt3,t1,t2,t3 |
---|
451 | Variable v1,v2,v3,g1,g11,gd,pq,g2,g22,fBETA |
---|
452 | |
---|
453 | ZP1 = zz + 1 |
---|
454 | ZP2 = zz + 2 |
---|
455 | ZP3 = zz + 3 |
---|
456 | ZP4 = zz + 4 |
---|
457 | ZP5 = zz + 5 |
---|
458 | ZP6 = zz + 6 |
---|
459 | ZP7 = zz + 7 |
---|
460 | |
---|
461 | // //small QR limit - use Guinier approx |
---|
462 | // Variable i_zero,Rg2,zp8 |
---|
463 | // zp8 = zz+8 |
---|
464 | // if(x*ravg < 0.1) |
---|
465 | // i_zero = scale*delrho*delrho*1e8*4*Pi/3*ravg^3 |
---|
466 | // i_zero *= zp6*zp5*zp4/zp1/zp1/zp1 //6th moment / 3rd moment |
---|
467 | // Rg2 = 3*zp8*zp7/5/(zp1^2)*ravg*ravg |
---|
468 | // pq = i_zero*exp(-x*x*Rg2/3) |
---|
469 | // pq += bkg |
---|
470 | // return(pq) |
---|
471 | // endif |
---|
472 | // |
---|
473 | aa = (zz+1)/x/Ravg |
---|
474 | |
---|
475 | AT1 = atan(1/aa) |
---|
476 | AT2 = atan(2/aa) |
---|
477 | // |
---|
478 | // CALCULATIONS ARE PERFORMED TO AVOID LARGE # ERRORS |
---|
479 | // - trick is to propogate the a^(z+7) term through the G1 |
---|
480 | // |
---|
481 | T1 = ZP7*log(aa) - zp1/2*log(aa*aa+4) |
---|
482 | T2 = ZP7*log(aa) - zp3/2*log(aa*aa+4) |
---|
483 | T3 = ZP7*log(aa) - zp2/2*log(aa*aa+4) |
---|
484 | // Print T1,T2,T3 |
---|
485 | RT1 = alog(T1) |
---|
486 | RT2 = alog(T2) |
---|
487 | RT3 = alog(T3) |
---|
488 | V1 = aa^6 - RT1*cos(zp1*at2) |
---|
489 | V2 = ZP1*ZP2*( aa^4 + RT2*cos(zp3*at2) ) |
---|
490 | V3 = -2*ZP1*RT3*SIN(zp2*at2) |
---|
491 | G1 = (V1+V2+V3) |
---|
492 | |
---|
493 | Pq = log(G1) - 6*log(ZP1) + 6*log(Ravg) |
---|
494 | Pq = alog(Pq)*8*PI*PI*delrho*delrho |
---|
495 | |
---|
496 | // |
---|
497 | // beta factor is not used here, but could be for the |
---|
498 | // decoupling approximation |
---|
499 | // |
---|
500 | G11 = G1 |
---|
501 | GD = -ZP7*log(aa) |
---|
502 | G1 = log(G11) + GD |
---|
503 | |
---|
504 | T1 = ZP1*at1 |
---|
505 | T2 = ZP2*at1 |
---|
506 | G2 = SIN( T1 ) - ZP1/SQRT(aa*aa+1)*COS( T2 ) |
---|
507 | G22 = G2*G2 |
---|
508 | fBETA = ZP1*log(aa) - ZP1*log(aa*aa+1) - G1 + log(G22) |
---|
509 | fBETA = 2*alog(fBETA) |
---|
510 | |
---|
511 | return(fBETA) |
---|
512 | |
---|
513 | //re-normalize by the average volume |
---|
514 | // vpoly = 4*Pi/3*zp3*zp2/zp1/zp1*(ravg)^3 |
---|
515 | // Pq /= vpoly |
---|
516 | //scale, convert to cm^-1 |
---|
517 | // Pq *= scale * 1e8 |
---|
518 | // add in the background |
---|
519 | // Pq += bkg |
---|
520 | |
---|
521 | //return (g1) |
---|
522 | // Return (Pq) |
---|
523 | End |
---|