# 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//
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
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:
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
119        scale = w[0]
121        len = w[2]
123        sldc = w[4]
124        slds = w[5]
125        bkg = w[6]
126
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
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
154        scale = w[0]
156        len = w[2]
158        sldc = w[4]
159        slds = w[5]
160        bkg = w[6]
161
163
164        NVAR dTheta = root:gDumTheta
165        NVAR dt = root:gDumT
166        dTheta = dum
168
169        Variable arg1,arg2
170        arg1 = x*len/2*cos(dum)
172
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
188        scale = w[0]
190        len = w[2]
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
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
214        scale = w[0]
216        len = w[2]
218        sldc = w[4]
219        slds = w[5]
220        bkg = w[6]
221
223
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.