source: sans/Release/trunk/NCNR_User_Procedures/SANS/Analysis/Models_v3.00/NewModels_2006/BinaryHardSpheres.ipf @ 338

Last change on this file since 338 was 338, checked in by srkline, 15 years ago

merging structural changes from Dev/trunk into Release/trunk

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