source: sans/Release/trunk/NCNR_User_Procedures/SANS/Analysis/Models/NewModels_2006/StickyHardSphereStruct_v40.ipf @ 381

Last change on this file since 381 was 381, checked in by srkline, 14 years ago

Merging Dev/trunk revision 374+ into Release/trunk for version 6.004

File size: 3.2 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 (A^1) for model: "
28        Prompt qmax "Enter maximum q-value (A^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 (A\\S-1\\M) "
42        Label left "S(q)"
43        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
44       
45        AddModelToStrings("StickyHS_Struct","coef_shsSQ","shsSQ")
46//
47End
48
49//AAO function
50Function StickyHS_Struct(cw,yw,xw) : FitFunc
51        Wave cw,yw,xw
52
53#if exists("StickyHS_StructX")
54        yw = StickyHS_StructX(cw,xw)
55#else
56        yw = fStickyHS_Struct(cw,xw)
57#endif
58        return(0)
59End
60
61Function fStickyHS_Struct(w,x) : FitFunc
62        Wave w
63        Variable x
64//       Input (fitting) variables are:
65        //radius = w[0]
66        //volume fraction = w[1]
67        //epsilon (perurbation param) = w[2]
68        //tau (stickiness) = w[3]
69        Variable rad,phi,eps,tau,eta
70        Variable sig,aa,etam1,qa,qb,qc,radic
71        Variable lam,lam2,test,mu,alpha,beta
72        Variable qv,kk,k2,k3,ds,dc,aq1,aq2,aq3,aq,bq1,bq2,bq3,bq,sq
73        rad = w[0]
74        phi = w[1]
75        eps = w[2]
76        tau = w[3]
77       
78        eta = phi/(1.0-eps)/(1.0-eps)/(1.0-eps)
79       
80        sig = 2.0 * rad
81        aa = sig/(1.0 - eps)
82        etam1 = 1.0 - eta
83//C
84//C  SOLVE QUADRATIC FOR LAMBDA
85//C
86        qa = eta/12.0
87        qb = -1.0*(tau + eta/(etam1))
88        qc = (1.0 + eta/2.0)/(etam1*etam1)
89        radic = qb*qb - 4.0*qa*qc
90        if(radic<0)
91                if(x>0.01 && x<0.015)
92                        Print "Lambda unphysical - both roots imaginary"
93                endif
94                return(-1)
95        endif
96//C   KEEP THE SMALLER ROOT, THE LARGER ONE IS UNPHYSICAL
97        lam = (-1.0*qb-sqrt(radic))/(2.0*qa)
98        lam2 = (-1.0*qb+sqrt(radic))/(2.0*qa)
99        if(lam2<lam)
100                lam = lam2
101        endif
102        test = 1.0 + 2.0*eta
103        mu = lam*eta*etam1
104        if(mu>test)
105                if(x>0.01 && x<0.015)
106                 Print "Lambda unphysical mu>test"
107                endif
108                return(-1)
109        endif
110        alpha = (1.0 + 2.0*eta - mu)/(etam1*etam1)
111        beta = (mu - 3.0*eta)/(2.0*etam1*etam1)
112       
113//C
114//C   CALCULATE THE STRUCTURE FACTOR
115//C
116
117        qv = x
118        kk = qv*aa
119        k2 = kk*kk
120        k3 = kk*k2
121        ds = sin(kk)
122        dc = cos(kk)
123        aq1 = ((ds - kk*dc)*alpha)/k3
124        aq2 = (beta*(1.0-dc))/k2
125        aq3 = (lam*ds)/(12.0*kk)
126        aq = 1.0 + 12.0*eta*(aq1+aq2-aq3)
127//
128        bq1 = alpha*(0.5/kk - ds/k2 + (1.0 - dc)/k3)
129        bq2 = beta*(1.0/kk - ds/k2)
130        bq3 = (lam/12.0)*((1.0 - dc)/kk)
131        bq = 12.0*eta*(bq1+bq2-bq3)
132//
133        sq = 1.0/(aq*aq +bq*bq)
134
135        Return (sq)
136End
137
Note: See TracBrowser for help on using the repository browser.