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

Revision 570, 7.2 KB checked in by srkline, 4 years ago (diff)

Change (1):
In preparation for release, updated pragma IgorVersion?=6.1 in all procedures

Change (2):
As a side benefit of requiring 6.1, we can use the MultiThread? keyword to thread any model function we like. The speed benefit is only noticeable on functions that require at least one integration and at least 100 points (resolution smearing is NOT threaded, too many threadSafe issues, too little benefit). I have chosen to use the MultiThread? only on the XOP assignment. In the Igor code there are too many functions that are not explicitly declared threadsafe, making for a mess.

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//
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//      Variable t1=StopMSTimer(-2)
93
94#if exists("ConvexLensX")
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
122        scale = w[0]
124//      len = w[2]
126        sldc = w[3]
127        slds = w[4]
128        bkg = w[5]
129
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
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
166        scale = w[0]
168        len = w[2]
170        sldc = w[4]
171        slds = w[5]
172        bkg = w[6]
173
175
176        NVAR dTheta = root:gDumTheta
177        NVAR dt = root:gDumT
178        dTheta = dum
180
181        Variable arg1,arg2
182        arg1 = x*len/2*cos(dum)
184
186
187        retVal *= retval*sin(dum)               // = |A(q)|^2*sin(theta)
188
189        return(retVal)
190End
191
192//returns the value of the integrand of the inner integral
193Function ConvLens_Inner(w,x,dum)
194        Wave w
195        Variable x,dum
196
197        Variable retVal
198        Variable scale,contr,bkg,inten,sldc,slds
200        scale = w[0]
202        len = w[2]
204        sldc = w[4]
205        slds = w[5]
206        bkg = w[6]
207
208        NVAR dTheta = root:gDumTheta
209        NVAR dt = root:gDumT
210        dt = dum
211
212        retVal = ConvLens(w,x,dt,dTheta)
213
215
216        return(retVal)
217End
218
219Function ConvLens(w,x,tt,Theta)
220        Wave w
221        Variable x,tt,Theta
222
223        Variable val,arg1,arg2
224        Variable scale,contr,bkg,inten,sldc,slds
226        scale = w[0]
228        len = w[2]
230        sldc = w[4]
231        slds = w[5]
232        bkg = w[6]
233
235
238
239        val = cos(arg1)*(1-tt*tt)*Besselj(1,arg2)/arg2
240
241        return(val)
242end
243
244//wrapper to calculate the smeared model as an AAO-Struct
245// fills the struct and calls the ususal function with the STRUCT parameter
246//
247// used only for the dependency, not for fitting
248//
249Function fSmearedConvexLens(coefW,yW,xW)
250        Wave coefW,yW,xW
251
252        String str = getWavesDataFolder(yW,0)
253        String DF="root:"+str+":"
254
255        WAVE resW = \$(DF+str+"_res")
256
257        STRUCT ResSmearAAOStruct fs
258        WAVE fs.coefW = coefW
259        WAVE fs.yW = yW
260        WAVE fs.xW = xW
261        WAVE fs.resW = resW
262
263        Variable err
264        err = SmearedConvexLens(fs)
265
266        return (0)
267End
268
269// this is all there is to the smeared calculation!
270//
271// 20 points should be fine here. This function is not much different than cylinders, where 20 is sufficient
272Function SmearedConvexLens(s) :FitFunc
273        Struct ResSmearAAOStruct &s
274
275//      the name of your unsmeared model (AAO) is the first argument
276        Smear_Model_20(ConvexLens,s.coefW,s.xW,s.yW,s.resW)
277
278        return(0)
279End
Note: See TracBrowser for help on using the repository browser.