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

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

LOTS of changes to the analysis ipf files:

-- see sphere.ipf for the simplest example of the changes --

  • #pragma Igor 6
  • #if directive to look for XOP
  • AAO unsmeared functions
  • STRUCT functions for smearing (also AAO)
  • new wrappers for dependencies to struct functions

(2006 models have NOT been completed yet, only the old models)

  • loading data files into data folders (PlotUtils?) + some streamlining of the loaders
  • Smear_Model_N is now AAO + some streamlining of the quadrature code

-- SHS and SW structure factor XOPs are crashing (need DP wave, I may have old XOP)
-- this breaks fitting of the smeared models until wrappers can be devised
-- all packages will be broken due to the new data folder structure
-- multiple instances of functions will now cause problems (MSA)
-- RPA model is a problem with its odd functional form (extra wave)

-- lots of other carnage to follow as the bugs and typos are shaken out

24 JUL 2007 SRK

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 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//
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.