# source:sans/Dev/trunk/NCNR_User_Procedures/SANS/Analysis/Models/NewModels_2006/TriaxialEllipsoid_v40.ipf@379

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

Removed Angstrom symbol (Mac code) and replaced it with simply a capital "A" so strange symbols won't appear anymore on Win.

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