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

Last change on this file since 1242 was 1242, checked in by srkline, 3 years ago

updating the IgorVersion? pragma to v7.0 for all files to be consistent.

File size: 17.8 KB
Line 
1#pragma TextEncoding = "MacRoman"
2#pragma rtGlobals=3             // Use modern global access method and strict wave access.
3#pragma IgorVersion = 7.00
4
5//
6//
7// testing routines to compare various integration methods and approximations
8// for calculating the resolution smearing from the white beam wavelength distribution
9//
10// none of these functions are used in the final version of the resolution smearing
11// for white beam or super white beam. Functions and definitions of WB and SWB are contatined
12// in the file: V_WhiteBeamDistribution.ipf
13//
14//
15//
16//
17//
18// IntegrateFn_N is something that I wrote (in GaussUtils) for quadrature with any number of
19//  points (user-selected)
20//
21// 2018:
22// my quadrature and the built-in function are equivalent. Romberg may be useful in some cases
23// especially for multiple integrals. then number of points and timing can be optimized. But either
24// method can be used. 
25//     
26//              answer = IntegrateFn_N(V_WB_testKernel,loLim,upLim,cw,qVals,nord)
27//      answer_Rom_WB = Integrate_BuiltIn(cw,loLim,upLim,qVals)
28
29// using a matrix multiplication for this calculation of the white beam wavelength smearing is NOT
30// recommended -- the calculation is not nearly accurate enough.
31//
32//
33// Using my built-in quadrature routines (see V_TestWavelengthIntegral) may be of use when
34// writing fitting functions for all of these cases. The built-in Integrate may be limited
35//
36// NOTE: -- beware what might happen to the calculations since there is a single global string
37//   containing the function name.
38//
39// NOTE:
40// -- a significant problem with using the coef waves that are used in the wrapper are that
41//   they are set up with a dependency, so doing the WB calculation also does the "regular"
42//   smeared calculation, doubling the time required...
43
44
45//
46// needs V_DummyFunctions for the FUNCREF to work - since it fails if I simply call the XOP
47//
48//
49// SANSModel_proto(w,x)         is in GaussUtils_v40.ipf
50//
51// FUNCREF SANSModel_proto fcn
52//
53
54
55
56// call the calculation
57// see DoTheFitButton in Wrapper_v40.ipf
58//
59//
60Proc V_Calc_WB_Smearing_top()
61
62        String folderStr,funcStr,coefStr
63       
64        ControlInfo/W=WrapperPanel popup_0
65        folderStr=S_Value
66       
67        ControlInfo/W=WrapperPanel popup_1
68        funcStr=S_Value
69       
70        ControlInfo/W=WrapperPanel popup_2
71        coefStr=S_Value
72
73        V_DoWavelengthIntegral_top(folderStr,funcStr,coefStr)
74
75        SetDataFolder root:
76End
77
78
79Proc V_Calc_WB_Smearing_mid()
80
81        String folderStr,funcStr,coefStr
82       
83        ControlInfo/W=WrapperPanel popup_0
84        folderStr=S_Value
85       
86        ControlInfo/W=WrapperPanel popup_1
87        funcStr=S_Value
88       
89        ControlInfo/W=WrapperPanel popup_2
90        coefStr=S_Value
91
92        V_DoWavelengthIntegral_mid(folderStr,funcStr,coefStr)
93
94        SetDataFolder root:
95End
96
97Proc V_Calc_WB_Smearing_interp()
98
99        String folderStr,funcStr,coefStr
100       
101        ControlInfo/W=WrapperPanel popup_0
102        folderStr=S_Value
103       
104        ControlInfo/W=WrapperPanel popup_1
105        funcStr=S_Value
106       
107        ControlInfo/W=WrapperPanel popup_2
108        coefStr=S_Value
109
110        V_DoWavelengthIntegral_interp(folderStr,funcStr,coefStr)
111
112        SetDataFolder root:
113End
114
115Proc V_Calc_WB_Smearing_triang()
116
117        String folderStr,funcStr,coefStr
118       
119        ControlInfo/W=WrapperPanel popup_0
120        folderStr=S_Value
121       
122        ControlInfo/W=WrapperPanel popup_1
123        funcStr=S_Value
124       
125        ControlInfo/W=WrapperPanel popup_2
126        coefStr=S_Value
127
128        V_DoWavelengthIntegral_triang(folderStr,funcStr,coefStr)
129
130        SetDataFolder root:
131End
132
133
134
135// uses built-in Integrate1d()
136//
137Function V_DoWavelengthIntegral_top(folderStr,funcStr,coefStr)
138        String folderStr,funcStr,coefStr
139
140        SetDataFolder $("root:"+folderStr)
141               
142        // gather the input waves
143        WAVE qVals = $(folderStr+"_q")
144//      WAVE cw = smear_coef_BroadPeak
145        WAVE cw = $coefStr
146       
147        funcStr = V_getXFuncStrFromCoef(cw)+"_"         //get the modelX name, tag on "_"
148        String/G root:gFunctionString = funcStr         // need a global reference to pass to Integrate1D
149       
150        // make a wave for the answer
151        Duplicate/O qvals answer_Rom_WB_top
152       
153        // do the integration
154        Variable loLim,upLim
155               
156        // define limits based on lo/mean, hi/mean of the wavelength distribution
157        // using the empirical definition, "top" of the peaks
158        loLim = 3.37/5.3
159        upLim = 8.25/5.3
160       
161//      // using the "middle"
162//      loLim = 3.37/5.3
163//      upLim = 8.37/5.3       
164//     
165//      // using the interpolated distribution (must change the function call)
166//      lolim = 3/5.3
167//      uplim = 9/5.3
168
169// using the "triangular" distribution (must change the function call)
170//      loLim = 4/5.3
171//      upLim = 8/5.3
172       
173        answer_Rom_WB_top = V_Integrate_BuiltIn_top(cw,loLim,upLim,qVals)
174
175// why do I need this? Is this because this is defined as the mean of the distribution
176//  and is needed to normalize the integral? verify this on paper.     
177        answer_Rom_WB_top *= 5.3
178
179// normalize the integral       
180        answer_Rom_WB_top /= 20926              // "top"  of peaks
181//      answer_Rom_WB /= 19933          // "middle"  of peaks
182//      answer_Rom_WB /= 20051          // interpolated distribution
183//      answer_Rom_WB /= 1              // triangular distribution (it's already normalized)
184
185// additional normalization???
186        answer_Rom_WB_top /= 1.05               //
187       
188        SetDataFolder root:
189       
190        return 0
191End
192
193
194// uses built-in Integrate1d()
195//
196Function V_DoWavelengthIntegral_mid(folderStr,funcStr,coefStr)
197        String folderStr,funcStr,coefStr
198
199        SetDataFolder $("root:"+folderStr)
200               
201        // gather the input waves
202        WAVE qVals = $(folderStr+"_q")
203//      WAVE cw = smear_coef_BroadPeak
204        WAVE cw = $coefStr
205       
206        funcStr = V_getXFuncStrFromCoef(cw)+"_"         //get the modelX name, tag on "_"
207        String/G root:gFunctionString = funcStr         // need a global reference to pass to Integrate1D
208       
209        // make a wave for the answer
210        Duplicate/O qvals answer_Rom_WB_mid
211       
212        // do the integration
213        Variable loLim,upLim
214               
215        // define limits based on lo/mean, hi/mean of the wavelength distribution
216        // using the empirical definition, "top" of the peaks
217//      loLim = 3.37/5.3
218//      upLim = 8.25/5.3
219       
220//      // using the "middle"
221        loLim = 3.37/5.3
222        upLim = 8.37/5.3       
223//     
224//      // using the interpolated distribution (must change the function call)
225//      lolim = 3/5.3
226//      uplim = 9/5.3
227
228// using the "triangular" distribution (must change the function call)
229//      loLim = 4/5.3
230//      upLim = 8/5.3
231       
232        answer_Rom_WB_mid = V_Integrate_BuiltIn_mid(cw,loLim,upLim,qVals)
233
234// why do I need this? Is this because this is defined as the mean of the distribution
235//  and is needed to normalize the integral? verify this on paper.     
236        answer_Rom_WB_mid *= 5.3
237
238// normalize the integral       
239//      answer_Rom_WB /= 20926          // "top"  of peaks
240        answer_Rom_WB_mid /= 19933              // "middle"  of peaks
241//      answer_Rom_WB /= 20051          // interpolated distribution
242//      answer_Rom_WB /= 1              // triangular distribution (it's already normalized)
243
244// additional normalization???
245        answer_Rom_WB_mid /= 1.05               // "middle"  of peaks
246       
247        SetDataFolder root:
248       
249        return 0
250End
251
252// uses built-in Integrate1d()
253//
254Function V_DoWavelengthIntegral_interp(folderStr,funcStr,coefStr)
255        String folderStr,funcStr,coefStr
256
257        SetDataFolder $("root:"+folderStr)
258               
259        // gather the input waves
260        WAVE qVals = $(folderStr+"_q")
261//      WAVE cw = smear_coef_BroadPeak
262        WAVE cw = $coefStr
263       
264        funcStr = V_getXFuncStrFromCoef(cw)+"_"         //get the modelX name, tag on "_"
265        String/G root:gFunctionString = funcStr         // need a global reference to pass to Integrate1D
266       
267        // make a wave for the answer
268        Duplicate/O qvals answer_Rom_WB_interp
269       
270        // do the integration
271        Variable loLim,upLim
272               
273        // define limits based on lo/mean, hi/mean of the wavelength distribution
274        // using the empirical definition, "top" of the peaks
275//      loLim = 3.37/5.3
276//      upLim = 8.25/5.3
277       
278//      // using the "middle"
279//      loLim = 3.37/5.3
280//      upLim = 8.37/5.3       
281//     
282//      // using the interpolated distribution (must change the function call)
283        lolim = 3/5.3
284        uplim = 9/5.3
285
286// using the "triangular" distribution (must change the function call)
287//      loLim = 4/5.3
288//      upLim = 8/5.3
289       
290        answer_Rom_WB_interp = V_Integrate_BuiltIn_interp(cw,loLim,upLim,qVals)
291
292// why do I need this? Is this because this is defined as the mean of the distribution
293//  and is needed to normalize the integral? verify this on paper.     
294        answer_Rom_WB_interp *= 5.3
295
296// normalize the integral       
297//      answer_Rom_WB /= 20926          // "top"  of peaks
298//      answer_Rom_WB /= 19933          // "middle"  of peaks
299        answer_Rom_WB_interp /= 20051           // interpolated distribution
300//      answer_Rom_WB /= 1              // triangular distribution (it's already normalized)
301
302// additional normalization???
303        answer_Rom_WB_interp /= 1.05            // "middle"  of peaks
304       
305        SetDataFolder root:
306       
307        return 0
308End
309
310// uses built-in Integrate1d()
311//
312Function V_DoWavelengthIntegral_triang(folderStr,funcStr,coefStr)
313        String folderStr,funcStr,coefStr
314
315        SetDataFolder $("root:"+folderStr)
316               
317        // gather the input waves
318        WAVE qVals = $(folderStr+"_q")
319//      WAVE cw = smear_coef_BroadPeak
320        WAVE cw = $coefStr
321       
322        funcStr = V_getXFuncStrFromCoef(cw)+"_"         //get the modelX name, tag on "_"
323        String/G root:gFunctionString = funcStr         // need a global reference to pass to Integrate1D
324       
325        // make a wave for the answer
326        Duplicate/O qvals answer_Rom_WB_triang
327       
328        // do the integration
329        Variable loLim,upLim
330               
331        // define limits based on lo/mean, hi/mean of the wavelength distribution
332        // using the empirical definition, "top" of the peaks
333//      loLim = 3.37/5.3
334//      upLim = 8.25/5.3
335       
336//      // using the "middle"
337//      loLim = 3.37/5.3
338//      upLim = 8.37/5.3       
339//     
340//      // using the interpolated distribution (must change the function call)
341//      lolim = 3/5.3
342//      uplim = 9/5.3
343
344// using the "triangular" distribution (must change the function call)
345        loLim = 4/5.3
346        upLim = 8/5.3
347       
348        answer_Rom_WB_triang = V_Integrate_BuiltIn_triangle(cw,loLim,upLim,qVals)
349
350// why do I need this? Is this because this is defined as the mean of the distribution
351//  and is needed to normalize the integral? verify this on paper.     
352        answer_Rom_WB_triang *= 5.3
353
354// normalize the integral       
355//      answer_Rom_WB /= 20926          // "top"  of peaks
356//      answer_Rom_WB /= 19933          // "middle"  of peaks
357//      answer_Rom_WB /= 20051          // interpolated distribution
358        answer_Rom_WB_triang /= 1               // triangular distribution (it's already normalized)
359
360// additional normalization???
361        answer_Rom_WB_triang /= 1.1             //
362       
363        SetDataFolder root:
364       
365        return 0
366End
367
368
369//
370// not used anymore - the built-in works fine, but this
371// may be of use if I convert all of these to fitting functions.
372//
373Function V_TestWavelengthIntegral(folderStr)
374        String folderStr
375
376        SetDataFolder $("root:"+folderStr)
377       
378       
379        // gather the input waves
380        WAVE qVals = $(folderStr+"_q")
381//      WAVE cw = smear_coef_sf
382//      WAVE cw = smear_coef_pgs
383        WAVE cw = smear_coef_BroadPeak
384       
385        // make a wave for the answer
386//      Duplicate/O qvals answer, answer_builtIn
387        Duplicate/O qvals answer_Quad
388       
389        // do the integration
390        // Function IntegrateFn_N(fcn,loLim,upLim,w,x,nord)                             
391
392        Variable loLim,upLim,nord
393       
394        nord = 76       // 20 quadrature points not enough for white beam (especially AgBeh test)
395       
396        loLim = 4/5.3
397        upLim = 8/5.3
398
399// 2018:
400// my quadrature and the built-in function are equivalent. Romberg may be useful in some cases
401// especially for multiple integrals. then number of points and timing can be optimized. But either
402// method can be used. 
403        answer_Quad = IntegrateFn_N(V_WB_testKernel,loLim,upLim,cw,qVals,nord)
404       
405
406// why do I need this? Is this because this is defined as the mean of the distribution
407//  and is needed to normalize the integral? verify this on paper.     
408        answer_Quad *= 5.3
409//      answer_builtIn *= 5.3
410       
411        SetDataFolder root:
412       
413        return 0
414End
415
416
417
418Function V_WB_testKernel(cw,x,dum)
419        Wave cw
420        Variable x              // the q-value for the calculation
421        Variable dum    // the dummy integration variable
422
423        Variable val
424        SVAR funcStr = root:gFunctionString
425        FUNCREF SANSModel_proto func = $funcStr
426               
427//      val = (1-dum*5.3/8)*BroadPeakX(cw,x/dum)
428        val = (1-dum*5.3/8)*func(cw,x/dum)
429       
430//      val = V_WhiteBeamDist(dum*5.3)*BroadPeakX(cw,x/dum)
431//      val = V_WhiteBeamDist(dum*5.3)*func(cw,x/dum)
432
433        return (val)
434End
435
436Proc WBDistr()
437
438        make/O/D distr
439        SetScale/I x 0.755,1.509,"", distr
440        distr = (1-x*5.3/8)
441        display distr
442
443end
444
445// the trick here is that declaring the last qVal wave as a variable
446// since this is implicitly called N times in the wave assignment of the answer wave
447Function V_Integrate_BuiltIn_top(cw,loLim,upLim,qVal)
448        Wave cw
449        Variable loLim,upLim
450        Variable qVal
451       
452        Variable/G root:qq = qval
453        Variable ans
454       
455//      ans = Integrate1D(V_intgrnd_top,lolim,uplim,2,0,cw)             //adaptive quadrature
456        ans = Integrate1D(V_intgrnd_top,lolim,uplim,1,0,cw)             // Romberg integration
457       
458        return ans
459end
460
461// the trick here is that declaring the last qVal wave as a variable
462// since this is implicitly called N times in the wave assignment of the answer wave
463Function V_Integrate_BuiltIn_mid(cw,loLim,upLim,qVal)
464        Wave cw
465        Variable loLim,upLim
466        Variable qVal
467       
468        Variable/G root:qq = qval
469        Variable ans
470       
471//      ans = Integrate1D(V_intgrnd_mid,lolim,uplim,2,0,cw)             //adaptive quadrature
472        ans = Integrate1D(V_intgrnd_mid,lolim,uplim,1,0,cw)             // Romberg integration
473       
474        return ans
475end
476
477// the trick here is that declaring the last qVal wave as a variable
478// since this is implicitly called N times in the wave assignment of the answer wave
479Function V_Integrate_BuiltIn_triangle(cw,loLim,upLim,qVal)
480        Wave cw
481        Variable loLim,upLim
482        Variable qVal
483       
484        Variable/G root:qq = qval
485        Variable ans
486       
487//      ans = Integrate1D(V_intgrnd_triangle,lolim,uplim,2,0,cw)                //adaptive quadrature
488        ans = Integrate1D(V_intgrnd_triangle,lolim,uplim,1,0,cw)                // Romberg integration
489       
490        return ans
491end
492
493// the trick here is that declaring the last qVal wave as a variable
494// since this is implicitly called N times in the wave assignment of the answer wave
495Function V_Integrate_BuiltIn_interp(cw,loLim,upLim,qVal)
496        Wave cw
497        Variable loLim,upLim
498        Variable qVal
499       
500        Variable/G root:qq = qval
501        Variable ans
502       
503//      ans = Integrate1D(V_intgrnd_interp,lolim,uplim,2,0,cw)          //adaptive quadrature
504        ans = Integrate1D(V_intgrnd_interp,lolim,uplim,1,0,cw)          // Romberg integration
505       
506        return ans
507end
508
509//
510// See V_DummyFunctions.ipf for the full list
511//
512//Function BroadPeakX_(cw,x)
513//      Wave cw
514//      Variable x
515//     
516//      return(BroadPeakX(cw,x))
517//end
518
519Function V_intgrnd_top(cw,dum)
520        Wave cw
521        Variable dum            // the dummy of the integration
522
523        Variable val
524        NVAR qq = root:qq               //the q-value of the integration, not part of cw, so pass global
525        SVAR funcStr = root:gFunctionString
526        FUNCREF SANSModel_proto func = $funcStr
527
528//      val = (1-dum*5.3/8)*BroadPeakX(cw,qq/dum)       
529//      val = (1-dum*5.3/8)*func(cw,qq/dum)
530
531//      val = V_WhiteBeamDist(dum*5.3)*BroadPeakX(cw,qq/dum)
532        val = V_WhiteBeamDist_top(dum*5.3)*func(cw,qq/dum)
533       
534//      val = V_WhiteBeamInterp(dum*5.3)*func(cw,qq/dum)
535
536        return (val)
537End
538
539Function V_intgrnd_mid(cw,dum)
540        Wave cw
541        Variable dum            // the dummy of the integration
542
543        Variable val
544        NVAR qq = root:qq               //the q-value of the integration, not part of cw, so pass global
545        SVAR funcStr = root:gFunctionString
546        FUNCREF SANSModel_proto func = $funcStr
547
548//      val = (1-dum*5.3/8)*BroadPeakX(cw,qq/dum)       
549//      val = (1-dum*5.3/8)*func(cw,qq/dum)
550
551//      val = V_WhiteBeamDist(dum*5.3)*BroadPeakX(cw,qq/dum)
552        val = V_WhiteBeamDist_mid(dum*5.3)*func(cw,qq/dum)
553       
554//      val = V_WhiteBeamInterp(dum*5.3)*func(cw,qq/dum)
555
556        return (val)
557End
558
559Function V_intgrnd_triangle(cw,dum)
560        Wave cw
561        Variable dum            // the dummy of the integration
562
563        Variable val
564        NVAR qq = root:qq               //the q-value of the integration, not part of cw, so pass global
565        SVAR funcStr = root:gFunctionString
566        FUNCREF SANSModel_proto func = $funcStr
567
568//      val = (1-dum*5.3/8)*BroadPeakX(cw,qq/dum)       
569        val = (1-dum*5.3/8)*func(cw,qq/dum)
570
571//      val = V_WhiteBeamDist(dum*5.3)*BroadPeakX(cw,qq/dum)
572//      val = V_WhiteBeamDist(dum*5.3)*func(cw,qq/dum)
573       
574//      val = V_WhiteBeamInterp(dum*5.3)*func(cw,qq/dum)
575
576        return (val)
577End
578
579Function V_intgrnd_interp(cw,dum)
580        Wave cw
581        Variable dum            // the dummy of the integration
582
583        Variable val
584        NVAR qq = root:qq               //the q-value of the integration, not part of cw, so pass global
585        SVAR funcStr = root:gFunctionString
586        FUNCREF SANSModel_proto func = $funcStr
587
588//      val = (1-dum*5.3/8)*BroadPeakX(cw,qq/dum)       
589//      val = (1-dum*5.3/8)*func(cw,qq/dum)
590
591//      val = V_WhiteBeamDist(dum*5.3)*BroadPeakX(cw,qq/dum)
592//      val = V_WhiteBeamDist(dum*5.3)*func(cw,qq/dum)
593       
594        val = V_WhiteBeamInterp(dum*5.3)*func(cw,qq/dum)
595
596        return (val)
597End
598
599
600////////////////////////////
601
602// need a function to return the model function name
603// given the coefficient wave
604//
605// want the function NameX for use in the integration, not the AAO function
606//
607
608// from the name of the coefficient wave, get the function name
609// be sure that there is no "Smeared" at the beginning of the name
610// tag X to the end of the name string
611//
612// then the funcString must be passed in as a global to the built-in integration function.
613//
614Function/S V_getXFuncStrFromCoef(cw)
615        Wave cw
616       
617        String cwStr = NameOfWave(cw)
618        String outStr = "",extStr=""
619       
620//      String convStr = ReplaceString("_",cwStr,".")           // change the _ to .   
621//      extStr = ParseFilePath(4, convStr, ":", 0, 0)           // extracts the last .nnn, without the .
622
623// go through the list of coefKWStr pairs
624// look for the cwStr
625// take up to the = (that is the funcStr)
626// remove "Smeared" if needed   
627        SVAR coefList=root:Packages:NIST:coefKWStr
628
629        Variable ii,num
630        String item
631       
632        num=ItemsInList(coefList,";")
633        ii=0
634        do
635                item = StringFromList(ii, coefList, ";")
636               
637                if(strsearch(item,cwStr,0) != -1)               //match
638                        item = ReplaceString("=",item,".")              //replace the = with .
639                        outStr = ParseFilePath(3, item, ":", 0, 0)              // extract file name without extension
640                        outStr = ReplaceString("Smeared",outStr,"")             // replace "Smeared" with null, if it's there
641                        ii = num + 1
642                endif
643               
644                ii+=1
645        while(ii<num)
646       
647        return(outStr+"X")
648end
649
650//////////////////////////////////////////
651// generates dummy functions of the form:
652//
653//Function BroadPeakX_(cw,x)
654//      Wave cw
655//      Variable x
656//      return(BroadPeakX(cw,x))
657//End
658//
659// so that I can use the FUNCREF
660// which fails for some reason when I just use the XOP name?
661//
662//
663// not everything ending in X is a model function - trimmed list is in V_DummyFunctions.ipf
664//
665Function V_generateDummyFuncs()
666
667        String list = FunctionList("*X",";","KIND:4")
668        Variable ii,num
669        String item,str
670       
671        num=ItemsInList(list,";")
672
673        NewNotebook/N=Notebook1/F=0
674       
675       
676        for(ii=0;ii<num;ii+=1)
677                item = StringFromList(ii,list,";")
678                str = "\r"
679                str = "Function "+item+"_(cw,x)\r"
680                str += "\tWave cw\r"
681                str += "\tVariable x\r"
682                str += "\treturn("+item+"(cw,x))\r"
683                str += "End\r\r"
684               
685                //print str
686       
687                Notebook $"", text=str
688               
689        endfor
690        return(0)
691       
692End
Note: See TracBrowser for help on using the repository browser.