source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/BroadPeak_Pix_2D.ipf @ 956

Last change on this file since 956 was 956, checked in by srkline, 8 years ago

another iteration of the proposed Nexus file structure for SANS and VSANS

added 2D fit function to be able to fit the pixel beam center for an arc on the detector. this is still in a rudimentary form, but is somewhat functional.

File size: 7.4 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 PlotBroadPeak_Pix2D(xDim,yDim)                                             
30        Variable xDim=48, yDim=256
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,4, 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        FillPixTriplet(xwave_PeakPix2D, ywave_PeakPix2D,zwave_PeakPix2D,xDim,yDim)
42       
43       
44        Variable/G g_PeakPix2D=0
45        g_PeakPix2D := 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 := 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
71Function FillPixTriplet(xwave_PeakPix2D, ywave_PeakPix2D,zwave_PeakPix2D,xDim,yDim)
72        Wave xwave_PeakPix2D, ywave_PeakPix2D,zwave_PeakPix2D
73        Variable xDim,yDim
74               
75        Variable ii,jj 
76        ii=0
77        jj=0
78        do
79                do
80                        xwave_PeakPix2D[ii*yDim+ jj] = ii
81                        ywave_PeakPix2D[ii*yDim+ jj] = jj
82                        jj+=1
83                while(jj<yDim)
84                jj=0
85                ii+=1
86        while(ii<xDim)
87        return(0)
88End
89
90Function UpdatePix2Mat(Qx,Qy,inten,linMat,mat)
91        Wave Qx,Qy,inten,linMat,mat
92       
93        Variable xrows=DimSize(mat, 0 )                 
94        Variable yrows=DimSize(mat, 1 )                 
95               
96        String folderStr=GetWavesDataFolder(Qx,1)
97        NVAR/Z gIsLogScale=$(folderStr+"gIsLogScale")
98       
99//      linMat = inten[q*xrows+p]
100        linMat = inten[p*yrows+q]
101       
102        if(gIsLogScale)
103                mat = log(linMat)
104        else
105                mat = linMat
106        endif
107       
108        return(0)
109End
110
111//
112//  Fit function that is actually a wrapper to dispatch the calculation to N threads
113//
114// nthreads is 1 or an even number, typically 2
115// it doesn't matter if npt is odd. In this case, fractional point numbers are passed
116// and the wave indexing works just fine - I tested this with test waves of 7 and 8 points
117// and the points "2.5" and "3.5" evaluate correctly as 2 and 3
118//
119Function BroadPeak_Pix2D(cw,zw,xw,yw) : FitFunc
120        Wave cw,zw,xw,yw
121       
122#if exists("BroadPeak_Pix2DX")                  //to hide the function if XOP not installed
123        MultiThread zw = BroadPeak_Pix2DX(cw,xw,yw)
124#else
125        MultiThread zw = I_BroadPeak_Pix2D(cw,xw,yw)
126#endif
127       
128        return(0)
129End
130
131//threaded version of the function
132ThreadSafe Function BroadPeak_Pix2D_T(cw,zw,xw,yw,p1,p2)
133        WAVE cw,zw,xw,yw
134        Variable p1,p2
135       
136#if exists("BroadPeak_Pix2DX")                  //to hide the function if XOP not installed
137        zw[p1,p2]= BroadPeak_Pix2DX(cw,xw,yw)
138#else
139        zw[p1,p2]= I_BroadPeak_Pix2D(cw,xw,yw)
140#endif
141
142        return 0
143End
144
145//// technically, I'm passing a coefficient wave that's TOO LONG to the XOP
146//// BEWARE
147//ThreadSafe Function I_BroadPeak_Pix2D(w,x,y)
148//      Wave w
149//      Variable x,y
150//     
151//      Variable retVal,qval
152////    WAVE tmp = root:tmp_Pix2D
153////    tmp = w[p]
154//     
155//      Variable xCtr,yCtr
156//      xCtr = w[7]
157//      yCtr = w[8]
158//     
159////    qval = sqrt((x-xCtr)^2+(y-yCtr)^2)                      // use if the pixels are square
160//      qval = sqrt((x-xCtr)^2+(y-yCtr)^2/4)                    // use for LR panels where the y pixels are half the size of x
161////    qval = sqrt((x-xCtr)^2/4+(y-yCtr)^2)                    // use for TB panels where the y pixels are twice the size of x
162//
163//      if(qval< 0.001)
164//              retval = w[6]                   //bgd
165//      else           
166//              retval = BroadPeakX(w,qval)             //pass only what BroadPeak needs
167////            retval = BroadPeakX(tmp,qval)           //pass only what BroadPeak needs
168//      endif
169//             
170//      return(retVal)
171//End
172
173
174
175//
176// This is not an XOP, but is correct in what it is passing and speed seems to be just fine.
177//
178ThreadSafe Function I_BroadPeak_Pix2D(w,xw,yw)
179//ThreadSafe Function fBroadPeak_Pix2D(w,xw,yw)
180        Wave w
181        Variable xw,yw
182
183        // variables are:                                                       
184        //[0] Porod term scaling
185        //[1] Porod exponent
186        //[2] Lorentzian term scaling
187        //[3] Lorentzian screening length [A]
188        //[4] peak location [1/A]
189        //[5] Lorentzian exponent
190        //[6] background
191       
192        //[7] xSize
193        //[8] ySize
194       
195        //[9] xCtr
196        //[10] yCtr
197       
198        Variable aa,nn,cc,LL,Qzero,mm,bgd,xctr,yctr,xSize,ySize
199        aa = w[0]
200        nn = w[1]
201        cc = w[2]
202        LL=w[3]
203        Qzero=w[4]
204        mm=w[5]
205        bgd=w[6]
206        xSize = w[7]
207        ySize = w[8]
208        xCtr = w[9]
209        yCtr = w[10]
210       
211//      local variables
212        Variable inten, qval,ratio
213
214//      x is the q-value for the calculation
215//      qval = sqrt(xw^2+yw^2)
216
217// ASSUMPTION
218// TODO (change this)
219// base the scaling on the xSize
220
221        ratio = (xSize/ySize)^2
222        if(ratio > 1)
223        //      qval = sqrt((xw-xCtr)^2+(yw-yCtr)^2)                    // use if the pixels are square
224                qval = sqrt((xw-xCtr)^2+(yw-yCtr)^2/ratio)                      // use for LR panels where the y pixels are half the size of x 
225        else
226                qval = sqrt((xw-xCtr)^2*ratio+(yw-yCtr)^2)                      // use for TB panels where the y pixels are twice the size of x
227        endif
228
229//      qval = sqrt((xw-xCtr)^2+(yw-yCtr)^2)                    // use if the pixels are square
230//      qval = sqrt((xw-xCtr)^2+(yw-yCtr)^2/4)                  // use for LR panels where the y pixels are half the size of x
231//      qval = sqrt((xw-xCtr)^2/4+(yw-yCtr)^2)                  // use for TB panels where the y pixels are twice the size of x
232
233        if(qval<.001)
234                return(bgd)
235        endif   
236       
237//      do the calculation and return the function value
238       
239        inten = aa/(qval)^nn + cc/(1 + (abs(qval-Qzero)*LL)^mm) + bgd
240
241        Return (inten)
242
243
244End
245
246
247
248//non-threaded version of the function, necessary for the smearing calculation
249// -- the smearing calculation can only calculate (nord) points at a time.
250//
251ThreadSafe Function BroadPeak_Pix2D_noThread(cw,zw,xw,yw)
252        WAVE cw,zw, xw,yw
253       
254#if exists("BroadPeak_Pix2DX")                  //to hide the function if XOP not installed
255        zw = BroadPeak_Pix2DX(cw,xw,yw)
256#else
257        zw = I_BroadPeak_Pix2D(cw,xw,yw)
258#endif
259
260        return 0
261End
262
263
264
265
266
267
268
269
270
271
272
273
274
Note: See TracBrowser for help on using the repository browser.