source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2008/CappedCylinder_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.2 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 "capped cylinder" with "flat" spherical end caps
7// where the radius of the end cap is larger than the radius of the cylinder.
8// The center of the spherical end caps is within the length of the cylinder.
9//
10// a double integral is used, both using Gaussian quadrature
11// routines that are now included with GaussUtils
12//
13// 76 point quadrature is necessary for both quadrature calls.
14//
15//
16// REFERENCE:
17//              H. Kaya, J. Appl. Cryst. (2004) 37, 223-230.
18//              H. Kaya and N-R deSouza, J. Appl. Cryst. (2004) 37, 508-509. (addenda and errata)
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 PlotCappedCylinder(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 (^1) for model: "
29        Prompt qmax "Enter maximum q-value (^1) for model: "
30        //
31        Make/O/D/n=(num) xwave_CapCyl, ywave_CapCyl
32        xwave_CapCyl =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
33        Make/O/D coef_CapCyl = {1,20,400,40,1e-6,6.3e-6,0}                      //CH#2
34        make/o/t parameters_CapCyl = {"Scale Factor","cylinder radius rc (A)","cylinder length (A)","end cap radius R >= rc (A)","SLD cylinder (A^-2)","SLD solvent (A^-2)","Incoherent Bgd (cm-1)"}    //CH#3
35        Edit parameters_CapCyl, coef_CapCyl
36       
37        Variable/G root:g_CapCyl
38        g_CapCyl := CappedCylinder(coef_CapCyl, ywave_CapCyl, xwave_CapCyl)
39        Display ywave_CapCyl vs xwave_CapCyl
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       
46        AddModelToStrings("CappedCylinder","coef_CapCyl","parameters_CapCyl","CapCyl")
47//
48End
49
50
51// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
52Proc PlotSmearedCappedCylinder(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_CapCyl = {1,20,400,40,1e-6,6.3e-6,0}                //CH#4
65        make/o/t smear_parameters_CapCyl = {"Scale Factor","cylinder radius rc (A)","cylinder length (A)","end cap radius R >= rc (A)","SLD cylinder (A^-2)","SLD solvent (A^-2)","Incoherent Bgd (cm-1)"}
66        Edit smear_parameters_CapCyl,smear_coef_CapCyl                                  //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_CapCyl,smeared_qvals                            //
71        SetScale d,0,0,"1/cm",smeared_CapCyl                                                    //
72                                       
73        Variable/G gs_CapCyl=0
74        gs_CapCyl := fSmearedCappedCylinder(smear_coef_CapCyl,smeared_CapCyl,smeared_qvals)     //this wrapper fills the STRUCT
75       
76        Display smeared_CapCyl 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:
83        AddModelToStrings("SmearedCappedCylinder","smear_coef_CapCyl","smear_parameters_CapCyl","CapCyl")
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 CappedCylinder(cw,yw,xw) : FitFunc
92        Wave cw,yw,xw
93       
94#if exists("CappedCylinderX")
95        MultiThread yw = CappedCylinderX(cw,xw)
96#else
97        yw = fCappedCylinder(cw,xw)
98#endif
99        return(0)
100End
101
102//
103// - a double integral - choose points wisely - 76 for both...
104//
105Function fCappedCylinder(w,x) : FitFunc
106        Wave w
107        Variable x
108//       Input (fitting) variables are:
109        //[0] scale factor
110        //[1] cylinder radius (little r)
111        //[2] cylinder length (big L)
112        //[3] end cap radius (big R)
113        //[4] sld cylinder (A^-2)
114        //[5] sld solvent
115        //[6] incoherent background (cm^-1)
116//      give them nice names
117        Variable scale,contr,bkg,inten,sldc,slds
118        Variable len,rad,hDist,endRad
119        scale = w[0]
120        rad = w[1]
121        len = w[2]
122        endRad = w[3]
123        sldc = w[4]
124        slds = w[5]
125        bkg = w[6]
126
127        hDist = -1*sqrt(abs(endRad^2-rad^2))
128               
129        contr = sldc-slds
130       
131        Variable/G root:gDumTheta=0,root:gDumT=0
132       
133        inten = IntegrateFn76(CapCyl_Outer,0,pi/2,w,x)
134       
135        Variable hh=abs(hdist)          //need a positive h for the volume of the spherical section
136        inten /= pi*rad*rad*len + 2*(1/3*pi*(endRad-hh)^2*(2*endRad+hh))                //divide by volume
137        inten *= 1e8            //convert to cm^-1
138        inten *= contr*contr
139        inten *= scale
140        inten += bkg
141       
142        Return (inten)
143End
144
145// outer integral
146// x is the q-value
147Function CapCyl_Outer(w,x,dum)
148        Wave w
149        Variable x,dum
150       
151        Variable retVal
152        Variable scale,contr,bkg,inten,sldc,slds
153        Variable len,rad,hDist,endRad
154        scale = w[0]
155        rad = w[1]
156        len = w[2]
157        endRad = w[3]
158        sldc = w[4]
159        slds = w[5]
160        bkg = w[6]
161
162        hDist = -1*sqrt(abs(endRad^2-rad^2))
163                       
164        NVAR dTheta = root:gDumTheta
165        NVAR dt = root:gDumT
166        dTheta = dum
167        retval = IntegrateFn76(CapCyl_Inner,-hDist/endRad,1,w,x)
168       
169        Variable arg1,arg2,be
170        arg1 = x*len/2*cos(dum)
171        arg2 = x*rad*sin(dum)
172       
173        if(arg2 == 0)
174                be = 0.5
175        else
176                be = Besselj(1, arg2)/arg2
177        endif
178        retVal += pi*rad*rad*len*sinc(arg1)*2*be
179       
180        retVal *= retval*sin(dum)               // = |A(q)|^2*sin(theta)
181       
182        return(retVal)
183End
184
185//returns the value of the integrand of the inner integral
186Function CapCyl_Inner(w,x,dum)
187        Wave w
188        Variable x,dum
189       
190        Variable retVal
191        Variable scale,contr,bkg,inten,sldc,slds
192        Variable len,rad,hDist,endRad
193        scale = w[0]
194        rad = w[1]
195        len = w[2]
196        endRad = w[3]
197        sldc = w[4]
198        slds = w[5]
199        bkg = w[6]
200       
201        NVAR dTheta = root:gDumTheta
202        NVAR dt = root:gDumT
203        dt = dum
204       
205        retVal = CapCyl(w,x,dt,dTheta)
206       
207        retVal *= 4*pi*endRad^3
208       
209        return(retVal)
210End
211
212Function CapCyl(w,x,tt,Theta)
213        Wave w
214        Variable x,tt,Theta
215       
216        Variable val,arg1,arg2
217        Variable scale,contr,bkg,inten,sldc,slds
218        Variable len,rad,hDist,endRad,be
219        scale = w[0]
220        rad = w[1]
221        len = w[2]
222        endRad = w[3]
223        sldc = w[4]
224        slds = w[5]
225        bkg = w[6]
226
227        hDist = -1*sqrt(abs(endRad^2-rad^2))
228               
229        arg1 = x*cos(theta)*(endRad*tt+hDist+len/2)
230        arg2 = x*endRad*sin(theta)*sqrt(1-tt*tt)
231       
232        if(arg2 == 0)
233                be = 0.5
234        else
235                be = Besselj(1, arg2)/arg2
236        endif
237       
238        val = cos(arg1)*(1-tt*tt)*be
239       
240        return(val)
241end
242
243//wrapper to calculate the smeared model as an AAO-Struct
244// fills the struct and calls the ususal function with the STRUCT parameter
245//
246// used only for the dependency, not for fitting
247//
248Function fSmearedCappedCylinder(coefW,yW,xW)
249        Wave coefW,yW,xW
250       
251        String str = getWavesDataFolder(yW,0)
252        String DF="root:"+str+":"
253       
254        WAVE resW = $(DF+str+"_res")
255       
256        STRUCT ResSmearAAOStruct fs
257        WAVE fs.coefW = coefW   
258        WAVE fs.yW = yW
259        WAVE fs.xW = xW
260        WAVE fs.resW = resW
261       
262        Variable err
263        err = SmearedCappedCylinder(fs)
264       
265        return (0)
266End
267               
268// this is all there is to the smeared calculation!
269//
270// 20 points should be fine here. This function is not much different than cylinders, where 20 is sufficient
271Function SmearedCappedCylinder(s) :FitFunc
272        Struct ResSmearAAOStruct &s
273
274//      the name of your unsmeared model (AAO) is the first argument
275        Smear_Model_20(CappedCylinder,s.coefW,s.xW,s.yW,s.resW)
276
277        return(0)
278End
Note: See TracBrowser for help on using the repository browser.