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

Revision 633, 7.2 KB checked in by srkline, 3 years ago (diff)

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.

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