source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/Models/NewModels_2008/Spherocylinder_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 spherocylinder, that is a cylinder with spherical end caps
7// where the radius of the end caps is the same as the radius of the cylinder
8//
9// a double integral is used, both using Gaussian quadrature
10// routines that are now included with GaussUtils
11//
12// 76 point quadrature is necessary for both quadrature calls.
13//
14//
15// REFERENCE:
16//              H. Kaya, J. Appl. Cryst. (2004) 37, 223-230.
17//              H. Kaya and N-R deSouza, J. Appl. Cryst. (2004) 37, 508-509. (addenda and errata)
18//
19////////////////////////////////////////////////////
20
21//this macro sets up all the necessary parameters and waves that are
22//needed to calculate the model function.
23//
24Proc PlotSpherocylinder(num,qmin,qmax)
25        Variable num=100, qmin=.001, qmax=.7
26        Prompt num "Enter number of data points for model: "
27        Prompt qmin "Enter minimum q-value (^1) for model: "
28        Prompt qmax "Enter maximum q-value (^1) for model: "
29        //
30        Make/O/D/n=(num) xwave_SphCyl, ywave_SphCyl
31        xwave_SphCyl =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
32        Make/O/D coef_SphCyl = {1,20,400,1e-6,6.3e-6,0}                 //CH#2
33        make/o/t parameters_SphCyl = {"Scale Factor","cylinder radius rc (A)","cylinder length (A)","SLD cylinder (A^-2)","SLD solvent (A^-2)","Incoherent Bgd (cm-1)"} //CH#3
34        Edit parameters_SphCyl, coef_SphCyl
35       
36        Variable/G root:g_SphCyl
37        g_SphCyl := Spherocylinder(coef_SphCyl, ywave_SphCyl, xwave_SphCyl)
38        Display ywave_SphCyl vs xwave_SphCyl
39        ModifyGraph marker=29, msize=2, mode=4
40        ModifyGraph log=1
41        Label bottom "q (A\\S-1\\M)"
42        Label left "I(q) (cm\\S-1\\M)"
43        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
44       
45        AddModelToStrings("Spherocylinder","coef_SphCyl","parameters_SphCyl","SphCyl")
46//
47End
48
49
50// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
51Proc PlotSmearedSpherocylinder(str)                                                             
52        String str
53        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
54       
55        // if any of the resolution waves are missing => abort
56        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
57                Abort
58        endif
59       
60        SetDataFolder $("root:"+str)
61       
62        // Setup parameter table for model function
63        Make/O/D smear_coef_SphCyl = {1,20,400,1e-6,6.3e-6,0}           //CH#4
64        make/o/t smear_parameters_SphCyl = {"Scale Factor","cylinder radius rc (A)","cylinder length (A)","SLD cylinder (A^-2)","SLD solvent (A^-2)","Incoherent Bgd (cm-1)"}
65        Edit smear_parameters_SphCyl,smear_coef_SphCyl                                  //display parameters in a table
66       
67        // output smeared intensity wave, dimensions are identical to experimental QSIG values
68        // make extra copy of experimental q-values for easy plotting
69        Duplicate/O $(str+"_q") smeared_SphCyl,smeared_qvals                            //
70        SetScale d,0,0,"1/cm",smeared_SphCyl                                                    //
71                                       
72        Variable/G gs_SphCyl=0
73        gs_SphCyl := fSmearedSpherocylinder(smear_coef_SphCyl,smeared_SphCyl,smeared_qvals)     //this wrapper fills the STRUCT
74       
75        Display smeared_SphCyl vs smeared_qvals                                                                 //
76        ModifyGraph log=1,marker=29,msize=2,mode=4
77        Label bottom "q (A\\S-1\\M)"
78        Label left "I(q) (cm\\S-1\\M)"
79        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
80       
81        SetDataFolder root:
82        AddModelToStrings("SmearedSpherocylinder","smear_coef_SphCyl","smear_parameters_SphCyl","SphCyl")
83End
84
85       
86
87//AAO version, uses XOP if available
88// simply calls the original single point calculation with
89// a wave assignment (this will behave nicely if given point ranges)
90Function Spherocylinder(cw,yw,xw) : FitFunc
91        Wave cw,yw,xw
92       
93#if exists("SpherocylinderX")
94        MultiThread yw = SpherocylinderX(cw,xw)
95#else
96        yw = fSpherocylinder(cw,xw)
97#endif
98        return(0)
99End
100
101//
102// - a double integral - choose points wisely - 76 for both...
103//
104Function fSpherocylinder(w,x) : FitFunc
105        Wave w
106        Variable x
107//       Input (fitting) variables are:
108        //[0] scale factor
109        //[1] cylinder radius (little r)
110        //[2] cylinder length (big L)
111        //[3] end cap radius (big R)
112        //[4] sld cylinder (A^-2)
113        //[5] sld solvent
114        //[6] incoherent background (cm^-1)
115//      give them nice names
116        Variable scale,contr,bkg,inten,sldc,slds
117        Variable len,rad,hDist,endRad
118        scale = w[0]
119        rad = w[1]
120        len = w[2]
121//      endRad = w[3]
122        sldc = w[3]
123        slds = w[4]
124        bkg = w[5]
125       
126        endRad = rad
127       
128        Make/O/D/N=7 SphCyl_tmp
129        SphCyl_tmp[0] = w[0]
130        SphCyl_tmp[1] = w[1]
131        SphCyl_tmp[2] = w[2]
132        SphCyl_tmp[3] = w[1]            //end radius is same as cylinder radius
133        SphCyl_tmp[4] = w[3]
134        SphCyl_tmp[5] = w[4]
135        SphCyl_tmp[6] = w[5]
136       
137        hDist = 0               //by definition
138               
139        contr = sldc-slds
140       
141        Variable/G root:gDumTheta=0,root:gDumT=0
142       
143        inten = IntegrateFn76(SphCyl_Outer,0,pi/2,SphCyl_tmp,x)
144       
145        inten /= pi*rad*rad*len + pi*4*endRad^3/3               //divide by volume
146        inten *= 1e8            //convert to cm^-1
147        inten *= contr*contr
148        inten *= scale
149        inten += bkg
150       
151        Return (inten)
152End
153
154// outer integral
155// x is the q-value
156Function SphCyl_Outer(w,x,dum)
157        Wave w
158        Variable x,dum
159       
160        Variable retVal
161        Variable scale,contr,bkg,inten,sldc,slds
162        Variable len,rad,hDist,endRad,be
163        scale = w[0]
164        rad = w[1]
165        len = w[2]
166        endRad = w[3]
167        sldc = w[4]
168        slds = w[5]
169        bkg = w[6]
170
171        hDist = 0
172                               
173        NVAR dTheta = root:gDumTheta
174        NVAR dt = root:gDumT
175        dTheta = dum
176        retval = IntegrateFn76(SphCyl_Inner,-hDist/endRad,1,w,x)
177       
178        Variable arg1,arg2
179        arg1 = x*len/2*cos(dum)
180        arg2 = x*rad*sin(dum)
181       
182        if(arg2 == 0)
183                be = 0.5
184        else
185                be = Besselj(1, arg2)/arg2
186        endif
187       
188        retVal += pi*rad*rad*len*sinc(arg1)*2*be
189       
190        retVal *= retval*sin(dum)               // = |A(q)|^2*sin(theta)
191       
192        return(retVal)
193End
194
195//returns the value of the integrand of the inner integral
196Function SphCyl_Inner(w,x,dum)
197        Wave w
198        Variable x,dum
199       
200        Variable retVal
201        Variable scale,contr,bkg,inten,sldc,slds
202        Variable len,rad,hDist,endRad
203        scale = w[0]
204        rad = w[1]
205        len = w[2]
206        endRad = w[3]
207        sldc = w[4]
208        slds = w[5]
209        bkg = w[6]
210       
211        NVAR dTheta = root:gDumTheta
212        NVAR dt = root:gDumT
213        dt = dum
214       
215        retVal = SphCyl(w,x,dt,dTheta)
216       
217        retVal *= 4*pi*endRad^3
218       
219        return(retVal)
220End
221
222Function SphCyl(w,x,tt,Theta)
223        Wave w
224        Variable x,tt,Theta
225       
226        Variable val,arg1,arg2
227        Variable scale,contr,bkg,inten,sldc,slds
228        Variable len,rad,hDist,endRad,be
229        scale = w[0]
230        rad = w[1]
231        len = w[2]
232        endRad = w[3]
233        sldc = w[4]
234        slds = w[5]
235        bkg = w[6]
236       
237        hDist = 0
238               
239        arg1 = x*cos(theta)*(endRad*tt+hDist+len/2)
240        arg2 = x*endRad*sin(theta)*sqrt(1-tt*tt)
241       
242        if(arg2 == 0)
243                be = 0.5
244        else
245                be = Besselj(1, arg2)/arg2
246        endif
247       
248        val = cos(arg1)*(1-tt*tt)*be
249       
250        return(val)
251end
252
253//wrapper to calculate the smeared model as an AAO-Struct
254// fills the struct and calls the ususal function with the STRUCT parameter
255//
256// used only for the dependency, not for fitting
257//
258Function fSmearedSpherocylinder(coefW,yW,xW)
259        Wave coefW,yW,xW
260       
261        String str = getWavesDataFolder(yW,0)
262        String DF="root:"+str+":"
263       
264        WAVE resW = $(DF+str+"_res")
265       
266        STRUCT ResSmearAAOStruct fs
267        WAVE fs.coefW = coefW   
268        WAVE fs.yW = yW
269        WAVE fs.xW = xW
270        WAVE fs.resW = resW
271       
272        Variable err
273        err = SmearedSpherocylinder(fs)
274       
275        return (0)
276End
277               
278// this is all there is to the smeared calculation!
279//
280// 20 points should be fine here. This function is not much different than cylinders, where 20 is sufficient
281Function SmearedSpherocylinder(s) :FitFunc
282        Struct ResSmearAAOStruct &s
283
284//      the name of your unsmeared model (AAO) is the first argument
285        Smear_Model_20(Spherocylinder,s.coefW,s.xW,s.yW,s.resW)
286
287        return(0)
288End
Note: See TracBrowser for help on using the repository browser.