source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2006/StickyHardSphereStruct_v40.ipf @ 510

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

Simple change in all of the model function files to include the name of the parameter wave in the Keyword=list that is generated when a model is plotted. This is becoming an issue where the proper parameter wave can't be deduced from just the suffix, then there is nothing to put in the table.

I should have added this when I initially wrote the wrapper...

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","parameters_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
Note: See TracBrowser for help on using the repository browser.