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

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

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

File size: 6.8 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 "barbell" shape - a cylinder with spherical end caps
7// of larger radius than 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 PlotBarbell(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_Barbell, ywave_Barbell
31        xwave_Barbell =  alog(log(qmin) + x*((log(qmax)-log(qmin))/num))
32        Make/O/D coef_Barbell = {1,20,400,40,1e-6,6.3e-6,0}                     //CH#2
33        make/o/t parameters_Barbell = {"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
34        Edit parameters_Barbell, coef_Barbell
35       
36        Variable/G root:g_Barbell
37        g_Barbell := Barbell(coef_Barbell, ywave_Barbell, xwave_Barbell)
38        Display ywave_Barbell vs xwave_Barbell
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("Barbell","coef_Barbell","Barbell")
46//
47End
48
49
50// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
51Proc PlotSmearedBarbell(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_Barbell = {1,20,400,40,1e-6,6.3e-6,0}               //CH#4
64        make/o/t smear_parameters_Barbell = {"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)"}
65        Edit smear_parameters_Barbell,smear_coef_Barbell                                        //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_Barbell,smeared_qvals                           //
70        SetScale d,0,0,"1/cm",smeared_Barbell                                                   //
71                                       
72        Variable/G gs_Barbell=0
73        gs_Barbell := fSmearedBarbell(smear_coef_Barbell,smeared_Barbell,smeared_qvals) //this wrapper fills the STRUCT
74       
75        Display smeared_Barbell 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("SmearedBarbell","smear_coef_Barbell","Barbell")
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 Barbell(cw,yw,xw) : FitFunc
91        Wave cw,yw,xw
92       
93#if exists("BarbellX")
94        yw = BarbellX(cw,xw)
95#else
96        yw = fBarbell(cw,xw)
97#endif
98        return(0)
99End
100
101//
102// - a double integral - choose points wisely - 76 for both...
103//
104Function fBarbell(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[4]
123        slds = w[5]
124        bkg = w[6]
125       
126        hDist = sqrt(endRad^2-rad^2)   
127               
128        contr = sldc-slds
129       
130        Variable/G root:gDumTheta=0,root:gDumT=0
131       
132        inten = IntegrateFn76(Barbell_Outer,0,pi/2,w,x)
133       
134        inten /= pi*rad*rad*len + 2*pi*(2*endRad^3/3+endRad^2*hDist-hDist^3/3)          //divide by volume
135        inten *= 1e8            //convert to cm^-1
136        inten *= contr*contr
137        inten *= scale
138        inten += bkg
139       
140        Return (inten)
141End
142
143// outer integral
144// x is the q-value
145Function Barbell_Outer(w,x,dum)
146        Wave w
147        Variable x,dum
148       
149        Variable retVal
150        Variable scale,contr,bkg,inten,sldc,slds
151        Variable len,rad,hDist,endRad
152        scale = w[0]
153        rad = w[1]
154        len = w[2]
155        endRad = w[3]
156        sldc = w[4]
157        slds = w[5]
158        bkg = w[6]
159
160        hDist = sqrt(endRad^2-rad^2)   
161                       
162        NVAR dTheta = root:gDumTheta
163        NVAR dt = root:gDumT
164        dTheta = dum
165        retval = IntegrateFn76(Barbell_Inner,-hDist/endRad,1,w,x)
166       
167        Variable arg1,arg2
168        arg1 = x*len/2*cos(dum)
169        arg2 = x*rad*sin(dum)
170       
171        retVal += pi*rad*rad*len*sinc(arg1)*2*Besselj(1, arg2)/arg2
172       
173        retVal *= retval*sin(dum)               // = |A(q)|^2*sin(theta)
174       
175        return(retVal)
176End
177
178//returns the value of the integrand of the inner integral
179Function Barbell_Inner(w,x,dum)
180        Wave w
181        Variable x,dum
182       
183        Variable retVal
184        Variable scale,contr,bkg,inten,sldc,slds
185        Variable len,rad,hDist,endRad
186        scale = w[0]
187        rad = w[1]
188        len = w[2]
189        endRad = w[3]
190        sldc = w[4]
191        slds = w[5]
192        bkg = w[6]
193       
194        NVAR dTheta = root:gDumTheta
195        NVAR dt = root:gDumT
196        dt = dum
197       
198        retVal = Barbell_integrand(w,x,dt,dTheta)
199       
200        retVal *= 4*pi*endRad^3
201       
202        return(retVal)
203End
204
205Function Barbell_integrand(w,x,tt,Theta)
206        Wave w
207        Variable x,tt,Theta
208       
209        Variable val,arg1,arg2
210        Variable scale,contr,bkg,inten,sldc,slds
211        Variable len,rad,hDist,endRad
212        scale = w[0]
213        rad = w[1]
214        len = w[2]
215        endRad = w[3]
216        sldc = w[4]
217        slds = w[5]
218        bkg = w[6]
219
220        hDist = sqrt(endRad^2-rad^2)   
221               
222        arg1 = x*cos(theta)*(endRad*tt+hDist+len/2)
223        arg2 = x*endRad*sin(theta)*sqrt(1-tt*tt)
224       
225        val = cos(arg1)*(1-tt*tt)*Besselj(1,arg2)/arg2
226       
227        return(val)
228end
229
230//wrapper to calculate the smeared model as an AAO-Struct
231// fills the struct and calls the ususal function with the STRUCT parameter
232//
233// used only for the dependency, not for fitting
234//
235Function fSmearedBarbell(coefW,yW,xW)
236        Wave coefW,yW,xW
237       
238        String str = getWavesDataFolder(yW,0)
239        String DF="root:"+str+":"
240       
241        WAVE resW = $(DF+str+"_res")
242       
243        STRUCT ResSmearAAOStruct fs
244        WAVE fs.coefW = coefW   
245        WAVE fs.yW = yW
246        WAVE fs.xW = xW
247        WAVE fs.resW = resW
248       
249        Variable err
250        err = SmearedBarbell(fs)
251       
252        return (0)
253End
254               
255// this is all there is to the smeared calculation!
256//
257// 20 points should be fine here. This function is not much different than cylinders, where 20 is sufficient
258Function SmearedBarbell(s) :FitFunc
259        Struct ResSmearAAOStruct &s
260
261//      the name of your unsmeared model (AAO) is the first argument
262        Smear_Model_20(Barbell,s.coefW,s.xW,s.yW,s.resW)
263
264        return(0)
265End
Note: See TracBrowser for help on using the repository browser.