source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/NewModels_2006/BinaryHardSpheres.ipf @ 131

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

Typo in dialog in every model...

File size: 15.8 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 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)
68End
69
70// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
71Proc PlotSmeared_BinaryHS(str)                                                         
72        String str
73        Prompt str,"Pick the data folder containing 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:
103End
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)
110Function 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)
119End
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//
127Function 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)
188End
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)
193Function 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)
202End
203
204Function 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)
244End
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)
249Function 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)
258End
259
260Function 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)
300End
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)
305Function 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)
314End
315
316Function 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)
355End
356
357Function 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)
364End
365
366
367Function 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)
449End
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//
457Function 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)
475End
476
477// this is all there is to the smeared calculation!
478Function 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)
485End
486
487
488
489Macro 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//
518End
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
524Macro 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
540End
541
542//calculates number fractions of each population based on the
543//coef_BinaryHS parameters
544Macro 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       
566End
Note: See TracBrowser for help on using the repository browser.