source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_BroadPeak_Pix_2D.ipf @ 1041

Last change on this file since 1041 was 1041, checked in by srkline, 6 years ago

Added the angle dependent transmission correction to the data correction in the raw_to_work step, in 2D

added a testing file that can generate fake event data, read, write, and decode it. Read is based on GBLoadWave. Hoepfully I'll not need to write an XOP. manipulation of the 64 bit words are done with simple bit shifts and logic.

also added are a number of error checking routines to improve behavior when wave, folders, etc. are missing.

File size: 7.3 KB
Line 
1#pragma rtGlobals=3             // Use modern global access method and strict wave access.
2#pragma IgorVersion=6.1
3
4
5// technically, I'm passing a coefficient wave that's TOO LONG to the XOP
6// BEWARE: see
7//ThreadSafe Function I_BroadPeak_Pix2D(w,x,y)
8
9// ALSO -- the pixels are not square in general, so this will add more complications...
10//      qval = sqrt((x-xCtr)^2+(y-yCtr)^2)                      // use if the pixels are square
11//      qval = sqrt((x-xCtr)^2+(y-yCtr)^2/4)                    // use for LR panels where the y pixels are half the size of x
12//      qval = sqrt((x-xCtr)^2/4+(y-yCtr)^2)                    // use for TB panels where the y pixels are twice the size of x
13
14//
15//
16// WaveStats/Q data_FL
17// coef_peakPix2d[2] = V_max
18// coef_peakPix2d[0] = 1
19// then set the xy center to something somewhat close (could do this based on FL, etc.)
20// then set the peak position somewhat close (how to do this??)
21//
22// FuncFitMD/H="11000111100"/NTHR=0 BroadPeak_Pix2D coef_PeakPix2D  data_FT /D
23//
24
25//
26// the calculation is done as for the QxQy data set:
27// three waves XYZ, then converted to a matrix
28//
29Proc V_PlotBroadPeak_Pix2D(xDim,yDim)                                           
30        Variable xDim=48, yDim=128
31        Prompt xDim "Enter X dimension: "
32        Prompt yDim "Enter Y dimension: "
33               
34        Make/O/D coef_PeakPix2D = {10, 3, 10, 0.3, 10, 2, 0.1, 8, 8, 100, 100}
35//      Make/O/D tmp_Pix2D =    {10, 3, 10, 0.3, 10, 2, 0.1}            //without the pixel ctrs                                       
36        make/o/t parameters_PeakPix2D = {"Porod Scale", "Porod Exponent","Lorentzian Scale","Lor Screening Length","Peak position","Lorentzian Exponent","Bgd [1/cm]", "xPix size (mm)","yPix size (mm)", "xCtr (pixels)", "yCtr (pixels)"}             
37        Edit parameters_PeakPix2D,coef_PeakPix2D                               
38       
39        // generate the triplet representation
40        Make/O/D/N=(xDim*yDim) xwave_PeakPix2D, ywave_PeakPix2D,zwave_PeakPix2D
41        V_FillPixTriplet(xwave_PeakPix2D, ywave_PeakPix2D,zwave_PeakPix2D,xDim,yDim)
42       
43       
44        Variable/G g_PeakPix2D=0
45        g_PeakPix2D := V_BroadPeak_Pix2D(coef_PeakPix2D,zwave_PeakPix2D,xwave_PeakPix2D,ywave_PeakPix2D)        //AAO 2D calculation
46       
47        Display ywave_PeakPix2D vs xwave_PeakPix2D
48        modifygraph log=0
49        ModifyGraph mode=3,marker=16,zColor(ywave_PeakPix2D)={zwave_PeakPix2D,*,*,YellowHot,0}
50        ModifyGraph standoff=0
51        ModifyGraph width={Plan,1,bottom,left}
52        ModifyGraph lowTrip=0.001
53        Label bottom "X pixels"
54        Label left "Y pixels"
55        AutoPositionWindow/M=1/R=$(WinName(0,1)) $WinName(0,2)
56       
57        // generate the matrix representation
58        Make/O/D/N=(xDim,yDim) PeakPix2D_mat            // use the point scaling of the matrix (=pixels)
59        Duplicate/O $"PeakPix2D_mat",$"PeakPix2D_lin"           //keep a linear-scaled version of the data
60        // _mat is for display, _lin is the real calculation
61
62        // not a function evaluation - this simply keeps the matrix for display in sync with the triplet calculation
63        Variable/G g_PeakPix2Dmat=0
64        g_PeakPix2Dmat := V_UpdatePix2Mat(xwave_PeakPix2D,ywave_PeakPix2D,zwave_PeakPix2D,PeakPix2D_lin,PeakPix2D_mat)
65       
66       
67        SetDataFolder root:
68//      AddModelToStrings("BroadPeak_Pix2D","coef_PeakPix2D","parameters_PeakPix2D","PeakPix2D")
69End
70
71
72//
73// this sets the x and y waves of the triplet to be the pixel numbers
74//
75// TODO:
76// -- this will need to be changed if I want to fit based on real-space mm
77//
78Function V_FillPixTriplet(xwave_PeakPix2D, ywave_PeakPix2D,zwave_PeakPix2D,xDim,yDim)
79        Wave xwave_PeakPix2D, ywave_PeakPix2D,zwave_PeakPix2D
80        Variable xDim,yDim
81               
82        Variable ii,jj 
83        ii=0
84        jj=0
85        do
86                do
87                        xwave_PeakPix2D[ii*yDim+ jj] = ii
88                        ywave_PeakPix2D[ii*yDim+ jj] = jj
89                        jj+=1
90                while(jj<yDim)
91                jj=0
92                ii+=1
93        while(ii<xDim)
94        return(0)
95End
96
97Function V_UpdatePix2Mat(Qx,Qy,inten,linMat,mat)
98        Wave Qx,Qy,inten,linMat,mat
99       
100        Variable xrows=DimSize(mat, 0 )                 
101        Variable yrows=DimSize(mat, 1 )                 
102               
103        String folderStr=GetWavesDataFolder(Qx,1)
104        NVAR/Z gIsLogScale=$(folderStr+"gIsLogScale")
105       
106//      linMat = inten[q*xrows+p]
107        linMat = inten[p*yrows+q]
108       
109        if(gIsLogScale)
110                mat = log(linMat)
111        else
112                mat = linMat
113        endif
114       
115        return(0)
116End
117
118//
119//  Fit function that is actually a wrapper to dispatch the calculation to N threads
120//
121// nthreads is 1 or an even number, typically 2
122// it doesn't matter if npt is odd. In this case, fractional point numbers are passed
123// and the wave indexing works just fine - I tested this with test waves of 7 and 8 points
124// and the points "2.5" and "3.5" evaluate correctly as 2 and 3
125//
126Function V_BroadPeak_Pix2D(cw,zw,xw,yw) : FitFunc
127        Wave cw,zw,xw,yw
128       
129#if exists("BroadPeak_Pix2DX")                  //to hide the function if XOP not installed
130        MultiThread zw = BroadPeak_Pix2DX(cw,xw,yw)
131#else
132        MultiThread zw = V_I_BroadPeak_Pix2D(cw,xw,yw)
133#endif
134       
135        return(0)
136End
137
138//threaded version of the function
139ThreadSafe Function V_BroadPeak_Pix2D_T(cw,zw,xw,yw,p1,p2)
140        WAVE cw,zw,xw,yw
141        Variable p1,p2
142       
143#if exists("BroadPeak_Pix2DX")                  //to hide the function if XOP not installed
144        zw[p1,p2]= BroadPeak_Pix2DX(cw,xw,yw)
145#else
146        zw[p1,p2]= V_I_BroadPeak_Pix2D(cw,xw,yw)
147#endif
148
149        return 0
150End
151
152//// technically, I'm passing a coefficient wave that's TOO LONG to the XOP
153//// BEWARE
154//ThreadSafe Function I_BroadPeak_Pix2D(w,x,y)
155//      Wave w
156//      Variable x,y
157//     
158//      Variable retVal,qval
159////    WAVE tmp = root:tmp_Pix2D
160////    tmp = w[p]
161//     
162//      Variable xCtr,yCtr
163//      xCtr = w[7]
164//      yCtr = w[8]
165//     
166////    qval = sqrt((x-xCtr)^2+(y-yCtr)^2)                      // use if the pixels are square
167//      qval = sqrt((x-xCtr)^2+(y-yCtr)^2/4)                    // use for LR panels where the y pixels are half the size of x
168////    qval = sqrt((x-xCtr)^2/4+(y-yCtr)^2)                    // use for TB panels where the y pixels are twice the size of x
169//
170//      if(qval< 0.001)
171//              retval = w[6]                   //bgd
172//      else           
173//              retval = BroadPeakX(w,qval)             //pass only what BroadPeak needs
174////            retval = BroadPeakX(tmp,qval)           //pass only what BroadPeak needs
175//      endif
176//             
177//      return(retVal)
178//End
179
180
181
182//
183// This is not an XOP, but is correct in what it is passing and speed seems to be just fine.
184//
185ThreadSafe Function V_I_BroadPeak_Pix2D(w,xw,yw)
186//ThreadSafe Function fBroadPeak_Pix2D(w,xw,yw)
187        Wave w
188        Variable xw,yw
189
190        // variables are:                                                       
191        //[0] Porod term scaling
192        //[1] Porod exponent
193        //[2] Lorentzian term scaling
194        //[3] Lorentzian screening length [A]
195        //[4] peak location [1/A]
196        //[5] Lorentzian exponent
197        //[6] background
198       
199        //[7] xSize
200        //[8] ySize
201       
202        //[9] xCtr
203        //[10] yCtr
204       
205        Variable aa,nn,cc,LL,Qzero,mm,bgd,xctr,yctr,xSize,ySize
206        aa = w[0]
207        nn = w[1]
208        cc = w[2]
209        LL=w[3]
210        Qzero=w[4]
211        mm=w[5]
212        bgd=w[6]
213        xSize = w[7]
214        ySize = w[8]
215        xCtr = w[9]
216        yCtr = w[10]
217       
218//      local variables
219        Variable inten, qval,ratio
220
221//      x is the q-value for the calculation
222//      qval = sqrt(xw^2+yw^2)
223
224// ASSUMPTION
225// TODO (change this)
226// base the scaling on the xSize
227
228        ratio = (xSize/ySize)^2
229        if(ratio > 1)
230        //      qval = sqrt((xw-xCtr)^2+(yw-yCtr)^2)                    // use if the pixels are square
231                qval = sqrt((xw-xCtr)^2+(yw-yCtr)^2/ratio)                      // use for LR panels where the y pixels are half the size of x 
232        else
233                qval = sqrt((xw-xCtr)^2*ratio+(yw-yCtr)^2)                      // use for TB panels where the y pixels are twice the size of x
234        endif
235
236        if(qval<.001)
237                return(bgd)
238        endif   
239       
240//      do the calculation and return the function value
241       
242        inten = aa/(qval)^nn + cc/(1 + (abs(qval-Qzero)*LL)^mm) + bgd
243
244        Return (inten)
245
246
247End
248
249
250
251//non-threaded version of the function, necessary for the smearing calculation
252// -- the smearing calculation can only calculate (nord) points at a time.
253//
254ThreadSafe Function V_BroadPeak_Pix2D_noThread(cw,zw,xw,yw)
255        WAVE cw,zw, xw,yw
256       
257#if exists("BroadPeak_Pix2DX")                  //to hide the function if XOP not installed
258        zw = BroadPeak_Pix2DX(cw,xw,yw)
259#else
260        zw = V_I_BroadPeak_Pix2D(cw,xw,yw)
261#endif
262
263        return 0
264End
265
266
267
268
269
270
271
272
273
274
275
276
277
Note: See TracBrowser for help on using the repository browser.