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

Last change on this file since 1129 was 1119, checked in by srkline, 4 years ago

added procedures to output QxQy_ASCII data. Each panel is output into its own file. Output format is the same as for 2D SANS data, including the 2D resolution function. However, reading the data back in with the SANS analysis macros currently does not redimension the data back to the matrix correctly as it assumes a square detector.

I will need to add the X and Y dimensions of each panel into the header, and then make use of these values when they are read in. Also, writing the QxQy? data is quick for the M and F panels ( 0.8 s) but is extremely slow for the back, High-res panel (120 s) since there are 1.1.E6 points there vs. 6144 pts. I'll need to find a way to speed this operation up.

File size: 7.5 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//
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//
226// base the scaling on the xSize
227
228// *** NOTE ***
229// "qval" here is NOT q
230// qval is a real space distance in the units of PIXELS. Not mm, not q, PIXELS
231// Don't put any meaning in the fitted values, it's simply a functional shape
232
233        ratio = (xSize/ySize)^2
234        if(ratio > 1)
235        //      qval = sqrt((xw-xCtr)^2+(yw-yCtr)^2)                    // use if the pixels are square
236                qval = sqrt((xw-xCtr)^2+((yw-yCtr)^2)/ratio)                    // use for LR panels where the y pixels are half the size of x 
237        else
238                qval = sqrt(((xw-xCtr)^2)*ratio+(yw-yCtr)^2)                    // use for TB panels where the y pixels are twice the size of x
239        endif
240
241        if(qval<.001)
242                return(bgd)
243        endif   
244       
245//      do the calculation and return the function value
246       
247        inten = aa/(qval)^nn + cc/(1 + (abs(qval-Qzero)*LL)^mm) + bgd
248
249        Return (inten)
250
251
252End
253
254
255
256//non-threaded version of the function, necessary for the smearing calculation
257// -- the smearing calculation can only calculate (nord) points at a time.
258//
259ThreadSafe Function V_BroadPeak_Pix2D_noThread(cw,zw,xw,yw)
260        WAVE cw,zw, xw,yw
261       
262#if exists("BroadPeak_Pix2DX")                  //to hide the function if XOP not installed
263        zw = BroadPeak_Pix2DX(cw,xw,yw)
264#else
265        zw = V_I_BroadPeak_Pix2D(cw,xw,yw)
266#endif
267
268        return 0
269End
270
271
272
273
274
275
276
277
278
279
280
281
282
Note: See TracBrowser for help on using the repository browser.