source: sans/Release/trunk/NCNR_User_Procedures/SANS/Analysis/Models/NewModels_2006/BinaryHardSpheres_v40.ipf @ 345

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

Merging NewModels_2006 back into where it belongs

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