source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/HardSphereStruct.ipf @ 145

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

LOTS of changes to the analysis ipf files:

-- see sphere.ipf for the simplest example of the changes --

  • #pragma Igor 6
  • #if directive to look for XOP
  • AAO unsmeared functions
  • STRUCT functions for smearing (also AAO)
  • new wrappers for dependencies to struct functions

(2006 models have NOT been completed yet, only the old models)

  • loading data files into data folders (PlotUtils?) + some streamlining of the loaders
  • Smear_Model_N is now AAO + some streamlining of the quadrature code

-- SHS and SW structure factor XOPs are crashing (need DP wave, I may have old XOP)
-- this breaks fitting of the smeared models until wrappers can be devised
-- all packages will be broken due to the new data folder structure
-- multiple instances of functions will now cause problems (MSA)
-- RPA model is a problem with its odd functional form (extra wave)

-- lots of other carnage to follow as the bugs and typos are shaken out

24 JUL 2007 SRK

File size: 3.3 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)
33End
34
35//AAO version
36Function 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)
45End
46
47///////////////////////////////////////////////////////////
48// no smeared version is calculated - it is simply not appropriate
49
50///////////////////////////////////////////////////////////////
51// unsmeared model calculation
52///////////////////////////
53Function 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
108End
109// End of HardSphereStruct
110
Note: See TracBrowser for help on using the repository browser.