source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/NewModels_2006/StickyHardSphereStruct.ipf @ 153

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

Changed Plot* and PlotSmeared?* naming schemes to be all consistent prefixes for the actual function name, so that the macros can be constructed from the function name, or vice versa.

also some tweaks to the wrapper to make sure that plot and append really work

File size: 3.1 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
3
4////////////////////////////////////////////////////
5//
6// description
7//
8//
9// reference -
10// - errors in paper
11//
12// warning about unphysical regions
13//
14// + how to properly handle if result is unphysical????
15//
16//
17// S. Kline 15 JUL 2004
18//
19////////////////////////////////////////////////////
20
21//this macro sets up all the necessary parameters and waves that are
22//needed to calculate the model function.
23//
24Proc PlotStickyHS_Struct(num,qmin,qmax)
25        Variable num=256, qmin=.001, qmax=.5
26        Prompt num "Enter number of data points for model: "
27        Prompt qmin "Enter minimum q-value (^1) for model: "
28        Prompt qmax "Enter maximum q-value (^1) for model: "
29//
30        Make/O/D/n=(num) xwave_shsSQ, ywave_shsSQ
31        xwave_shsSQ =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
32        Make/O/D coef_shsSQ = {50,0.1,0.05,0.2}                 //CH#2
33        make/o/t parameters_shsSQ = {"Radius","volume fraction","perturbation parameter (0.1)","stickiness, tau"}       //CH#3
34        Edit parameters_shsSQ, coef_shsSQ
35        Variable/G root:g_shsSQ
36        g_shsSQ  := StickyHS_Struct(coef_shsSQ, ywave_shsSQ, xwave_shsSQ)
37//      ywave_shsSQ  := StickyHS_Struct(coef_shsSQ, xwave_shsSQ)
38        Display ywave_shsSQ vs xwave_shsSQ
39        ModifyGraph marker=29, msize=2, mode=4
40        ModifyGraph log=0
41        Label bottom "q (\\S-1\\M) "
42        Label left "S(q)"
43        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
44//
45End
46
47//AAO function
48Function StickyHS_Struct(cw,yw,xw) : FitFunc
49        Wave cw,yw,xw
50
51#if exists("StickyHS_StructX")
52        yw = StickyHS_StructX(cw,xw)
53#else
54        yw = fStickyHS_Struct(cw,xw)
55#endif
56        return(0)
57End
58
59Function fStickyHS_Struct(w,x) : FitFunc
60        Wave w
61        Variable x
62//       Input (fitting) variables are:
63        //radius = w[0]
64        //volume fraction = w[1]
65        //epsilon (perurbation param) = w[2]
66        //tau (stickiness) = w[3]
67        Variable rad,phi,eps,tau,eta
68        Variable sig,aa,etam1,qa,qb,qc,radic
69        Variable lam,lam2,test,mu,alpha,beta
70        Variable qv,kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,sq
71        rad = w[0]
72        phi = w[1]
73        eps = w[2]
74        tau = w[3]
75       
76        eta = phi/(1.0-eps)/(1.0-eps)/(1.0-eps)
77       
78        sig = 2.0 * rad
79        aa = sig/(1.0 - eps)
80        etam1 = 1.0 - eta
81//C
82//C  SOLVE QUADRATIC FOR LAMBDA
83//C
84        qa = eta/12.0
85        qb = -1.0*(tau + eta/(etam1))
86        qc = (1.0 + eta/2.0)/(etam1*etam1)
87        radic = qb*qb - 4.0*qa*qc
88        if(radic<0)
89                if(x>0.01 && x<0.015)
90                        Print "Lambda unphysical - both roots imaginary"
91                endif
92                return(-1)
93        endif
94//C   KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL
95        lam = (-1.0*qb-sqrt(radic))/(2.0*qa)
96        lam2 = (-1.0*qb+sqrt(radic))/(2.0*qa)
97        if(lam2<lam)
98                lam = lam2
99        endif
100        test = 1.0 + 2.0*eta
101        mu = lam*eta*etam1
102        if(mu>test)
103                if(x>0.01 && x<0.015)
104                 Print "Lambda unphysical mu>test"
105                endif
106                return(-1)
107        endif
108        alpha = (1.0 + 2.0*eta - mu)/(etam1*etam1)
109        beta = (mu - 3.0*eta)/(2.0*etam1*etam1)
110       
111//C
112//C   CALCULATE THE STRUCTURE FACTOR
113//C
114
115        qv = x
116        kk = qv*aa
117        k2 = kk*kk
118        k3 = kk*k2
119        ds = sin(kk)
120        dc = cos(kk)
121        aq1 = ((ds - kk*dc)*alpha)/k3
122        aq2 = (beta*(1.0-dc))/k2
123        aq3 = (lam*ds)/(12.0*kk)
124        aq = 1.0 + 12.0*eta*(aq1+aq2-aq3)
125//
126        bq1 = alpha*(0.5/kk - ds/k2 + (1.0 - dc)/k3)
127        bq2 = beta*(1.0/kk - ds/k2)
128        bq3 = (lam/12.0)*((1.0 - dc)/kk)
129        bq = 12.0*eta*(bq1+bq2-bq3)
130//
131        sq = 1.0/(aq*aq +bq*bq)
132
133        Return (sq)
134End
135
Note: See TracBrowser for help on using the repository browser.