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 | // |
---|
24 | Proc Plot_StickyHS_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 | // |
---|
45 | End |
---|
46 | |
---|
47 | //AAO function |
---|
48 | Function 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) |
---|
57 | End |
---|
58 | |
---|
59 | Function 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) |
---|
134 | End |
---|
135 | |
---|