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

Last change on this file since 633 was 633, checked in by srkline, 13 years ago

Corrected models to explicitly return proper values for I(q=0). There are some models that just can't be fixed, and these typically return NaN. Some, however, are simply numerically unstable at extreme conditions. Beware.

File size: 7.3 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion=6.1
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//
11// 76 point quadrature is necessary for both quadrature calls.
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       
44        AddModelToStrings("ConvexLens","coef_ConvLens","parameters_ConvLens","ConvLens")
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:
81        AddModelToStrings("SmearedConvexLens","smear_coef_ConvLens","smear_parameters_ConvLens","ConvLens")
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//      Variable t1=StopMSTimer(-2)
93
94#if exists("ConvexLensX")
95        MultiThread yw = ConvexLensX(cw,xw)
96#else
97        yw = fConvexLens(cw,xw)
98#endif
99
100//      Print "elapsed time = ",(StopMSTimer(-2) - t1)/1e6
101
102        return(0)
103End
104
105//
106// - a double integral - choose points wisely - 76 for both...
107//
108Function fConvexLens(w,x) : FitFunc
109        Wave w
110        Variable x
111//       Input (fitting) variables are:
112        //[0] scale factor
113        //[1] cylinder radius (little r)
114        //[2] cylinder length (big L)
115        //[3] end cap radius (big R)
116        //[4] sld cylinder (A^-2)
117        //[5] sld solvent
118        //[6] incoherent background (cm^-1)
119//      give them nice names
120        Variable scale,contr,bkg,inten,sldc,slds
121        Variable len,rad,hDist,endRad
122        scale = w[0]
123        rad = w[1]
124//      len = w[2]
125        endRad = w[2]
126        sldc = w[3]
127        slds = w[4]
128        bkg = w[5]
129
130        hDist = -1*sqrt(abs(endRad^2-rad^2))
131
132        Make/O/D/N=7 CLens_tmp
133        CLens_tmp[0] = w[0]
134        CLens_tmp[1] = w[1]
135        CLens_tmp[2] = 0.01             //length is some small number, essentially zero
136        CLens_tmp[3] = w[2]
137        CLens_tmp[4] = w[3]
138        CLens_tmp[5] = w[4]
139        CLens_tmp[6] = w[5]
140               
141        contr = sldc-slds
142       
143        Variable/G root:gDumTheta=0,root:gDumT=0
144       
145        inten = IntegrateFn76(ConvLens_Outer,0,pi/2,CLens_tmp,x)
146       
147        Variable hh=abs(hDist)          //need positive value for spherical cap volume
148        inten /= 2*(1/3*pi*(endRad-hh)^2*(2*endRad+hh))         //divide by volume
149        inten *= 1e8            //convert to cm^-1
150        inten *= contr*contr
151        inten *= scale
152        inten += bkg
153       
154        Return (inten)
155End
156
157// outer integral
158// x is the q-value
159Function ConvLens_Outer(w,x,dum)
160        Wave w
161        Variable x,dum
162       
163        Variable retVal
164        Variable scale,contr,bkg,inten,sldc,slds
165        Variable len,rad,hDist,endRad
166        scale = w[0]
167        rad = w[1]
168        len = w[2]
169        endRad = w[3]
170        sldc = w[4]
171        slds = w[5]
172        bkg = w[6]
173       
174        hDist = -1*sqrt(abs(endRad^2-rad^2))
175                       
176        NVAR dTheta = root:gDumTheta
177        NVAR dt = root:gDumT
178        dTheta = dum
179        retval = IntegrateFn76(ConvLens_Inner,-hDist/endRad,1,w,x)
180       
181        Variable arg1,arg2,be
182        arg1 = x*len/2*cos(dum)
183        arg2 = x*rad*sin(dum)
184       
185        if(arg2 == 0)
186                be = 0.5
187        else
188                be = Besselj(1, arg2)/arg2
189        endif
190       
191        retVal += pi*rad*rad*len*sinc(arg1)*2*be
192       
193        retVal *= retval*sin(dum)               // = |A(q)|^2*sin(theta)
194       
195        return(retVal)
196End
197
198//returns the value of the integrand of the inner integral
199Function ConvLens_Inner(w,x,dum)
200        Wave w
201        Variable x,dum
202       
203        Variable retVal
204        Variable scale,contr,bkg,inten,sldc,slds
205        Variable len,rad,hDist,endRad
206        scale = w[0]
207        rad = w[1]
208        len = w[2]
209        endRad = w[3]
210        sldc = w[4]
211        slds = w[5]
212        bkg = w[6]
213       
214        NVAR dTheta = root:gDumTheta
215        NVAR dt = root:gDumT
216        dt = dum
217       
218        retVal = ConvLens(w,x,dt,dTheta)
219       
220        retVal *= 4*pi*endRad^3
221       
222        return(retVal)
223End
224
225Function ConvLens(w,x,tt,Theta)
226        Wave w
227        Variable x,tt,Theta
228       
229        Variable val,arg1,arg2
230        Variable scale,contr,bkg,inten,sldc,slds
231        Variable len,rad,hDist,endRad,be
232        scale = w[0]
233        rad = w[1]
234        len = w[2]
235        endRad = w[3]
236        sldc = w[4]
237        slds = w[5]
238        bkg = w[6]
239
240        hDist = -1*sqrt(abs(endRad^2-rad^2))
241               
242        arg1 = x*cos(theta)*(endRad*tt+hDist+len/2)
243        arg2 = x*endRad*sin(theta)*sqrt(1-tt*tt)
244       
245        if(arg2 == 0)
246                be = 0.5
247        else
248                be = Besselj(1, arg2)/arg2
249        endif
250       
251        val = cos(arg1)*(1-tt*tt)*be
252       
253        return(val)
254end
255
256//wrapper to calculate the smeared model as an AAO-Struct
257// fills the struct and calls the ususal function with the STRUCT parameter
258//
259// used only for the dependency, not for fitting
260//
261Function fSmearedConvexLens(coefW,yW,xW)
262        Wave coefW,yW,xW
263       
264        String str = getWavesDataFolder(yW,0)
265        String DF="root:"+str+":"
266       
267        WAVE resW = $(DF+str+"_res")
268       
269        STRUCT ResSmearAAOStruct fs
270        WAVE fs.coefW = coefW   
271        WAVE fs.yW = yW
272        WAVE fs.xW = xW
273        WAVE fs.resW = resW
274       
275        Variable err
276        err = SmearedConvexLens(fs)
277       
278        return (0)
279End
280               
281// this is all there is to the smeared calculation!
282//
283// 20 points should be fine here. This function is not much different than cylinders, where 20 is sufficient
284Function SmearedConvexLens(s) :FitFunc
285        Struct ResSmearAAOStruct &s
286
287//      the name of your unsmeared model (AAO) is the first argument
288        Smear_Model_20(ConvexLens,s.coefW,s.xW,s.yW,s.resW)
289
290        return(0)
291End
Note: See TracBrowser for help on using the repository browser.