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