1 | #pragma rtGlobals=1 // Use modern global access method. |
---|
2 | #pragma IgorVersion = 6.0 |
---|
3 | |
---|
4 | |
---|
5 | //////////////////////////////////////////////////// |
---|
6 | // |
---|
7 | // Calculates the scattering from a binary mixture of |
---|
8 | // hard spheres |
---|
9 | // |
---|
10 | // there are some typographical errors in Ashcroft/Langreth's paper |
---|
11 | // Physical Review, v. 156 (1967) 685-692 |
---|
12 | // |
---|
13 | // Errata on Phys. Rev. 166 (1968) 934. |
---|
14 | // |
---|
15 | //(A5) - the entire term should be multiplied by 1/2 |
---|
16 | // |
---|
17 | //final equation for beta12 should be (1+a) rather than (1-a) |
---|
18 | // |
---|
19 | // |
---|
20 | //Definitions are consistent with notation in the paper: |
---|
21 | // |
---|
22 | // phi is total volume fraction |
---|
23 | // nf2 (x) is number density ratio as defined in paper |
---|
24 | // aa = alpha as defined in paper |
---|
25 | // r2 is the radius of the LARGER sphere (angstroms) |
---|
26 | // Sij are the partial structure factor output arrays |
---|
27 | // |
---|
28 | // S. Kline 15 JUL 2004 |
---|
29 | // see: bhs.c and ashcroft.f |
---|
30 | // |
---|
31 | //////////////////////////////////////////////////// |
---|
32 | |
---|
33 | //this macro sets up all the necessary parameters and waves that are |
---|
34 | //needed to calculate the model function. |
---|
35 | // |
---|
36 | // larger sphere radius(angstroms) = guess[0] |
---|
37 | // smaller sphere radius (A) = guess[1] |
---|
38 | // volume fraction of larger spheres = guess[2] |
---|
39 | // volume fraction of small spheres = guess[3] |
---|
40 | // size ratio, alpha(0<a<1) = derived |
---|
41 | // SLD(A-2) of larger particle = guess[4] |
---|
42 | // SLD(A-2) of smaller particle = guess[5] |
---|
43 | // SLD(A-2) of the solvent = guess[6] |
---|
44 | //background = guess[7] |
---|
45 | |
---|
46 | Proc Plot_BinaryHS(num,qmin,qmax) |
---|
47 | Variable num=256, qmin=.001, qmax=.7 |
---|
48 | Prompt num "Enter number of data points for model: " |
---|
49 | Prompt qmin "Enter minimum q-value (-1) for model: " |
---|
50 | Prompt qmax "Enter maximum q-value (-1) for model: " |
---|
51 | // |
---|
52 | Make/O/D/n=(num) xwave_BinaryHS, ywave_BinaryHS |
---|
53 | xwave_BinaryHS = alog(log(qmin) + x*((log(qmax)-log(qmin))/num)) |
---|
54 | Make/O/D coef_BinaryHS = {100,25,0.2,0.1,3.5e-6,0.5e-6,6.36e-6,0.001} //CH#2 |
---|
55 | make/o/t parameters_BinaryHS = {"large radius","small radius","volume fraction large spheres","volume fraction small spheres","large sphere SLD","small sphere SLD","solvent SLD","Incoherent Bgd (cm-1)"} |
---|
56 | Edit parameters_BinaryHS, coef_BinaryHS |
---|
57 | ModifyTable width(parameters_BinaryHS)=160 |
---|
58 | ModifyTable width(coef_BinaryHS)=90 |
---|
59 | |
---|
60 | Variable/G root:g_BinaryHS |
---|
61 | g_BinaryHS := BinaryHS(coef_BinaryHS, ywave_BinaryHS,xwave_BinaryHS) |
---|
62 | Display ywave_BinaryHS vs xwave_BinaryHS |
---|
63 | ModifyGraph marker=29, msize=2, mode=4 |
---|
64 | ModifyGraph log=1 |
---|
65 | Label bottom "q (\\S-1\\M) " |
---|
66 | Label left "I(q) (cm\\S-1\\M)" |
---|
67 | AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2) |
---|
68 | End |
---|
69 | |
---|
70 | // - sets up a dependency to a wrapper, not the actual SmearedModelFunction |
---|
71 | Proc PlotSmeared_BinaryHS(str) |
---|
72 | String str |
---|
73 | Prompt str,"Pick the data folder conatining the resolution you want",popup,getAList(4) |
---|
74 | |
---|
75 | // if any of the resolution waves are missing => abort |
---|
76 | if(ResolutionWavesMissingDF(str)) //updated to NOT use global strings (in GaussUtils) |
---|
77 | Abort |
---|
78 | endif |
---|
79 | |
---|
80 | SetDataFolder $("root:"+str) |
---|
81 | |
---|
82 | // Setup parameter table for model function |
---|
83 | Make/O/D smear_coef_BinaryHS = {100,25,0.2,0.1,3.5e-6,0.5e-6,6.36e-6,0.001} //CH#4 |
---|
84 | make/o/t smear_parameters_BinaryHS = {"large radius","small radius","volume fraction large spheres","volume fraction small spheres","large sphere SLD","small sphere SLD","solvent SLD","Incoherent Bgd (cm-1)"} |
---|
85 | Edit smear_parameters_BinaryHS,smear_coef_BinaryHS //display parameters in a table |
---|
86 | ModifyTable width(smear_parameters_BinaryHS)=160 |
---|
87 | ModifyTable width(smear_coef_BinaryHS)=90 |
---|
88 | // output smeared intensity wave, dimensions are identical to experimental QSIG values |
---|
89 | // make extra copy of experimental q-values for easy plotting |
---|
90 | Duplicate/O $(str+"_q") smeared_BinaryHS,smeared_qvals // |
---|
91 | SetScale d,0,0,"1/cm",smeared_BinaryHS // |
---|
92 | |
---|
93 | Variable/G gs_BinaryHS=0 |
---|
94 | gs_BinaryHS := fBinaryHS_Smeared(smear_coef_BinaryHS,smeared_BinaryHS,smeared_qvals) //this wrapper fills the STRUCT |
---|
95 | |
---|
96 | Display smeared_BinaryHS vs smeared_qvals // |
---|
97 | ModifyGraph log=1,marker=29,msize=2,mode=4 |
---|
98 | Label bottom "q (\\S-1\\M)" |
---|
99 | Label left "I(q) (cm\\S-1\\M)" |
---|
100 | AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2) |
---|
101 | |
---|
102 | SetDataFolder root: |
---|
103 | End |
---|
104 | |
---|
105 | |
---|
106 | |
---|
107 | //AAO version, uses XOP if available |
---|
108 | // simply calls the original single point calculation with |
---|
109 | // a wave assignment (this will behave nicely if given point ranges) |
---|
110 | Function BinaryHS(cw,yw,xw) : FitFunc |
---|
111 | Wave cw,yw,xw |
---|
112 | |
---|
113 | #if exists("BinaryHSX") |
---|
114 | yw = BinaryHSX(cw,xw) |
---|
115 | #else |
---|
116 | yw = fBinaryHS(cw,xw) |
---|
117 | #endif |
---|
118 | return(0) |
---|
119 | End |
---|
120 | |
---|
121 | |
---|
122 | //CH#1 |
---|
123 | // you should write your function to calculate the intensity |
---|
124 | // for a single q-value (that's the input parameter x) |
---|
125 | // based on the wave (array) of parameters that you send it (w) |
---|
126 | // |
---|
127 | Function fBinaryHS(w,x) : FitFunc |
---|
128 | Wave w |
---|
129 | Variable x |
---|
130 | // Input (fitting) variables are: |
---|
131 | // larger sphere radius(angstroms) = guess[0] |
---|
132 | // smaller sphere radius (A) = w[1] |
---|
133 | // number fraction of larger spheres = guess[2] |
---|
134 | // total volume fraction of spheres = guess[3] |
---|
135 | // size ratio, alpha(0<a<1) = derived |
---|
136 | // SLD(A-2) of larger particle = guess[4] |
---|
137 | // SLD(A-2) of smaller particle = guess[5] |
---|
138 | // SLD(A-2) of the solvent = guess[6] |
---|
139 | //background = guess[7] |
---|
140 | |
---|
141 | // give them nice names |
---|
142 | Variable r2,r1,nf2,phi,aa,rho2,rho1,rhos,inten,bgd |
---|
143 | Variable err,psf11,psf12,psf22 |
---|
144 | Variable phi1,phi2,phr,a3 |
---|
145 | |
---|
146 | r2 = w[0] |
---|
147 | r1 = w[1] |
---|
148 | phi2 = w[2] |
---|
149 | phi1 = w[3] |
---|
150 | rho2 = w[4] |
---|
151 | rho1 = w[5] |
---|
152 | rhos = w[6] |
---|
153 | bgd = w[7] |
---|
154 | |
---|
155 | phi = w[2] + w[3] //total volume fraction |
---|
156 | aa = r1/r2 //alpha(0<a<1) |
---|
157 | |
---|
158 | //calculate the number fraction of larger spheres (eqn 2 in reference) |
---|
159 | a3=aa^3 |
---|
160 | phr=phi2/phi |
---|
161 | nf2 = phr*a3/(1-phr+phr*a3) |
---|
162 | // calculate the PSF's here |
---|
163 | |
---|
164 | err = Ashcroft(x,r2,nf2,aa,phi,psf11,psf22,psf12) |
---|
165 | |
---|
166 | // /* do form factor calculations */ |
---|
167 | Variable v1,v2,n1,n2,qr1,qr2,b1,b2 |
---|
168 | v1 = 4.0*PI/3.0*r1*r1*r1 |
---|
169 | v2 = 4.0*PI/3.0*r2*r2*r2 |
---|
170 | // a3 = aa*aa*aa |
---|
171 | // phi1 = phi*(1.0-nf2)*a3/(nf2+(1.0-nf2)*a3) |
---|
172 | // phi2 = phi - phi1 |
---|
173 | n1 = phi1/v1 |
---|
174 | n2 = phi2/v2 |
---|
175 | |
---|
176 | qr1 = r1*x |
---|
177 | qr2 = r2*x |
---|
178 | b1 = r1*r1*r1*(rho1-rhos)*BHSbfunc(qr1) |
---|
179 | b2 = r2*r2*r2*(rho2-rhos)*BHSbfunc(qr2) |
---|
180 | inten = n1*b1*b1*psf11 |
---|
181 | inten += sqrt(n1*n2)*2.0*b1*b2*psf12 |
---|
182 | inten += n2*b2*b2*psf22 |
---|
183 | ///* convert I(1/A) to (1/cm) */ |
---|
184 | inten *= 1.0e8 |
---|
185 | |
---|
186 | inten += bgd |
---|
187 | Return (inten) |
---|
188 | End |
---|
189 | |
---|
190 | //AAO version, uses XOP if available |
---|
191 | // simply calls the original single point calculation with |
---|
192 | // a wave assignment (this will behave nicely if given point ranges) |
---|
193 | Function BinaryHS_PSF11(cw,yw,xw) : FitFunc |
---|
194 | Wave cw,yw,xw |
---|
195 | |
---|
196 | #if exists("BinaryHS_PSF11X") |
---|
197 | yw = BinaryHS_PSF11X(cw,xw) |
---|
198 | #else |
---|
199 | yw = fBinaryHS_PSF11(cw,xw) |
---|
200 | #endif |
---|
201 | return(0) |
---|
202 | End |
---|
203 | |
---|
204 | Function fBinaryHS_PSF11(w,x) : FitFunc |
---|
205 | Wave w |
---|
206 | Variable x |
---|
207 | // Input (fitting) variables are: |
---|
208 | // larger sphere radius(angstroms) = guess[0] |
---|
209 | // smaller sphere radius (A) = w[1] |
---|
210 | // number fraction of larger spheres = guess[2] |
---|
211 | // total volume fraction of spheres = guess[3] |
---|
212 | // size ratio, alpha(0<a<1) = derived |
---|
213 | // SLD(A-2) of larger particle = guess[4] |
---|
214 | // SLD(A-2) of smaller particle = guess[5] |
---|
215 | // SLD(A-2) of the solvent = guess[6] |
---|
216 | //background = guess[7] |
---|
217 | |
---|
218 | // give them nice names |
---|
219 | Variable r2,r1,nf2,phi,aa,rho2,rho1,rhos,inten,bgd |
---|
220 | Variable err,psf11,psf12,psf22 |
---|
221 | Variable phi1,phi2,a3,phr |
---|
222 | |
---|
223 | r2 = w[0] |
---|
224 | r1 = w[1] |
---|
225 | phi2 = w[2] |
---|
226 | phi1 = w[3] |
---|
227 | rho2 = w[4] |
---|
228 | rho1 = w[5] |
---|
229 | rhos = w[6] |
---|
230 | bgd = w[7] |
---|
231 | |
---|
232 | phi = w[2] + w[3] //total volume fraction |
---|
233 | aa = r1/r2 //alpha(0<a<1) |
---|
234 | |
---|
235 | //calculate the number fraction of larger spheres (eqn 2 in reference) |
---|
236 | a3=aa^3 |
---|
237 | phr=phi2/phi |
---|
238 | nf2 = phr*a3/(1-phr+phr*a3) |
---|
239 | |
---|
240 | // calculate the PSF's here |
---|
241 | |
---|
242 | err = Ashcroft(x,r2,nf2,aa,phi,psf11,psf22,psf12) |
---|
243 | return(psf11) |
---|
244 | End |
---|
245 | |
---|
246 | //AAO version, uses XOP if available |
---|
247 | // simply calls the original single point calculation with |
---|
248 | // a wave assignment (this will behave nicely if given point ranges) |
---|
249 | Function BinaryHS_PSF12(cw,yw,xw) : FitFunc |
---|
250 | Wave cw,yw,xw |
---|
251 | |
---|
252 | #if exists("BinaryHS_PSF12X") |
---|
253 | yw = BinaryHS_PSF12X(cw,xw) |
---|
254 | #else |
---|
255 | yw = fBinaryHS_PSF12(cw,xw) |
---|
256 | #endif |
---|
257 | return(0) |
---|
258 | End |
---|
259 | |
---|
260 | Function fBinaryHS_PSF12(w,x) : FitFunc |
---|
261 | Wave w |
---|
262 | Variable x |
---|
263 | // Input (fitting) variables are: |
---|
264 | // larger sphere radius(angstroms) = guess[0] |
---|
265 | // smaller sphere radius (A) = w[1] |
---|
266 | // number fraction of larger spheres = guess[2] |
---|
267 | // total volume fraction of spheres = guess[3] |
---|
268 | // size ratio, alpha(0<a<1) = derived |
---|
269 | // SLD(A-2) of larger particle = guess[4] |
---|
270 | // SLD(A-2) of smaller particle = guess[5] |
---|
271 | // SLD(A-2) of the solvent = guess[6] |
---|
272 | //background = guess[7] |
---|
273 | |
---|
274 | // give them nice names |
---|
275 | Variable r2,r1,nf2,phi,aa,rho2,rho1,rhos,inten,bgd |
---|
276 | Variable err,psf11,psf12,psf22 |
---|
277 | Variable phi1,phi2,a3,phr |
---|
278 | |
---|
279 | r2 = w[0] |
---|
280 | r1 = w[1] |
---|
281 | phi2 = w[2] |
---|
282 | phi1 = w[3] |
---|
283 | rho2 = w[4] |
---|
284 | rho1 = w[5] |
---|
285 | rhos = w[6] |
---|
286 | bgd = w[7] |
---|
287 | |
---|
288 | phi = w[2] + w[3] //total volume fraction |
---|
289 | aa = r1/r2 //alpha(0<a<1) |
---|
290 | |
---|
291 | //calculate the number fraction of larger spheres (eqn 2 in reference) |
---|
292 | a3=aa^3 |
---|
293 | phr=phi2/phi |
---|
294 | nf2 = phr*a3/(1-phr+phr*a3) |
---|
295 | |
---|
296 | // calculate the PSF's here |
---|
297 | |
---|
298 | err = Ashcroft(x,r2,nf2,aa,phi,psf11,psf22,psf12) |
---|
299 | return(psf12) |
---|
300 | End |
---|
301 | |
---|
302 | //AAO version, uses XOP if available |
---|
303 | // simply calls the original single point calculation with |
---|
304 | // a wave assignment (this will behave nicely if given point ranges) |
---|
305 | Function BinaryHS_PSF22(cw,yw,xw) : FitFunc |
---|
306 | Wave cw,yw,xw |
---|
307 | |
---|
308 | #if exists("BinaryHS_PSF22X") |
---|
309 | yw = BinaryHS_PSF22X(cw,xw) |
---|
310 | #else |
---|
311 | yw = fBinaryHS_PSF22(cw,xw) |
---|
312 | #endif |
---|
313 | return(0) |
---|
314 | End |
---|
315 | |
---|
316 | Function fBinaryHS_PSF22(w,x) : FitFunc |
---|
317 | Wave w |
---|
318 | Variable x |
---|
319 | // Input (fitting) variables are: |
---|
320 | // larger sphere radius(angstroms) = guess[0] |
---|
321 | // smaller sphere radius (A) = w[1] |
---|
322 | // number fraction of larger spheres = guess[2] |
---|
323 | // total volume fraction of spheres = guess[3] |
---|
324 | // size ratio, alpha(0<a<1) = derived |
---|
325 | // SLD(A-2) of larger particle = guess[4] |
---|
326 | // SLD(A-2) of smaller particle = guess[5] |
---|
327 | // SLD(A-2) of the solvent = guess[6] |
---|
328 | //background = guess[7] |
---|
329 | |
---|
330 | // give them nice names |
---|
331 | Variable r2,r1,nf2,phi,aa,rho2,rho1,rhos,inten,bgd |
---|
332 | Variable err,psf11,psf12,psf22 |
---|
333 | Variable phi1,phi2,phr,a3 |
---|
334 | |
---|
335 | r2 = w[0] |
---|
336 | r1 = w[1] |
---|
337 | phi2 = w[2] |
---|
338 | phi1 = w[3] |
---|
339 | rho2 = w[4] |
---|
340 | rho1 = w[5] |
---|
341 | rhos = w[6] |
---|
342 | bgd = w[7] |
---|
343 | |
---|
344 | phi = w[2] + w[3] //total volume fraction |
---|
345 | aa = r1/r2 //alpha(0<a<1) |
---|
346 | |
---|
347 | //calculate the number fraction of larger spheres (eqn 2 in reference) |
---|
348 | a3=aa^3 |
---|
349 | phr=phi2/phi |
---|
350 | nf2 = phr*a3/(1-phr+phr*a3) |
---|
351 | // calculate the PSF's here |
---|
352 | |
---|
353 | err = Ashcroft(x,r2,nf2,aa,phi,psf11,psf22,psf12) |
---|
354 | return(psf22) |
---|
355 | End |
---|
356 | |
---|
357 | Function BHSbfunc(qr) |
---|
358 | Variable qr |
---|
359 | |
---|
360 | Variable ans |
---|
361 | |
---|
362 | ans = 4.0*pi*(sin(qr)-qr*cos(qr))/qr/qr/qr |
---|
363 | return(ans) |
---|
364 | End |
---|
365 | |
---|
366 | |
---|
367 | Function Ashcroft(qval,r2,nf2,aa,phi,s11,s22,s12) |
---|
368 | Variable qval,r2,nf2,aa,phi,&s11,&s22,&s12 |
---|
369 | |
---|
370 | // CALCULATE CONSTANT TERMS |
---|
371 | Variable s1,s2,v,a3,v1,v2,g11,g12,g22,wmv,wmv3,wmv4 |
---|
372 | Variable a1,a2i,a2,b1,b2,b12,gm1,gm12 |
---|
373 | Variable err,yy,ay,ay2,ay3,t1,t2,t3,f11,y2,y3,tt1,tt2,tt3 |
---|
374 | Variable c11,c22,c12,f12,f22,ttt1,ttt2,ttt3,ttt4,yl,y13 |
---|
375 | Variable t21,t22,t23,t31,t32,t33,t41,t42,yl3,wma3,y1 |
---|
376 | |
---|
377 | s2 = 2.0*r2 |
---|
378 | s1 = aa*s2 |
---|
379 | v = phi |
---|
380 | a3 = aa*aa*aa |
---|
381 | V1=((1.-nf2)*A3/(nf2+(1.-nf2)*A3))*V |
---|
382 | V2=(nf2/(nf2+(1.-nf2)*A3))*V |
---|
383 | G11=((1.+.5*V)+1.5*V2*(aa-1.))/(1.-V)/(1.-V) |
---|
384 | G22=((1.+.5*V)+1.5*V1*(1./aa-1.))/(1.-V)/(1.-v) |
---|
385 | G12=((1.+.5*V)+1.5*(1.-aa)*(V1-V2)/(1.+aa))/(1.-V)/(1.-v) |
---|
386 | wmv = 1/(1.-v) |
---|
387 | wmv3 = wmv*wmv*wmv |
---|
388 | wmv4 = wmv*wmv3 |
---|
389 | A1=3.*wmv4*((V1+A3*V2)*(1.+V+V*v)-3.*V1*V2*(1.-aa)*(1.-aa)*(1.+V1+aa*(1.+V2))) + ((V1+A3*V2)*(1.+2.*V)+(1.+V+V*v)-3.*V1*V2*(1.-aa)*(1.-aa)-3.*V2*(1.-aa)*(1.-aa)*(1.+V1+aa*(1.+V2)))*wmv3 |
---|
390 | A2I=((V1+A3*V2)*(1.+V+V*v)-3.*V1*V2*(1.-aa)*(1.-aa)*(1.+V1+aa*(1.+V2)))*3*wmv4 + ((V1+A3*V2)*(1.+2.*V)+A3*(1.+V+V*v)-3.*V1*V2*(1.-aa)*(1.-aa)*aa-3.*V1*(1.-aa)*(1.-aa)*(1.+V1+aa*(1.+V2)))*wmv3 |
---|
391 | A2=A2I/a3 |
---|
392 | B1=-6.*(V1*G11*g11+.25*V2*(1.+aa)*(1.+aa)*aa*G12*g12) |
---|
393 | B2=-6.*(V2*G22*g22+.25*V1/A3*(1.+aa)*(1.+aa)*G12*g12) |
---|
394 | B12=-3.*aa*(1.+aa)*(V1*G11/aa/aa+V2*G22)*G12 |
---|
395 | GM1=(V1*A1+A3*V2*A2)*.5 |
---|
396 | GM12=2.*GM1*(1.-aa)/aa |
---|
397 | //C |
---|
398 | //C CALCULATE THE DIRECT CORRELATION FUNCTIONS AND PRINT RESULTS |
---|
399 | //C |
---|
400 | // DO 20 J=1,npts |
---|
401 | |
---|
402 | yy=qval*s2 |
---|
403 | //c calculate direct correlation functions |
---|
404 | //c ----c11 |
---|
405 | AY=aa*yy |
---|
406 | ay2 = ay*ay |
---|
407 | ay3 = ay*ay*ay |
---|
408 | T1=A1*(SIN(AY)-AY*COS(AY)) |
---|
409 | T2=B1*(2.*AY*sin(AY)-(AY2-2.)*cos(AY)-2.)/AY |
---|
410 | T3=GM1*((4.*AY*ay2-24.*AY)*sin(AY)-(AY2*ay2-12.*AY2+24.)*cos(AY)+24.)/AY3 |
---|
411 | F11=24.*V1*(T1+T2+T3)/AY3 |
---|
412 | |
---|
413 | //c ------c22 |
---|
414 | y2=yy*yy |
---|
415 | y3=yy*y2 |
---|
416 | TT1=A2*(sin(yy)-yy*cos(yy)) |
---|
417 | TT2=B2*(2.*yy*sin(yy)-(Y2-2.)*cos(yy)-2.)/yy |
---|
418 | TT3=GM1*((4.*Y3-24.*yy)*sin(yy)-(Y2*y2-12.*Y2+24.)*cos(yy)+24.)/ay3 |
---|
419 | F22=24.*V2*(TT1+TT2+TT3)/Y3 |
---|
420 | |
---|
421 | //c -----c12 |
---|
422 | YL=.5*yy*(1.-aa) |
---|
423 | yl3=yl*yl*yl |
---|
424 | wma3 = (1.-aa)*(1.-aa)*(1.-aa) |
---|
425 | Y1=aa*yy |
---|
426 | y13 = y1*y1*y1 |
---|
427 | TTT1=3.*wma3*V*sqrt(nf2)*sqrt(1.-nf2)*A1*(sin(YL)-YL*cos(YL))/((nf2+(1.-nf2)*A3)*YL3) |
---|
428 | T21=B12*(2.*Y1*cos(Y1)+(Y1^2-2.)*sin(Y1)) |
---|
429 | T22=GM12*((3.*Y1*y1-6.)*cos(Y1)+(Y1^3-6.*Y1)*sin(Y1)+6.)/Y1 |
---|
430 | T23=GM1*((4.*Y13-24.*Y1)*cos(Y1)+(Y13*y1-12.*Y1*y1+24.)*sin(Y1))/(Y1*y1) |
---|
431 | T31=B12*(2.*Y1*sin(Y1)-(Y1^2-2.)*cos(Y1)-2.) |
---|
432 | T32=GM12*((3.*Y1^2-6.)*sin(Y1)-(Y1^3-6.*Y1)*cos(Y1))/Y1 |
---|
433 | T33=GM1*((4.*Y13-24.*Y1)*sin(Y1)-(Y13*y1-12.*Y1*y1+24.)*cos(Y1)+24.)/(y1*y1) |
---|
434 | T41=cos(YL)*((sin(Y1)-Y1*cos(Y1))/(Y1*y1) + (1.-aa)/(2.*aa)*(1.-cos(Y1))/Y1) |
---|
435 | T42=sin(YL)*((cos(Y1)+Y1*sin(Y1)-1.)/(Y1*y1) + (1.-aa)/(2.*aa)*sin(Y1)/Y1) |
---|
436 | TTT2=sin(YL)*(T21+T22+T23)/(y13*y1) |
---|
437 | TTT3=cos(YL)*(T31+T32+T33)/(y13*y1) |
---|
438 | TTT4=A1*(T41+T42)/Y1 |
---|
439 | F12=TTT1+24.*V*sqrt(nf2)*sqrt(1.-nf2)*A3*(TTT2+TTT3+TTT4)/(nf2+(1.-nf2)*A3) |
---|
440 | |
---|
441 | C11=F11 |
---|
442 | C22=F22 |
---|
443 | C12=F12 |
---|
444 | S11=1./(1.+C11-(C12)*c12/(1.+C22)) |
---|
445 | S22=1./(1.+C22-(C12)*c12/(1.+C11)) |
---|
446 | S12=-C12/((1.+C11)*(1.+C22)-(C12)*(c12)) |
---|
447 | |
---|
448 | return(err) |
---|
449 | End |
---|
450 | |
---|
451 | |
---|
452 | //wrapper to calculate the smeared model as an AAO-Struct |
---|
453 | // fills the struct and calls the ususal function with the STRUCT parameter |
---|
454 | // |
---|
455 | // used only for the dependency, not for fitting |
---|
456 | // |
---|
457 | Function fBinaryHS_Smeared(coefW,yW,xW) |
---|
458 | Wave coefW,yW,xW |
---|
459 | |
---|
460 | String str = getWavesDataFolder(yW,0) |
---|
461 | String DF="root:"+str+":" |
---|
462 | |
---|
463 | WAVE resW = $(DF+str+"_res") |
---|
464 | |
---|
465 | STRUCT ResSmearAAOStruct fs |
---|
466 | WAVE fs.coefW = coefW |
---|
467 | WAVE fs.yW = yW |
---|
468 | WAVE fs.xW = xW |
---|
469 | WAVE fs.resW = resW |
---|
470 | |
---|
471 | Variable err |
---|
472 | err = BinaryHS_Smeared(fs) |
---|
473 | |
---|
474 | return (0) |
---|
475 | End |
---|
476 | |
---|
477 | // this is all there is to the smeared calculation! |
---|
478 | Function BinaryHS_Smeared(s) :FitFunc |
---|
479 | Struct ResSmearAAOStruct &s |
---|
480 | |
---|
481 | // the name of your unsmeared model (AAO) is the first argument |
---|
482 | s.yW = Smear_Model_20(BinaryHS,s.coefW,s.xW,s.resW) |
---|
483 | |
---|
484 | return(0) |
---|
485 | End |
---|
486 | |
---|
487 | |
---|
488 | |
---|
489 | Macro Plot_BinaryHS_PSF() |
---|
490 | |
---|
491 | if(WaveExists(coef_BinaryHS)==0) |
---|
492 | abort "You need to plot the model first to create the coefficient table" |
---|
493 | Endif |
---|
494 | |
---|
495 | Make/O/D/n=(numpnts(xwave_BinaryHS)) psf11_BinaryHS,psf12_BinaryHS,psf22_BinaryHS,QD2_BinaryHS |
---|
496 | |
---|
497 | psf11_BinaryHS := BinaryHS_psf11(coef_BinaryHS, xwave_BinaryHS) |
---|
498 | psf12_BinaryHS := BinaryHS_psf12(coef_BinaryHS, xwave_BinaryHS) |
---|
499 | psf22_BinaryHS := BinaryHS_psf22(coef_BinaryHS, xwave_BinaryHS) |
---|
500 | QD2_BinaryHS := xwave_BinaryHS*coef_BinaryHS[0]*2 |
---|
501 | // Display psf11_BinaryHS vs xwave_BinaryHS |
---|
502 | // AppendtoGraph psf12_BinaryHS vs xwave_BinaryHS |
---|
503 | // AppendtoGraph psf22_BinaryHS vs xwave_BinaryHS |
---|
504 | |
---|
505 | Display psf11_BinaryHS vs QD2_BinaryHS |
---|
506 | AppendtoGraph psf12_BinaryHS vs QD2_BinaryHS |
---|
507 | AppendtoGraph psf22_BinaryHS vs QD2_BinaryHS |
---|
508 | |
---|
509 | ModifyGraph marker=19, msize=2, mode=4 |
---|
510 | ModifyGraph lsize=2,rgb(psf12_BinaryHS)=(2,39321,1) |
---|
511 | ModifyGraph rgb(psf22_BinaryHS)=(0,0,65535) |
---|
512 | ModifyGraph log=0,grid=1,mirror=2 |
---|
513 | SetAxis bottom 0,30 |
---|
514 | Label bottom "q*LargeDiameter" |
---|
515 | Label left "Sij(q)" |
---|
516 | Legend |
---|
517 | // |
---|
518 | End |
---|
519 | |
---|
520 | |
---|
521 | //useful for finding the parameters that duplicate the |
---|
522 | //figures in the original reference (uses the same notation) |
---|
523 | //automatically changes the coefficient wave |
---|
524 | Macro Duplicate_AL_Parameters(eta,xx,alpha,Rlarge) |
---|
525 | Variable eta=0.45,xx=0.4,alpha=0.7,Rlarge=100 |
---|
526 | |
---|
527 | Variable r1,phi1,phi2,a3 |
---|
528 | r1 = alpha*Rlarge |
---|
529 | a3 = alpha*alpha*alpha |
---|
530 | phi1 = eta*(1.0-xx)*a3/(xx+(1.0-xx)*a3) //eqn [2] |
---|
531 | phi2 = eta - phi1 |
---|
532 | |
---|
533 | Print "phi (larger) = ",phi2 |
---|
534 | Print "phi (smaller) = ",phi1 |
---|
535 | Print "Radius (smaller) = ",r1 |
---|
536 | |
---|
537 | coef_BinaryHS[2] = phi2 |
---|
538 | coef_BinaryHS[3] = phi1 |
---|
539 | coef_BinaryHS[1] = r1 |
---|
540 | End |
---|
541 | |
---|
542 | //calculates number fractions of each population based on the |
---|
543 | //coef_BinaryHS parameters |
---|
544 | Macro Calculate_BHS_Parameters() |
---|
545 | |
---|
546 | if(exists("coef_BinaryHS") != 1) |
---|
547 | Abort "You must plot the BHS model first to create the coefficient wave" |
---|
548 | endif |
---|
549 | Variable r1,r2,phi1,phi2,aa //same notation as paper - r2 is LARGER |
---|
550 | Variable a3,xx,phi |
---|
551 | r1 = coef_BinaryHS[1] |
---|
552 | r2 = coef_BinaryHS[0] |
---|
553 | phi1 = coef_BinaryHS[3] |
---|
554 | phi2 = coef_BinaryHS[2] |
---|
555 | |
---|
556 | phi = phi1+phi2 |
---|
557 | aa = r1/r2 |
---|
558 | a3 = aa^3 |
---|
559 | |
---|
560 | xx = phi2/phi*a3 |
---|
561 | xx /= (1-(phi2/phi)+(phi2/phi)*a3) |
---|
562 | |
---|
563 | Print "Number fraction (larger) = ",xx |
---|
564 | Print "Number fraction (smaller) = ",1-xx |
---|
565 | |
---|
566 | End |
---|