source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v4.00/HardSphereStruct_v40.ipf @ 278

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

changed variable name in hsStruct calculation to avoid conflict with new Igor function (beta)

File size: 3.4 KB
Line 
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
13Proc 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       
34        AddModelToStrings("HardSphereStruct","coef_hss","hss")
35End
36
37//AAO version
38Function HardSphereStruct(cw,yw,xw) : FitFunc
39        Wave cw,yw,xw
40
41#if exists("HardSphereStructX")
42        yw = HardSphereStructX(cw,xw)
43#else
44        yw = fHardSphereStruct(cw,xw)
45#endif
46        return(0)
47End
48
49///////////////////////////////////////////////////////////
50// no smeared version is calculated - it is simply not appropriate
51
52///////////////////////////////////////////////////////////////
53// unsmeared model calculation
54///////////////////////////
55Function fHardSphereStruct(w,x) : FitFunc
56        Wave w
57        Variable x
58                       
59//     SUBROUTINE HSSTRCT: CALCULATES THE STRUCTURE FACTOR FOR A
60//                         DISPERSION OF MONODISPERSE HARD SPHERES
61//                         IN THE PERCUS-YEVICK APPROXIMATION
62//
63//     REFS:  PERCUS,YEVICK PHYS. REV. 110 1 (1958)
64//            THIELE J. CHEM PHYS. 39 474 (1968)
65//            WERTHEIM  PHYS. REV. LETT. 47 1462 (1981)
66//
67// Input variables are:
68        //[0] radius
69        //[1] volume fraction
70        //Variable timer=StartMSTimer
71       
72        Variable r,phi,struc
73        r = w[0]
74        phi = w[1]
75
76// Local variables
77        Variable denom,dnum,alpha,fbeta,gamm,q,a,asq,ath,afor,rca,rsa
78        Variable calp,cbeta,cgam,prefac,c,vstruc
79//  COMPUTE CONSTANTS
80//
81      DENOM = (1.0-PHI)^4
82      DNUM = (1.0 + 2.0*PHI)^2
83      ALPHA = DNUM/DENOM
84      fBETA = -6.0*PHI*((1.0 + PHI/2.0)^2)/DENOM
85      GAMM = 0.50*PHI*DNUM/DENOM
86//
87//  CALCULATE THE STRUCTURE FACTOR
88//     
89// loop over q-values used to be here
90//      DO 10 I=1,NPTSM
91        Q = x           // q-value for the calculation is passed in as variable x
92        A = 2.0*Q*R
93//        IF(A.LT.1.0D-10)  A = 1.0D-10
94        ASQ = A*A
95        ATH = ASQ*A
96        AFOR = ATH*A
97        RCA = COS(A)
98        RSA = SIN(A)
99        CALP = ALPHA*(RSA/ASQ - RCA/A)
100        CBETA = fBETA*(2.0*RSA/ASQ - (ASQ - 2.0)*RCA/ATH - 2.0/ATH)
101        CGAM = GAMM*(-RCA/A + (4.0/A)*((3.0*ASQ - 6.0)*RCA/AFOR + (ASQ - 6.0)*RSA/ATH + 6.0/AFOR))
102        PREFAC = -24.0*PHI/A
103        C = PREFAC*(CALP + CBETA + CGAM)
104        VSTRUC = 1.0/(1.0-C)
105        STRUC = VSTRUC
106//   10 CONTINUE
107        //Variable elapse=StopMSTimer(timer)
108      //Print "HS struct eval time (s) = ",elapse
109      RETURN Struc
110End
111// End of HardSphereStruct
112
Note: See TracBrowser for help on using the repository browser.