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 | // |
---|
29 | Proc 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") |
---|
69 | End |
---|
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 | // |
---|
78 | Function 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) |
---|
95 | End |
---|
96 | |
---|
97 | Function 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) |
---|
116 | End |
---|
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 | // |
---|
126 | Function 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) |
---|
136 | End |
---|
137 | |
---|
138 | //threaded version of the function |
---|
139 | ThreadSafe 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 |
---|
150 | End |
---|
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 | // |
---|
185 | ThreadSafe 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 | // *** 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 | |
---|
252 | End |
---|
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 | // |
---|
259 | ThreadSafe 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 |
---|
269 | End |
---|
270 | |
---|
271 | |
---|
272 | |
---|
273 | |
---|
274 | |
---|
275 | |
---|
276 | |
---|
277 | |
---|
278 | |
---|
279 | |
---|
280 | |
---|
281 | |
---|
282 | |
---|