source: sans/Dev/trunk/NCNR_User_Procedures/Analysis/VSANS/V_WB_GaussSpheres.ipf @ 1098

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

renamed white beam smearing models so they would be grouped together on the file list.

Re-worked the logic and flow of the averaging/plotting/saving steps of the reduction protocol so that it would flow cleanly and leave room for changes for the multitude of different collimation conditions. The averaging routines are now aware of the collimation conditions so that the appropriate resolution can be calculated. The collimation string is also written out to the averaged data file as element[9] of the protocol. The hope is that one could key on this collimation string to decide how to proceed with the analysis.

File size: 5.9 KB
Line 
1#pragma TextEncoding = "MacRoman"
2#pragma rtGlobals=3             // Use modern global access method and strict wave access.
3
4
5// turn the White Beam resolution smearing into a fitting function
6// so that the wavelength smeared function can then be smeared by
7// a gaussian resolution function that has geometry only.
8//
9// The geometry only resolution is generated by passing dl/l=0 to the resolution
10// calculation.
11//
12// This representation uses the "middle" of the distribution
13
14//
15//
16
17//#include "sphere_v40"
18// plots the form factor of  spheres with a Gaussian radius distribution
19//
20// also can plot the distribution itself, based on the current model parameters
21//
22// integration is currently done using 20-pt quadrature, but may benefit from
23//switching to an adaptive integration.
24//
25
26Proc PlotGaussSpheresWB(num,qmin,qmax)
27        Variable num=128,qmin=0.001,qmax=0.7
28        Prompt num "Enter number of data points for model: "
29        Prompt qmin "Enter minimum q-value (A^-1) for model: "
30        Prompt qmax "Enter maximum q-value (A^-1) for model: "
31       
32        Make/O/D/N=(num) xwave_pgsWB,ywave_pgsWB
33        xwave_pgsWB = alog( log(qmin) + x*((log(qmax)-log(qmin))/num) )
34        Make/O/D coef_pgsWB = {0.01,60,0.2,1e-6,3e-6,0.001}
35        make/O/T parameters_pgsWB = {"Volume Fraction (scale)","mean radius (A)","polydisp (sig/avg)","SLD sphere (A-2)","SLD solvent (A-2)","bkg (cm-1 sr-1)"}
36        Edit parameters_pgsWB,coef_pgsWB
37       
38        Variable/G root:g_pgsWB
39        g_pgsWB := GaussSpheresWB(coef_pgsWB,ywave_pgsWB,xwave_pgsWB)
40        Display ywave_pgsWB vs xwave_pgsWB
41        ModifyGraph log=1,marker=29,msize=2,mode=4
42        Label bottom "q (A\\S-1\\M)"
43        Label left "Intensity (cm\\S-1\\M)"
44        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
45       
46        AddModelToStrings("GaussSpheresWB","coef_pgsWB","parameters_pgsWB","pgsWB")
47End
48
49// - sets up a dependency to a wrapper, not the actual SmearedModelFunction
50Proc PlotSmearedGaussSpheresWB(str)                                                             
51        String str
52        Prompt str,"Pick the data folder containing the resolution you want",popup,getAList(4)
53       
54        // if any of the resolution waves are missing => abort
55        if(ResolutionWavesMissingDF(str))               //updated to NOT use global strings (in GaussUtils)
56                Abort
57        endif
58       
59        SetDataFolder $("root:"+str)
60       
61        // Setup parameter table for model function
62        Make/O/D smear_coef_pgsWB = {0.01,60,0.2,1e-6,3e-6,0.001}                                       
63        make/o/t smear_parameters_pgsWB = {"Volume Fraction (scale)","mean radius (A)","polydisp (sig/avg)","SLD sphere (A-2)","SLD solvent (A-2)","bkg (cm-1 sr-1)"}   
64        Edit smear_parameters_pgsWB,smear_coef_pgsWB                                   
65       
66        // output smeared intensity wave, dimensions are identical to experimental QSIG values
67        // make extra copy of experimental q-values for easy plotting
68        Duplicate/O $(str+"_q") smeared_pgsWB,smeared_qvals                             
69        SetScale d,0,0,"1/cm",smeared_pgsWB                                                     
70                                       
71        Variable/G gs_pgsWB=0
72        gs_pgsWB := fSmearedGaussSpheresWB(smear_coef_pgsWB,smeared_pgsWB,smeared_qvals)        //this wrapper fills the STRUCT
73       
74        Display smeared_pgsWB vs smeared_qvals                                                                 
75        ModifyGraph log=1,marker=29,msize=2,mode=4
76        Label bottom "q (A\\S-1\\M)"
77        Label left "Intensity (cm\\S-1\\M)"
78        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
79       
80        SetDataFolder root:
81        AddModelToStrings("SmearedGaussSpheresWB","smear_coef_pgsWB","smear_parameters_pgsWB","pgsWB")
82End
83       
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 GaussSpheresWB(cw,yw,xw) : FitFunc
91        Wave cw,yw,xw
92       
93#if exists("GaussSpheresX")
94//      MultiThread yw = GaussSpheresX(cw,xw)
95        yw = V_fGaussSpheresWB(cw,xw)
96
97#else
98//      yw = fGaussSpheresWB(cw,xw)
99        yw = 1
100#endif
101        return(0)
102End
103
104Function V_fGaussSpheresWB(w,xx) : FitFunc
105        wave w
106        variable xx
107       
108        Variable scale,rad,pd,sig,rho,rhos,bkg,delrho,inten,loLim,upLim
109       
110        //the coefficient values
111//      scale=w[0]
112//      rad=w[1]
113//      pd=w[2]
114//      sig=pd*rad
115//      rho=w[3]
116//      rhos=w[4]
117//      delrho=rho-rhos
118//      bkg=w[5]
119       
120       
121        // define limits based on lo/mean, hi/mean of the wavelength distribution
122        // using the empirical definition, "middle" of the peaks
123        loLim = 3.37/5.3
124        upLim = 8.37/5.3
125       
126        inten = V_IntegrGaussSphereWB_mid(w,loLim,upLim,xx)
127
128// why do I need this? Is this because this is defined as the mean of the distribution
129//  and is needed to normalize the integral? verify this on paper.     
130        inten *= 5.3
131
132// normalize the integral       
133        inten /= 19933          // "middle"  of peaks
134
135// additional normalization???
136        inten /= 1.05           //
137       
138        Return(inten)
139End
140
141
142// the trick here is that declaring the last qVal wave as a variable
143// since this is implicitly called N times in the wave assignment of the answer wave
144Function V_IntegrGaussSphereWB_mid(cw,loLim,upLim,qVal)
145        Wave cw
146        Variable loLim,upLim
147        Variable qVal
148       
149        Variable/G root:qq = qval
150        Variable ans
151       
152//      ans = Integrate1D(V_intgrnd_top,lolim,uplim,2,0,cw)             //adaptive quadrature
153        ans = Integrate1D(V_integrand_pgsWB,lolim,uplim,1,0,cw)         // Romberg integration
154       
155        return ans
156end
157
158Function V_integrand_pgsWB(cw,dum)
159        Wave cw
160        Variable dum            // the dummy of the integration
161
162        Variable val
163        NVAR qq = root:qq               //the q-value of the integration, not part of cw, so pass global
164//      SVAR funcStr = root:gFunctionString
165//      FUNCREF SANSModel_proto func = $funcStr
166
167        val = V_WhiteBeamDist_mid(dum*5.3)*GaussSpheresX(cw,qq/dum)
168       
169        return (val)
170End
171
172//wrapper to calculate the smeared model as an AAO-Struct
173// fills the struct and calls the ususal function with the STRUCT parameter
174//
175// used only for the dependency, not for fitting
176//
177Function fSmearedGaussSpheresWB(coefW,yW,xW)
178        Wave coefW,yW,xW
179       
180        String str = getWavesDataFolder(yW,0)
181        String DF="root:"+str+":"
182       
183        WAVE resW = $(DF+str+"_res")
184       
185        STRUCT ResSmearAAOStruct fs
186        WAVE fs.coefW = coefW   
187        WAVE fs.yW = yW
188        WAVE fs.xW = xW
189        WAVE fs.resW = resW
190       
191        Variable err
192        err = SmearedGaussSpheresWB(fs)
193       
194        return (0)
195End
196
197// this is all there is to the smeared calculation!
198Function SmearedGaussSpheresWB(s) :FitFunc
199        Struct ResSmearAAOStruct &s
200
201//      the name of your unsmeared model (AAO) is the first argument
202        Smear_Model_20(GaussSpheresWB,s.coefW,s.xW,s.yW,s.resW)
203
204        return(0)
205End
206
207
Note: See TracBrowser for help on using the repository browser.