source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/BinaryHardSpheres_v40.ipf @ 510

Last change on this file since 510 was 510, checked in by srkline, 14 years ago

Simple change in all of the model function files to include the name of the parameter wave in the Keyword=list that is generated when a model is plotted. This is becoming an issue where the proper parameter wave can't be deduced from just the suffix, then there is nothing to put in the table.

I should have added this when I initially wrote the wrapper...

File size: 16.2 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 (A-1) for model: "
50        Prompt qmax "Enter maximum q-value (A-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 (A\\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","parameters_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 (A\\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","smear_parameters_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.