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

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

Simple change in all of the model function files to include the name of the parameter wave in the Keyword=list that is generated when a model is plotted. This is becoming an issue where the proper parameter wave can't be deduced from just the suffix, then there is nothing to put in the table.

I should have added this when I initially wrote the wrapper...

File size: 7.1 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 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//
15// 76 point quadrature is necessary for both quadrature calls.
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       
48        AddModelToStrings("Dumbbell","coef_Dumb","parameters_Dumb","Dumb")
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:
85        AddModelToStrings("SmearedDumbbell","smear_coef_Dumb","smear_parameters_Dumb","Dumb")
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
120        Variable len,rad,hDist,endRad
121        scale = w[0]
122        rad = w[1]
123//      len = w[2]
124        endRad = w[2]
125        sldc = w[3]
126        slds = w[4]
127        bkg = w[5]
128       
129        hDist = sqrt(endRad^2-rad^2)   
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       
147        inten /= 2*pi*(2*endRad^3/3+endRad^2*hDist-hDist^3/3)           //divide by volume
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
164        Variable len,rad,hDist,endRad
165        scale = w[0]
166        rad = w[1]
167        len = w[2]
168        endRad = w[3]
169        sldc = w[4]
170        slds = w[5]
171        bkg = w[6]
172
173        hDist = sqrt(endRad^2-rad^2)
174                       
175        NVAR dTheta = root:gDumTheta
176        NVAR dt = root:gDumT
177        dTheta = dum
178        retval = IntegrateFn76(Dumb_Inner,-hDist/endRad,1,w,x)
179       
180        Variable arg1,arg2
181        arg1 = x*len/2*cos(dum)
182        arg2 = x*rad*sin(dum)
183       
184        retVal += pi*rad*rad*len*sinc(arg1)*2*Besselj(1, arg2)/arg2
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
198        Variable len,rad,hDist,endRad
199        scale = w[0]
200        rad = w[1]
201        len = w[2]
202        endRad = w[3]
203        sldc = w[4]
204        slds = w[5]
205        bkg = w[6]
206//      hDist = sqrt(endRad^2-rad^2)   
207       
208        NVAR dTheta = root:gDumTheta
209        NVAR dt = root:gDumT
210        dt = dum
211       
212        retVal = Dumb(w,x,dt,dTheta)
213       
214        retVal *= 4*pi*endRad^3
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
225        Variable len,rad,hDist,endRad
226        scale = w[0]
227        rad = w[1]
228        len = w[2]
229        endRad = w[3]
230        sldc = w[4]
231        slds = w[5]
232        bkg = w[6]
233
234        hDist = sqrt(endRad^2-rad^2)   
235               
236        arg1 = x*cos(theta)*(endRad*tt+hDist+len/2)
237        arg2 = x*endRad*sin(theta)*sqrt(1-tt*tt)
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.