1 | #pragma rtGlobals=1 // Use modern global access method. |
---|
2 | #pragma IgorVersion = 6.0 |
---|
3 | |
---|
4 | //////////////////////////////////////////////// |
---|
5 | // GaussUtils.ipf is not required (no smearing calculation done) PlotUtils.ipf should be included |
---|
6 | //////////////////////////////////////////////// |
---|
7 | // this function is for the hard sphere structure factor |
---|
8 | // |
---|
9 | // |
---|
10 | // 06 NOV 98 SRK |
---|
11 | //////////////////////////////////////////////// |
---|
12 | |
---|
13 | Proc PlotHardSphereStruct(num,qmin,qmax) |
---|
14 | Variable num=128,qmin=0.001,qmax=0.3 |
---|
15 | Prompt num "Enter number of data points for model: " |
---|
16 | Prompt qmin "Enter minimum q-value (^-1) for model: " |
---|
17 | Prompt qmax "Enter maximum q-value (^-1) for model: " |
---|
18 | |
---|
19 | Make/O/D/n=(num) xwave_hss,ywave_hss |
---|
20 | //xwave_hss = (x+1)*((qmax-qmin)/num) |
---|
21 | xwave_hss = alog(log(qmin) + x*((log(qmax)-log(qmin))/num)) |
---|
22 | Make/O/D coef_hss = {50,0.2} |
---|
23 | make/o/t parameters_hss = {"Radius (A)","vol fraction"} |
---|
24 | Edit parameters_hss,coef_hss |
---|
25 | Variable/G root:g_hss |
---|
26 | g_hss := HardSphereStruct(coef_hss,ywave_hss,xwave_hss) |
---|
27 | // ywave_hss := HardSphereStruct(coef_hss,xwave_hss) |
---|
28 | Display ywave_hss vs xwave_hss |
---|
29 | ModifyGraph marker=29,msize=2,mode=4 |
---|
30 | Label bottom "q (\\S-1\\M)" |
---|
31 | Label left "Structure Factor" |
---|
32 | AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2) |
---|
33 | End |
---|
34 | |
---|
35 | //AAO version |
---|
36 | Function HardSphereStruct(cw,yw,xw) : FitFunc |
---|
37 | Wave cw,yw,xw |
---|
38 | |
---|
39 | #if exists("HardSphereStructX") |
---|
40 | yw = HardSphereStructX(cw,xw) |
---|
41 | #else |
---|
42 | yw = fHardSphereStruct(cw,xw) |
---|
43 | #endif |
---|
44 | return(0) |
---|
45 | End |
---|
46 | |
---|
47 | /////////////////////////////////////////////////////////// |
---|
48 | // no smeared version is calculated - it is simply not appropriate |
---|
49 | |
---|
50 | /////////////////////////////////////////////////////////////// |
---|
51 | // unsmeared model calculation |
---|
52 | /////////////////////////// |
---|
53 | Function fHardSphereStruct(w,x) : FitFunc |
---|
54 | Wave w |
---|
55 | Variable x |
---|
56 | |
---|
57 | // SUBROUTINE HSSTRCT: CALCULATES THE STRUCTURE FACTOR FOR A |
---|
58 | // DISPERSION OF MONODISPERSE HARD SPHERES |
---|
59 | // IN THE PERCUS-YEVICK APPROXIMATION |
---|
60 | // |
---|
61 | // REFS: PERCUS,YEVICK PHYS. REV. 110 1 (1958) |
---|
62 | // THIELE J. CHEM PHYS. 39 474 (1968) |
---|
63 | // WERTHEIM PHYS. REV. LETT. 47 1462 (1981) |
---|
64 | // |
---|
65 | // Input variables are: |
---|
66 | //[0] radius |
---|
67 | //[1] volume fraction |
---|
68 | //Variable timer=StartMSTimer |
---|
69 | |
---|
70 | Variable r,phi,struc |
---|
71 | r = w[0] |
---|
72 | phi = w[1] |
---|
73 | |
---|
74 | // Local variables |
---|
75 | Variable denom,dnum,alpha,beta,gamm,q,a,asq,ath,afor,rca,rsa |
---|
76 | Variable calp,cbeta,cgam,prefac,c,vstruc |
---|
77 | // COMPUTE CONSTANTS |
---|
78 | // |
---|
79 | DENOM = (1.0-PHI)^4 |
---|
80 | DNUM = (1.0 + 2.0*PHI)^2 |
---|
81 | ALPHA = DNUM/DENOM |
---|
82 | BETA = -6.0*PHI*((1.0 + PHI/2.0)^2)/DENOM |
---|
83 | GAMM = 0.50*PHI*DNUM/DENOM |
---|
84 | // |
---|
85 | // CALCULATE THE STRUCTURE FACTOR |
---|
86 | // |
---|
87 | // loop over q-values used to be here |
---|
88 | // DO 10 I=1,NPTSM |
---|
89 | Q = x // q-value for the calculation is passed in as variable x |
---|
90 | A = 2.0*Q*R |
---|
91 | // IF(A.LT.1.0D-10) A = 1.0D-10 |
---|
92 | ASQ = A*A |
---|
93 | ATH = ASQ*A |
---|
94 | AFOR = ATH*A |
---|
95 | RCA = COS(A) |
---|
96 | RSA = SIN(A) |
---|
97 | CALP = ALPHA*(RSA/ASQ - RCA/A) |
---|
98 | CBETA = BETA*(2.0*RSA/ASQ - (ASQ - 2.0)*RCA/ATH - 2.0/ATH) |
---|
99 | CGAM = GAMM*(-RCA/A + (4.0/A)*((3.0*ASQ - 6.0)*RCA/AFOR + (ASQ - 6.0)*RSA/ATH + 6.0/AFOR)) |
---|
100 | PREFAC = -24.0*PHI/A |
---|
101 | C = PREFAC*(CALP + CBETA + CGAM) |
---|
102 | VSTRUC = 1.0/(1.0-C) |
---|
103 | STRUC = VSTRUC |
---|
104 | // 10 CONTINUE |
---|
105 | //Variable elapse=StopMSTimer(timer) |
---|
106 | //Print "HS struct eval time (s) = ",elapse |
---|
107 | RETURN Struc |
---|
108 | End |
---|
109 | // End of HardSphereStruct |
---|
110 | |
---|