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

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

Bug fixes for some of the newly added (2008) Analysis model functions. All are tested now and should be correct.

Also included in these changes are a few bug fixes in the package loader and Correct

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//
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","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","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#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
116        Variable len,rad,hDist,endRad
117        scale = w[0]
118        rad = w[1]
119//      len = w[2]
120        endRad = w[2]
121        sldc = w[3]
122        slds = w[4]
123        bkg = w[5]
124
125        hDist = -1*sqrt(abs(endRad^2-rad^2))
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
143        inten /= 2*(1/3*pi*(endRad-hh)^2*(2*endRad+hh))         //divide by 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
160        Variable len,rad,hDist,endRad
161        scale = w[0]
162        rad = w[1]
163        len = w[2]
164        endRad = w[3]
165        sldc = w[4]
166        slds = w[5]
167        bkg = w[6]
168       
169        hDist = -1*sqrt(abs(endRad^2-rad^2))
170                       
171        NVAR dTheta = root:gDumTheta
172        NVAR dt = root:gDumT
173        dTheta = dum
174        retval = IntegrateFn76(ConvLens_Inner,-hDist/endRad,1,w,x)
175       
176        Variable arg1,arg2
177        arg1 = x*len/2*cos(dum)
178        arg2 = x*rad*sin(dum)
179       
180        retVal += pi*rad*rad*len*sinc(arg1)*2*Besselj(1, arg2)/arg2
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
194        Variable len,rad,hDist,endRad
195        scale = w[0]
196        rad = w[1]
197        len = w[2]
198        endRad = w[3]
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       
209        retVal *= 4*pi*endRad^3
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
220        Variable len,rad,hDist,endRad
221        scale = w[0]
222        rad = w[1]
223        len = w[2]
224        endRad = w[3]
225        sldc = w[4]
226        slds = w[5]
227        bkg = w[6]
228
229        hDist = -1*sqrt(abs(endRad^2-rad^2))
230               
231        arg1 = x*cos(theta)*(endRad*tt+hDist+len/2)
232        arg2 = x*endRad*sin(theta)*sqrt(1-tt*tt)
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.