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

Last change on this file since 128 was 128, checked in by srkline, 16 years ago

Second round of changes to ipf files to (hopefully) complete the changes to bring all of the model functions to be:

  • this set of changes is for the "2006" set of models
  • AAO
  • new dependencies
  • smearing using structures
  • #if conditional compile for XOPS

-this needs to be FULLY tested for correctness and bugs
-packages are certainly still broken (data folders)

SRK 25 JUL 07
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 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:
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.