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

Revision 449, 7.0 KB checked in by srkline, 5 years ago (diff)

Adding the model of a cylinder with spherical end caps. By taking all of the limiting cases, 5 models are added. Since the geometries are sufficiently different, I thought it best to add as separate models. Corresponding changes have been made to the Wrapper and the modelPicker

Line
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma IgorVersion = 6.0
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")
97        yw = DumbbellX(cw,xw)
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
185
186        retVal *= retval*sin(dum)               // = |A(q)|^2*sin(theta)
187
188        return(retVal)
189End
190
191//returns the value of the integrand of the inner integral
192Function Dumb_Inner(w,x,dum)
193        Wave w
194        Variable x,dum
195
196        Variable retVal
197        Variable scale,contr,bkg,inten,sldc,slds
199        scale = w[0]
201        len = w[2]
203        sldc = w[4]
204        slds = w[5]
205        bkg = w[6]
207
208        NVAR dTheta = root:gDumTheta
209        NVAR dt = root:gDumT
210        dt = dum
211
212        retVal = Dumb(w,x,dt,dTheta)
213
215
216        return(retVal)
217End
218
219Function Dumb(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 fSmearedDumbbell(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 = SmearedDumbbell(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 SmearedDumbbell(s) :FitFunc
273        Struct ResSmearAAOStruct &s
274
275//      the name of your unsmeared model (AAO) is the first argument
276        Smear_Model_20(Dumbbell,s.coefW,s.xW,s.yW,s.resW)
277
278        return(0)
279End
Note: See TracBrowser for help on using the repository browser.