source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_WhiteBeamSmear.ipf @ 1090

Last change on this file since 1090 was 1090, checked in by srkline, 5 years ago

Adding proceures - not complete yet- to allow smearing of data under white beam conditions. currently only wavelength effects taken into account, not the geometry. Unknown if geometry is a significant correction given the width of the qhite beam distribution. further testing is needed.

File size: 8.5 KB
Line 
1#pragma TextEncoding = "MacRoman"
2#pragma rtGlobals=3             // Use modern global access method and strict wave access.
3
4//
5//
6// testing routines to compare various integration methods and approximations
7// for calculating the resolution smearing from the white beam wavelength distribution
8//
9//
10//
11// IntegrateFn_N is something that I wrote (in GaussUtils) for quadrature with any number of
12//  points (user-selected)
13//
14// 2018:
15// my quadrature and the built-in function are equivalent. Romberg may be useful in some cases
16// especially for multiple integrals. then number of points and timing can be optimized. But either
17// method can be used. 
18//     
19//              answer = IntegrateFn_N(V_WB_testKernel,loLim,upLim,cw,qVals,nord)
20//      answer_Rom_WB = Integrate_BuiltIn(cw,loLim,upLim,qVals)
21
22// using a matrix multiplication for this calculation of the white beam wavelength smearing is NOT
23// recommended -- the calculation is not nearly accurate enough.
24//
25//
26// Using my built-in quadrature routines (see V_TestWavelengthIntegral) may be of use when
27// writing fitting functions for all of these cases. The built-in Integrate may be limited
28//
29// TODO -- beware what might happen to the calculations since there is a single global string
30//   containing the function name.
31//
32// TODO:
33// -- a significant problem with using the coef waves that are used in the wrapper are that
34//   they are set up with a dependency, so doing the WB calculation also does the "regular"
35//   calculation, doubling the time required...
36
37
38//
39// SANSModel_proto(w,x)         is in GaussUtils_v40.ipf
40//
41// FUNCREF SANSModel_proto fcn
42//
43
44
45
46// call the calculation
47// see DoTheFitButton in Wrapper_v40.ipf
48//
49//
50Macro V_Calc_WB_Smearing()
51
52        String folderStr,funcStr,coefStr
53       
54        ControlInfo/W=WrapperPanel popup_0
55        folderStr=S_Value
56       
57        ControlInfo/W=WrapperPanel popup_1
58        funcStr=S_Value
59       
60        ControlInfo/W=WrapperPanel popup_2
61        coefStr=S_Value
62
63        V_DoWavelengthIntegral(folderStr,funcStr,coefStr)
64
65        SetDataFolder root:
66End
67
68
69// uses built-in Integrate1d()
70//
71Function V_DoWavelengthIntegral(folderStr,funcStr,coefStr)
72        String folderStr,funcStr,coefStr
73
74        SetDataFolder $("root:"+folderStr)
75               
76        // gather the input waves
77        WAVE qVals = $(folderStr+"_q")
78//      WAVE cw = smear_coef_BroadPeak
79        WAVE cw = $coefStr
80       
81        funcStr = V_getXFuncStrFromCoef(cw)+"_"         //get the modelX name, tag on "_"
82        String/G root:gFunctionString = funcStr         // need a global reference to pass to Integrate1D
83       
84        // make a wave for the answer
85        Duplicate/O qvals answer_Rom_WB
86       
87        // do the integration
88        Variable loLim,upLim
89               
90        // define limits based on lo/mean, hi/mean of the wavelength distribution
91        // using the empirical definition, "top" of the peaks
92        loLim = 3.37/5.3
93        upLim = 8.25/5.3
94       
95//      // using the "middle"
96//      loLim = 3.37/5.3
97//      upLim = 8.37/5.3       
98//     
99//      // using the interpolated distribution (change the function call)
100//      lolim = 3/5.3
101//      uplim = 9/5.3
102
103// using the "trangular" distribution (cange the function call)
104//      loLim = 4/5.3
105//      upLim = 8/5.3
106       
107        answer_Rom_WB = Integrate_BuiltIn(cw,loLim,upLim,qVals)
108
109// why do I need this? Is this because this is defined as the mean of the distribution
110//  and is needed to normalize the integral? verify this on paper.     
111        answer_Rom_WB *= 5.3
112
113// normalize the integral       
114        answer_Rom_WB /= 20926          // "top"  of peaks
115//      answer_Rom_WB /= 19933          // "middle"  of peaks
116//      answer_Rom_WB /= 20051          // interpolated distribution
117//      answer_Rom_WB /= 1              // triangular distribution (it's already normalized)
118       
119        SetDataFolder root:
120       
121        return 0
122End
123
124//
125// not used anymore - the built-in works fine, but this
126// may be of use if I convert all of these to fitting functions.
127//
128Function V_TestWavelengthIntegral(folderStr)
129        String folderStr
130
131        SetDataFolder $("root:"+folderStr)
132       
133       
134        // gather the input waves
135        WAVE qVals = $(folderStr+"_q")
136//      WAVE cw = smear_coef_sf
137//      WAVE cw = smear_coef_pgs
138        WAVE cw = smear_coef_BroadPeak
139       
140        // make a wave for the answer
141//      Duplicate/O qvals answer, answer_builtIn
142        Duplicate/O qvals answer_Quad
143       
144        // do the integration
145        // Function IntegrateFn_N(fcn,loLim,upLim,w,x,nord)                             
146
147        Variable loLim,upLim,nord
148       
149        nord = 76       // 20 quadrature points not enough for white beam (especially AgBeh test)
150       
151        loLim = 4/5.3
152        upLim = 8/5.3
153
154// 2018:
155// my quadrature and the built-in function are equivalent. Romberg may be useful in some cases
156// especially for multiple integrals. then number of points and timing can be optimized. But either
157// method can be used. 
158        answer_Quad = IntegrateFn_N(V_WB_testKernel,loLim,upLim,cw,qVals,nord)
159       
160
161// why do I need this? Is this because this is defined as the mean of the distribution
162//  and is needed to normalize the integral? verify this on paper.     
163        answer_Quad *= 5.3
164//      answer_builtIn *= 5.3
165       
166        SetDataFolder root:
167       
168        return 0
169End
170
171
172
173Function V_WB_testKernel(cw,x,dum)
174        Wave cw
175        Variable x              // the q-value for the calculation
176        Variable dum    // the dummy integration variable
177
178        Variable val
179        SVAR funcStr = root:gFunctionString
180        FUNCREF SANSModel_proto func = $funcStr
181               
182//      val = (1-dum*5.3/8)*BroadPeakX(cw,x/dum)
183        val = (1-dum*5.3/8)*func(cw,x/dum)
184       
185//      val = V_WhiteBeamDist(dum*5.3)*BroadPeakX(cw,x/dum)
186        val = V_WhiteBeamDist(dum*5.3)*func(cw,x/dum)
187
188        return (val)
189End
190
191Proc WBDistr()
192
193        make/O/D distr
194        SetScale/I x 0.755,1.509,"", distr
195        distr = (1-x*5.3/8)
196        display distr
197
198end
199
200// the trick here is that declaring the last qVal wave as a variable
201// since this is implicitly called N times in the wave assignment of the answer wave
202Function Integrate_BuiltIn(cw,loLim,upLim,qVal)
203        Wave cw
204        Variable loLim,upLim
205        Variable qVal
206       
207        Variable/G root:qq = qval
208        Variable ans
209       
210//      ans = Integrate1D(intgrnd,lolim,uplim,2,0,cw)           //adaptive quadrature
211        ans = Integrate1D(intgrnd,lolim,uplim,1,0,cw)           // Romberg integration
212       
213        return ans
214end
215
216
217//
218// See V_DunnyFunctions.ipf for the full list
219//
220//Function BroadPeakX_(cw,x)
221//      Wave cw
222//      Variable x
223//     
224//      return(BroadPeakX(cw,x))
225//end
226
227Function intgrnd(cw,dum)
228        Wave cw
229        Variable dum            // the dummy of the integration
230
231        Variable val
232        NVAR qq = root:qq               //the q-value of the integration, not part of cw, so pass global
233        SVAR funcStr = root:gFunctionString
234        FUNCREF SANSModel_proto func = $funcStr
235
236//      val = (1-dum*5.3/8)*BroadPeakX(cw,qq/dum)       
237//      val = (1-dum*5.3/8)*func(cw,qq/dum)
238
239//      val = V_WhiteBeamDist(dum*5.3)*BroadPeakX(cw,qq/dum)
240        val = V_WhiteBeamDist(dum*5.3)*func(cw,qq/dum)
241       
242//      val = V_WhiteBeamInterp(dum*5.3)*func(cw,qq/dum)
243
244        return (val)
245End
246
247
248////////////////////////////
249
250// need a function to return the model function name
251// given the coefficient wave
252//
253// want the function NameX for use in the integration, not the AAO function
254//
255
256// from the name of the coefficient wave, get the function name
257// be sure that there is no "Smeared" at the beginning of the name
258// tag X to the end of the name string
259//
260// then the funcString must be passed in as a global to the built-in integration function.
261//
262Function/S V_getXFuncStrFromCoef(cw)
263        Wave cw
264       
265        String cwStr = NameOfWave(cw)
266        String outStr = "",extStr=""
267       
268//      String convStr = ReplaceString("_",cwStr,".")           // change the _ to .   
269//      extStr = ParseFilePath(4, convStr, ":", 0, 0)           // extracts the last .nnn, without the .
270
271// go through the list of coefKWStr pairs
272// look for the cwStr
273// take up to the = (that is the funcStr)
274// remove "Smeared" if needed   
275        SVAR coefList=root:Packages:NIST:coefKWStr
276
277        Variable ii,num
278        String item
279       
280        num=ItemsInList(coefList,";")
281        ii=0
282        do
283                item = StringFromList(ii, coefList, ";")
284               
285                if(strsearch(item,cwStr,0) != -1)               //match
286                        item = ReplaceString("=",item,".")              //replace the = with .
287                        outStr = ParseFilePath(3, item, ":", 0, 0)              // extract file name without extension
288                        outStr = ReplaceString("Smeared",outStr,"")             // replace "Smeared" with null, if it's there
289                        ii = num + 1
290                endif
291               
292                ii+=1
293        while(ii<num)
294       
295        return(outStr+"X")
296end
297
298//////////////////////////////////////////
299// generates dummy functions of the form:
300//
301//Function BroadPeakX_(cw,x)
302//      Wave cw
303//      Variable x
304//      return(BroadPeakX(cw,x))
305//End
306//
307// so that I can use the FUNCREF
308// which fails for some reason when I just use the XOP name?
309//
310//
311// not everything ending in X is a model function - trimmed list is in V_DummyFunctions.ipf
312//
313Function V_generateDummyFuncs()
314
315        String list = FunctionList("*X",";","KIND:4")
316        Variable ii,num
317        String item,str
318       
319        num=ItemsInList(list,";")
320
321        NewNotebook/N=Notebook1/F=0
322       
323       
324        for(ii=0;ii<num;ii+=1)
325                item = StringFromList(ii,list,";")
326                str = "\r"
327                str = "Function "+item+"_(cw,x)\r"
328                str += "\tWave cw\r"
329                str += "\tVariable x\r"
330                str += "\treturn("+item+"(cw,x))\r"
331                str += "End\r\r"
332               
333                //print str
334       
335                Notebook $"", text=str
336               
337        endfor
338        return(0)
339       
340End
Note: See TracBrowser for help on using the repository browser.