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

Last change on this file since 252 was 252, checked in by ajj, 15 years ago

Renaming files with _v40 - see trac ticket #108

File size: 3.8 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
3
4////////////////////////////////////////////////
5// GaussUtils.ipf is not necessary (no smearing calculations done) PlotUtils.ipf is recommended
6////////////////////////////////////////////////
7//
8// this function calculates the interparticle structure factor for spherical particles interacting
9// through a square well potential.
10//
11// 06 NOV 98 SRK
12////////////////////////////////////////////////
13
14Proc PlotSquareWellStruct(num,qmin,qmax)
15        Variable num=128,qmin=0.001,qmax=0.3
16        Prompt num "Enter number of data points for model: "
17        Prompt qmin "Enter minimum q-value (^-1) for model: "
18        Prompt qmax "Enter maximum q-value (^-1) for model: "
19       
20        Make/O/D/n=(num) xwave_sws,ywave_sws
21        xwave_sws = alog(log(qmin) + x*((log(qmax)-log(qmin))/num))             
22        Make/O/D coef_sws = {50,0.04,1.5,1.2}
23        make/o/t parameters_sws = {"Radius (A)","vol fraction","well depth (kT)","well width (diameters)"}
24        Edit parameters_sws,coef_sws
25        Variable/G root:g_sws
26        g_sws := SquareWellStruct(coef_sws,ywave_sws,xwave_sws)
27//      ywave_sws := SquareWellStruct(coef_sws,xwave_sws)
28        Display ywave_sws vs xwave_sws
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("SquareWellStruct","coef_sws","sws")
35End
36
37//AAO version
38Function SquareWellStruct(cw,yw,xw) : FitFunc
39        Wave cw,yw,xw
40       
41#if exists("SquareWellStructX")
42        yw = SquareWellStructX(cw,xw)
43#else
44        yw = fSquareWellStruct(cw,xw)
45#endif
46        return(0)
47End
48///////////////////////////////////////////////////////////
49// smearing a structure factor is not appropriate, so it's not done
50///////////////////////////////////////////////////////////////
51// unsmeared model calculation
52///////////////////////////
53Function fSquareWellStruct(w,x) : FitFunc
54        Wave w
55        Variable x
56
57//     SUBROUTINE SQWELL: CALCULATES THE STRUCTURE FACTOR FOR A
58//                        DISPERSION OF MONODISPERSE HARD SPHERES
59//     IN THE Mean Spherical APPROXIMATION ASSUMING THE SPHERES
60//     INTERACT THROUGH A SQUARE WELL POTENTIAL.
61//** not the best choice of closure ** see note below
62//     REFS:  SHARMA,SHARMA, PHYSICA 89A,(1977),212
63//
64//     
65 
66// NOTE - depths >1.5kT and volume fractions > 0.08 give UNPHYSICAL RESULTS
67// when compared to Monte Carlo simulations
68
69// Input variables are:
70        //[0] radius
71        //[1] volume fraction
72        //[2] well depth e/kT, dimensionless, +ve depths are attractive
73        //[3] well width, multiples of diameter
74       
75        Variable req,phis,edibkb,lambda,struc
76        req = w[0]
77        phis = w[1]
78        edibkb = w[2]
79        lambda = w[3]
80       
81//  COMPUTE CONSTANTS
82//  local variables
83        Variable sigma,eta,eta2,eta3,eta4,etam1,etam14,alpha,beta,gamma
84        Variable qvs,sk,sk2,sk3,sk4,t1,t2,t3,t4,ck
85         
86      SIGMA = req*2.
87      ETA = phis
88      ETA2 = ETA*ETA
89      ETA3 = ETA*ETA2
90      ETA4 = ETA*ETA3       
91      ETAM1 = 1. - ETA
92      ETAM14 = ETAM1*ETAM1*ETAM1*ETAM1
93      ALPHA = (  ( (1. + 2.*ETA)^2 ) + ETA3*( ETA-4.0 )  )/ETAM14
94      BETA = -(ETA/3.0) * ( 18. + 20.*ETA - 12.*ETA2 + ETA4 )/ETAM14
95      GAMMA = 0.5*ETA*( (1. + 2.*ETA)^2 + ETA3*(ETA-4.) )/ETAM14
96//
97//  CALCULATE THE STRUCTURE FACTOR
98//
99// the loop over q-values used to be here     
100//      DO 20 I=1,NPTSM
101        QVS = x
102        SK = x*SIGMA
103        SK2 = SK*SK
104        SK3 = SK*SK2
105        SK4 = SK3*SK
106        T1 = ALPHA * SK3 * ( SIN(SK) - SK * COS(SK) )
107        T2 = BETA * SK2 * ( 2.*SK*SIN(SK) - (SK2-2.)*COS(SK) - 2.0 )
108        T3 =   ( 4.0*SK3 - 24.*SK ) * SIN(SK) 
109        T3 = T3 - ( SK4 - 12.0*SK2 + 24.0 )*COS(SK) + 24.0   
110        T3 = GAMMA*T3
111        T4 = -EDIBKB*SK3*(SIN(LAMBDA*SK) - LAMBDA*SK*COS(LAMBDA*SK)+ SK*COS(SK) - SIN(SK) )
112        CK =  -24.0*ETA*( T1 + T2 + T3 + T4 )/SK3/SK3
113        STRUC  = 1./(1.-CK)
114//   20 CONTINUE
115      Return struc
116End
117
Note: See TracBrowser for help on using the repository browser.