source: sans/Analysis/trunk/Put in User Procedures/SANS_Models_v3.00/SquareWellStruct.ipf @ 56

Last change on this file since 56 was 42, checked in by srkline, 16 years ago

initial checkin of Analysis v.3.00 files

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