# source:sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2008/ConvexLens_v40.ipf@449

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

Adding the model of a cylinder with spherical end caps. By taking all of the limiting cases, 5 models are added. Since the geometries are sufficiently different, I thought it best to add as separate models. Corresponding changes have been made to the Wrapper and the modelPicker

File size: 7.0 KB
Line
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
3
4////////////////////////////////////////////////////
5//
6// calculates the scattering of a convex lens.
7//
8// a double integral is used, both using Gaussian quadrature
9// routines that are now included with GaussUtils
10//
12//
13//
14// REFERENCE:
15//              H. Kaya, J. Appl. Cryst. (2004) 37, 223-230.
16//              H. Kaya and N-R deSouza, J. Appl. Cryst. (2004) 37, 508-509. (addenda and errata)
17//
18////////////////////////////////////////////////////
19
20//this macro sets up all the necessary parameters and waves that are
21//needed to calculate the model function.
22//
23Proc PlotConvexLens(num,qmin,qmax)
24        Variable num=100, qmin=.001, qmax=.7
25        Prompt num "Enter number of data points for model: "
26        Prompt qmin "Enter minimum q-value (^1) for model: "
27        Prompt qmax "Enter maximum q-value (^1) for model: "
28        //
29        Make/O/D/n=(num) xwave_ConvLens, ywave_ConvLens
30        xwave_ConvLens =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
31        Make/O/D coef_ConvLens = {1,20,40,1e-6,6.3e-6,0}                        //CH#2
32        make/o/t parameters_ConvLens = {"Scale Factor","cylinder radius rc (A)","end cap radius R >= rc (A)","SLD cylinder (A^-2)","SLD solvent (A^-2)","Incoherent Bgd (cm-1)"}        //CH#3
33        Edit parameters_ConvLens, coef_ConvLens
34
35        Variable/G root:g_ConvLens
36        g_ConvLens := ConvexLens(coef_ConvLens, ywave_ConvLens, xwave_ConvLens)
37        Display ywave_ConvLens vs xwave_ConvLens
38        ModifyGraph marker=29, msize=2, mode=4
39        ModifyGraph log=1
40        Label bottom "q (A\\S-1\\M)"
41        Label left "I(q) (cm\\S-1\\M)"
42        AutoPositionWindow/M=1/R=\$(WinName(0,1)) \$WinName(0,2)
43
45//
46End
47
48
49// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
50Proc PlotSmearedConvexLens(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_ConvLens = {1,20,40,1e-6,6.3e-6,0}          //CH#4
63        make/o/t smear_parameters_ConvLens = {"Scale Factor","cylinder radius rc (A)","end cap radius R >= rc (A)","SLD cylinder (A^-2)","SLD solvent (A^-2)","Incoherent Bgd (cm-1)"}
64        Edit smear_parameters_ConvLens,smear_coef_ConvLens                                      //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_ConvLens,smeared_qvals                          //
69        SetScale d,0,0,"1/cm",smeared_ConvLens                                                  //
70
71        Variable/G gs_ConvLens=0
72        gs_ConvLens := fSmearedConvexLens(smear_coef_ConvLens,smeared_ConvLens,smeared_qvals)   //this wrapper fills the STRUCT
73
74        Display smeared_ConvLens vs smeared_qvals                                                                       //
75        ModifyGraph log=1,marker=29,msize=2,mode=4
76        Label bottom "q (A\\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:
82End
83
84
85
86//AAO version, uses XOP if available
87// simply calls the original single point calculation with
88// a wave assignment (this will behave nicely if given point ranges)
89Function ConvexLens(cw,yw,xw) : FitFunc
90        Wave cw,yw,xw
91
92#if exists("ConvexLensX")
93        yw = ConvexLensX(cw,xw)
94#else
95        yw = fConvexLens(cw,xw)
96#endif
97        return(0)
98End
99
100//
101// - a double integral - choose points wisely - 76 for both...
102//
103Function fConvexLens(w,x) : FitFunc
104        Wave w
105        Variable x
106//       Input (fitting) variables are:
107        //[0] scale factor
108        //[1] cylinder radius (little r)
109        //[2] cylinder length (big L)
110        //[3] end cap radius (big R)
111        //[4] sld cylinder (A^-2)
112        //[5] sld solvent
113        //[6] incoherent background (cm^-1)
114//      give them nice names
115        Variable scale,contr,bkg,inten,sldc,slds
117        scale = w[0]
119//      len = w[2]
121        sldc = w[4]
122        slds = w[5]
123        bkg = w[6]
124
126
127        Make/O/D/N=7 CLens_tmp
128        CLens_tmp[0] = w[0]
129        CLens_tmp[1] = w[1]
130        CLens_tmp[2] = 0.01             //length is some small number, essentially zero
131        CLens_tmp[3] = w[2]
132        CLens_tmp[4] = w[3]
133        CLens_tmp[5] = w[4]
134        CLens_tmp[6] = w[5]
135
136        contr = sldc-slds
137
138        Variable/G root:gDumTheta=0,root:gDumT=0
139
140        inten = IntegrateFn76(ConvLens_Outer,0,pi/2,CLens_tmp,x)
141
142        Variable hh=abs(hDist)          //need positive value for spherical cap volume
144        inten *= 1e8            //convert to cm^-1
145        inten *= contr*contr
146        inten *= scale
147        inten += bkg
148
149        Return (inten)
150End
151
152// outer integral
153// x is the q-value
154Function ConvLens_Outer(w,x,dum)
155        Wave w
156        Variable x,dum
157
158        Variable retVal
159        Variable scale,contr,bkg,inten,sldc,slds
161        scale = w[0]
163        len = w[2]
165        sldc = w[4]
166        slds = w[5]
167        bkg = w[6]
168
170
171        NVAR dTheta = root:gDumTheta
172        NVAR dt = root:gDumT
173        dTheta = dum
175
176        Variable arg1,arg2
177        arg1 = x*len/2*cos(dum)
179
181
182        retVal *= retval*sin(dum)               // = |A(q)|^2*sin(theta)
183
184        return(retVal)
185End
186
187//returns the value of the integrand of the inner integral
188Function ConvLens_Inner(w,x,dum)
189        Wave w
190        Variable x,dum
191
192        Variable retVal
193        Variable scale,contr,bkg,inten,sldc,slds
195        scale = w[0]
197        len = w[2]
199        sldc = w[4]
200        slds = w[5]
201        bkg = w[6]
202
203        NVAR dTheta = root:gDumTheta
204        NVAR dt = root:gDumT
205        dt = dum
206
207        retVal = ConvLens(w,x,dt,dTheta)
208
210
211        return(retVal)
212End
213
214Function ConvLens(w,x,tt,Theta)
215        Wave w
216        Variable x,tt,Theta
217
218        Variable val,arg1,arg2
219        Variable scale,contr,bkg,inten,sldc,slds
221        scale = w[0]
223        len = w[2]
225        sldc = w[4]
226        slds = w[5]
227        bkg = w[6]
228
230
233
234        val = cos(arg1)*(1-tt*tt)*Besselj(1,arg2)/arg2
235
236        return(val)
237end
238
239//wrapper to calculate the smeared model as an AAO-Struct
240// fills the struct and calls the ususal function with the STRUCT parameter
241//
242// used only for the dependency, not for fitting
243//
244Function fSmearedConvexLens(coefW,yW,xW)
245        Wave coefW,yW,xW
246
247        String str = getWavesDataFolder(yW,0)
248        String DF="root:"+str+":"
249
250        WAVE resW = \$(DF+str+"_res")
251
252        STRUCT ResSmearAAOStruct fs
253        WAVE fs.coefW = coefW
254        WAVE fs.yW = yW
255        WAVE fs.xW = xW
256        WAVE fs.resW = resW
257
258        Variable err
259        err = SmearedConvexLens(fs)
260
261        return (0)
262End
263
264// this is all there is to the smeared calculation!
265//
266// 20 points should be fine here. This function is not much different than cylinders, where 20 is sufficient
267Function SmearedConvexLens(s) :FitFunc
268        Struct ResSmearAAOStruct &s
269
270//      the name of your unsmeared model (AAO) is the first argument
271        Smear_Model_20(ConvexLens,s.coefW,s.xW,s.yW,s.resW)
272
273        return(0)
274End
Note: See TracBrowser for help on using the repository browser.