source: sans/Analysis/branches/ajj_23APR07/IGOR_Package_Files/Put in User Procedures/SANS_Models_v3.00/NewModels_2006/TriaxialEllipsoid.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: 5.8 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
3
4//#include "GaussUtils"
5//#include "PlotUtils"
6
7////////////////////////////////////////////////////
8//
9// calculates the scattering from a triaxial ellipsoid
10// with semi-axes a <= b <= c
11//
12// - the user must make sure that the constraints are not violated
13// otherwise the calculation will not be correct
14//
15// a double integral is used, both using Gaussian quadrature
16// routines that are now included with GaussUtils
17// 20-pt quadrature appears to be enough, 76 pt is available
18// by changing the function calls
19//
20////////////////////////////////////////////////////
21
22//this macro sets up all the necessary parameters and waves that are
23//needed to calculate the model function.
24//
25Proc PlotTriaxialEllipsoid(num,qmin,qmax)
26        Variable num=100, qmin=.001, qmax=.7
27        Prompt num "Enter number of data points for model: "
28        Prompt qmin "Enter minimum q-value (^1) for model: "
29        Prompt qmax "Enter maximum q-value (^1) for model: "
30        //
31        Make/O/D/n=(num) xwave_triax, ywave_triax
32        xwave_triax =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
33        Make/O/D coef_triax = {1,35,100,400,6e-6,0}                     //CH#2
34        make/o/t parameters_triax = {"Scale Factor","Semi-axis A [smallest]()","Semi-axis B ()","Semi-axis C [largest]()","Contrast (^-2)","Incoherent Bgd (cm-1)"} //CH#3
35        Edit parameters_triax, coef_triax
36       
37        Variable/G root:g_triax
38        g_triax := TriaxialEllipsoid(coef_triax, ywave_triax, xwave_triax)
39        Display ywave_triax vs xwave_triax
40        ModifyGraph marker=29, msize=2, mode=4
41        ModifyGraph log=1
42        Label bottom "q (\\S-1\\M) "
43        Label left "I(q) (cm\\S-1\\M)"
44        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
45//
46End
47
48
49// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
50Proc PlotSmearedTriaxialEllipsoid(str)                                                         
51        String str
52        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
53       
54        // if any of the resolution waves are missing => abort
55        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
56                Abort
57        endif
58       
59        SetDataFolder $("root:"+str)
60       
61        // Setup parameter table for model function
62        Make/O/D smear_coef_triax = {1,35,100,400,6e-6,0}               //CH#4
63        make/o/t smear_parameters_triax = {"Scale Factor","A ()","B ()","C ()","Contrast (^-2)","Incoherent Bgd (cm-1)"}
64        Edit smear_parameters_triax,smear_coef_triax                                    //display parameters in a table
65       
66        // output smeared intensity wave, dimensions are identical to experimental QSIG values
67        // make extra copy of experimental q-values for easy plotting
68        Duplicate/O $(str+"_q") smeared_triax,smeared_qvals                             //
69        SetScale d,0,0,"1/cm",smeared_triax                                                     //
70                                       
71        Variable/G gs_triax=0
72        gs_triax := fSmearedTriaxialEllipsoid(smear_coef_triax,smeared_triax,smeared_qvals)     //this wrapper fills the STRUCT
73       
74        Display smeared_triax vs smeared_qvals                                                                  //
75        ModifyGraph log=1,marker=29,msize=2,mode=4
76        Label bottom "q (\\S-1\\M)"
77        Label left "I(q) (cm\\S-1\\M)"
78        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
79       
80        SetDataFolder root:
81End
82
83       
84
85//AAO version, uses XOP if available
86// simply calls the original single point calculation with
87// a wave assignment (this will behave nicely if given point ranges)
88Function TriaxialEllipsoid(cw,yw,xw) : FitFunc
89        Wave cw,yw,xw
90       
91#if exists("TriaxialEllipsoidX")
92        yw = TriaxialEllipsoidX(cw,xw)
93#else
94        yw = fTriaxialEllipsoid(cw,xw)
95#endif
96        return(0)
97End
98
99// calculates the form factor of an ellipsoidal solid
100// with semi-axes of a,b,c
101// - a double integral - choose points wisely
102//
103Function fTriaxialEllipsoid(w,x) : FitFunc
104        Wave w
105        Variable x
106//       Input (fitting) variables are:
107        //[0] scale factor
108        //[1] semi-axis A (A)
109        //[2] semi-axis B (A)
110        //[3] semi-axis C (A)
111        //[4] contrast (A^-2)
112        //[5] incoherent background (cm^-1)
113//      give them nice names
114        Variable scale,aa,bb,cc,contr,bkg,inten,qq,ii,arg,mu
115        scale = w[0]
116        aa = w[1]
117        bb = w[2]
118        cc = w[3]
119        contr = w[4]
120        bkg = w[5]
121       
122        Variable/G root:gDumY=0,root:gDumX=0
123       
124        inten = IntegrateFn20(TaE_Outer,0,1,w,x)
125       
126        inten *= 4*Pi/3*aa*cc*bb                //multiply by volume
127        inten *= 1e8            //convert to cm^-1
128        inten *= contr*contr
129        inten *= scale
130        inten += bkg
131       
132        Return (inten)
133End
134
135// outer integral
136// x is the q-value - remember that "mu" in the notation = B*Q
137Function TaE_Outer(w,x,dum)
138        Wave w
139        Variable x,dum
140       
141        Variable retVal,mu,aa,bb,cc,mudum,arg
142        aa = w[1]
143        bb = w[2]
144        cc = w[3]
145       
146        NVAR dy = root:gDumY
147        NVAR dx = root:gDumX
148        dy = dum
149        retval = IntegrateFn20(TaE_inner,0,1,w,x)
150       
151        return(retVal)
152End
153
154//returns the value of the integrand of the inner integral
155Function TaE_Inner(w,x,dum)
156        Wave w
157        Variable x,dum
158       
159        Variable aa,bb,cc,retVal
160       
161        NVAR dy = root:gDumY
162        NVAR dx = root:gDumX
163        dx = dum
164        aa = w[1]
165        bb = w[2]
166        cc = w[3]
167        retVal = TaE(x,aa,bb,cc,dx,dy)
168       
169        return(retVal)
170End
171
172Function TaE(qq,aa,bb,cc,dx,dy)
173        Variable qq,aa,bb,cc,dx,dy
174       
175        Variable val,arg
176        arg = aa*aa*cos(pi*dx/2)*cos(pi*dx/2)
177        arg += bb*bb*sin(pi*dx/2)*sin(pi*dx/2)*(1-dy*dy)
178        arg += cc*cc*dy*dy
179        arg = qq*sqrt(arg)
180       
181        val = 9*((sin(arg) - arg*cos(arg))/arg^3 )^2
182       
183        return(val)
184end
185
186//wrapper to calculate the smeared model as an AAO-Struct
187// fills the struct and calls the ususal function with the STRUCT parameter
188//
189// used only for the dependency, not for fitting
190//
191Function fSmearedTriaxialEllipsoid(coefW,yW,xW)
192        Wave coefW,yW,xW
193       
194        String str = getWavesDataFolder(yW,0)
195        String DF="root:"+str+":"
196       
197        WAVE resW = $(DF+str+"_res")
198       
199        STRUCT ResSmearAAOStruct fs
200        WAVE fs.coefW = coefW   
201        WAVE fs.yW = yW
202        WAVE fs.xW = xW
203        WAVE fs.resW = resW
204       
205        Variable err
206        err = SmearedTriaxialEllipsoid(fs)
207       
208        return (0)
209End
210               
211// this is all there is to the smeared calculation!
212Function SmearedTriaxialEllipsoid(s) :FitFunc
213        Struct ResSmearAAOStruct &s
214
215//      the name of your unsmeared model (AAO) is the first argument
216        Smear_Model_20(TriaxialEllipsoid,s.coefW,s.xW,s.yW,s.resW)
217
218        return(0)
219End
Note: See TracBrowser for help on using the repository browser.