1 | #pragma rtGlobals=1 // Use modern global access method. |
---|
2 | #pragma IgorVersion=6.1 |
---|
3 | |
---|
4 | |
---|
5 | |
---|
6 | // there are utlity functions here (commented out) for plotting the 2D resolution function |
---|
7 | // at any point on the detector |
---|
8 | |
---|
9 | |
---|
10 | // I can only calculate (nord) points in an AAO fashion, so threading the function calculation |
---|
11 | // doesn't help. Even if I can rewrite this to calculate nord*nord AAO, that will typically be 10*10=100 |
---|
12 | // and that is not enough points to thread with any benefit |
---|
13 | |
---|
14 | // but the calculation of each q value is independent, so I can split the total number of q-values between processors |
---|
15 | // |
---|
16 | // -- but I must be careful to pass this a function that is not already threaded! |
---|
17 | // --- BUT - I can't pass either a function reference, OR a structure to a thread! |
---|
18 | // ---- So sadly, each 2D resolution calculation must be threaded by hand. |
---|
19 | // |
---|
20 | // See the Sphere_2D function for an example of how to thread the smearing |
---|
21 | // |
---|
22 | // |
---|
23 | // SRK May 2010 |
---|
24 | |
---|
25 | // |
---|
26 | // NOTE: there is a definition of FindTheta = -1* FindPhi() that is a duplicate of what is in RawWindowHook.ipf |
---|
27 | // and should eventually be in a common location for analysis and reduction packages |
---|
28 | // |
---|
29 | |
---|
30 | |
---|
31 | //// this is the completely generic 2D smearing, not threaded, but takes FUNCREF and STRUCT parameters |
---|
32 | // |
---|
33 | // uses resolution ellipse defined perpendicular "Y" and parallel "X", |
---|
34 | // rotating the ellipse into its proper orientaiton based on qx,qy |
---|
35 | // |
---|
36 | // 5 gauss points is not enough - it gives artifacts as a function of phi |
---|
37 | // 10 gauss points is minimally sufficient |
---|
38 | // 20 gauss points are needed if lots of oscillations (just like in 1D) |
---|
39 | // even more may be necessary for highly peaked functions |
---|
40 | // |
---|
41 | // |
---|
42 | Function Smear_2DModel_PP(fcn,s,nord) |
---|
43 | FUNCREF SANS_2D_ModelAAO_proto fcn |
---|
44 | Struct ResSmear_2D_AAOStruct &s |
---|
45 | Variable nord |
---|
46 | |
---|
47 | String weightStr,zStr |
---|
48 | |
---|
49 | Variable ii,jj,kk,num |
---|
50 | Variable qx,qy,qz,qval,fs |
---|
51 | Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut |
---|
52 | |
---|
53 | Variable a,b,c,normFactor,phi,theta,maxSig,numStdDev=3 |
---|
54 | |
---|
55 | switch(nord) |
---|
56 | case 5: |
---|
57 | weightStr="gauss5wt" |
---|
58 | zStr="gauss5z" |
---|
59 | if (WaveExists($weightStr) == 0) |
---|
60 | Make/O/D/N=(nord) $weightStr,$zStr |
---|
61 | Make5GaussPoints($weightStr,$zStr) |
---|
62 | endif |
---|
63 | break |
---|
64 | case 10: |
---|
65 | weightStr="gauss10wt" |
---|
66 | zStr="gauss10z" |
---|
67 | if (WaveExists($weightStr) == 0) |
---|
68 | Make/O/D/N=(nord) $weightStr,$zStr |
---|
69 | Make10GaussPoints($weightStr,$zStr) |
---|
70 | endif |
---|
71 | break |
---|
72 | case 20: |
---|
73 | weightStr="gauss20wt" |
---|
74 | zStr="gauss20z" |
---|
75 | if (WaveExists($weightStr) == 0) |
---|
76 | Make/O/D/N=(nord) $weightStr,$zStr |
---|
77 | Make20GaussPoints($weightStr,$zStr) |
---|
78 | endif |
---|
79 | break |
---|
80 | default: |
---|
81 | Abort "Smear_2DModel_PP called with invalid nord value" |
---|
82 | endswitch |
---|
83 | |
---|
84 | Wave/Z wt = $weightStr |
---|
85 | Wave/Z xi = $zStr |
---|
86 | |
---|
87 | /// keep these waves local |
---|
88 | // Make/O/D/N=1 yPtW |
---|
89 | Make/O/D/N=(nord) fcnRet,xptW,res_tot,yptW |
---|
90 | |
---|
91 | // now just loop over the points as specified |
---|
92 | num=numpnts(s.xw[0]) |
---|
93 | |
---|
94 | answer=0 |
---|
95 | |
---|
96 | Variable spl,spp,apl,app,bpl,bpp,phi_pt,qpl_pt |
---|
97 | Variable qperp_pt,phi_prime,q_prime |
---|
98 | |
---|
99 | Variable t1=StopMSTimer(-2) |
---|
100 | |
---|
101 | //loop over q-values |
---|
102 | for(ii=0;ii<num;ii+=1) |
---|
103 | |
---|
104 | // if(mod(ii, 1000 ) == 0) |
---|
105 | // Print "ii= ",ii |
---|
106 | // endif |
---|
107 | |
---|
108 | qx = s.xw[0][ii] |
---|
109 | qy = s.xw[1][ii] |
---|
110 | qz = s.qz[ii] |
---|
111 | qval = sqrt(qx^2+qy^2+qz^2) |
---|
112 | spl = s.sQpl[ii] |
---|
113 | spp = s.sQpp[ii] |
---|
114 | fs = s.fs[ii] |
---|
115 | |
---|
116 | normFactor = 2*pi*spl*spp |
---|
117 | |
---|
118 | phi = -1*FindTheta(qx,qy) //Findtheta is an exact duplicate of FindPhi() * -1 |
---|
119 | |
---|
120 | apl = -numStdDev*spl + qval //parallel = q integration limits |
---|
121 | bpl = numStdDev*spl + qval |
---|
122 | /// app = -numStdDev*spp + phi //perpendicular = phi integration limits (WRONG) |
---|
123 | /// bpp = numStdDev*spp + phi |
---|
124 | app = -numStdDev*spp + 0 //q_perp = 0 |
---|
125 | bpp = numStdDev*spp + 0 |
---|
126 | |
---|
127 | //make sure the limits are reasonable. |
---|
128 | if(apl < 0) |
---|
129 | apl = 0 |
---|
130 | endif |
---|
131 | // do I need to specially handle limits when phi ~ 0? |
---|
132 | |
---|
133 | |
---|
134 | sumOut = 0 |
---|
135 | for(jj=0;jj<nord;jj+=1) // call phi the "outer' |
---|
136 | /// phi_pt = (xi[jj]*(bpp-app)+app+bpp)/2 |
---|
137 | qperp_pt = (xi[jj]*(bpp-app)+app+bpp)/2 //this is now q_perp |
---|
138 | |
---|
139 | sumIn=0 |
---|
140 | for(kk=0;kk<nord;kk+=1) //at phi, integrate over Qpl |
---|
141 | |
---|
142 | qpl_pt = (xi[kk]*(bpl-apl)+apl+bpl)/2 |
---|
143 | |
---|
144 | /// FindQxQy(qpl_pt,phi_pt,qx_pt,qy_pt) //find the corresponding QxQy to the Q,phi |
---|
145 | |
---|
146 | // find QxQy given Qpl and Qperp on the grid |
---|
147 | // |
---|
148 | q_prime = sqrt(qpl_pt^2+qperp_pt^2) |
---|
149 | phi_prime = phi + qperp_pt/qpl_pt |
---|
150 | FindQxQy(q_prime,phi_prime,qx_pt,qy_pt) |
---|
151 | |
---|
152 | yPtw[kk] = qy_pt //phi is the same in this loop, but qy is not |
---|
153 | xPtW[kk] = qx_pt //qx is different here too, as we're varying Qpl |
---|
154 | |
---|
155 | res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (qperp_pt)^2/spp/spp ) ) |
---|
156 | /// res_tot[kk] = exp(-0.5*( (qpl_pt-qval)^2/spl/spl + (phi_pt-phi)^2/spp/spp ) ) |
---|
157 | res_tot[kk] /= normFactor |
---|
158 | // res_tot[kk] *= fs |
---|
159 | |
---|
160 | endfor |
---|
161 | |
---|
162 | fcn(s.coefW,fcnRet,xptw,yptw) //calculate nord pts at a time |
---|
163 | |
---|
164 | //sumIn += wt[jj]*wt[kk]*res_tot*fcnRet[0] |
---|
165 | fcnRet *= wt[jj]*wt*res_tot |
---|
166 | // |
---|
167 | answer += (bpl-apl)/2.0*sum(fcnRet) // |
---|
168 | endfor |
---|
169 | |
---|
170 | answer *= (bpp-app)/2.0 |
---|
171 | s.zw[ii] = answer |
---|
172 | endfor |
---|
173 | |
---|
174 | Variable elap = (StopMSTimer(-2) - t1)/1e6 |
---|
175 | Print "elapsed time = ",elap |
---|
176 | |
---|
177 | return(0) |
---|
178 | end |
---|
179 | |
---|
180 | |
---|
181 | // this is generic, but I need to declare the Cylinder2D function threadsafe |
---|
182 | // and this calculation is significantly slower than the manually threaded calculation |
---|
183 | // if the function is fast to calculate. Once the function has polydispersity on 2 or more parameters |
---|
184 | // then this AAO calculation and the manual threading are both painfully slow, and more similar in execution time |
---|
185 | // |
---|
186 | // For 128x128 data, and 10x10 smearing, using 25 points for each polydispersity (using DANSE code) |
---|
187 | // Successively making more things polydisperse gives the following timing (in seconds): |
---|
188 | // |
---|
189 | // monodisp sigR sigTheta sigPhi |
---|
190 | // manual THR 1.1 3.9 74 1844 |
---|
191 | // this AAO 8.6 11.4 104 1930 |
---|
192 | // |
---|
193 | // and using 5 points for each polydispersity: (I see no visual difference) |
---|
194 | // manual THR 1.1 1.6 3.9 16.1 |
---|
195 | // |
---|
196 | // so clearly -- use 5 points for the polydispersities, unless there's a good reason not to - and |
---|
197 | // certainly for a survey, it's the way to go. |
---|
198 | // |
---|
199 | Function Smear_2DModel_PP_AAO(fcn,s,nord) |
---|
200 | FUNCREF SANS_2D_ModelAAO_proto fcn |
---|
201 | Struct ResSmear_2D_AAOStruct &s |
---|
202 | Variable nord |
---|
203 | |
---|
204 | String weightStr,zStr |
---|
205 | |
---|
206 | Variable ii,jj,kk,num |
---|
207 | Variable qx,qy,qz,qval,fs |
---|
208 | Variable qy_pt,qx_pt,res_x,res_y,answer,sumIn,sumOut |
---|
209 | |
---|
210 | Variable a,b,c,normFactor,phi,theta,maxSig,numStdDev=3 |
---|
211 | |
---|
212 | switch(nord) |
---|
213 | case 5: |
---|
214 | weightStr="gauss5wt" |
---|
215 | zStr="gauss5z" |
---|
216 | if (WaveExists($weightStr) == 0) |
---|
217 | Make/O/D/N=(nord) $weightStr,$zStr |
---|
218 | Make5GaussPoints($weightStr,$zStr) |
---|
219 | endif |
---|
220 | break |
---|
221 | case 10: |
---|
222 | weightStr="gauss10wt" |
---|
223 | zStr="gauss10z" |
---|
224 | if (WaveExists($weightStr) == 0) |
---|
225 | Make/O/D/N=(nord) $weightStr,$zStr |
---|
226 | Make10GaussPoints($weightStr,$zStr) |
---|
227 | endif |
---|
228 | break |
---|
229 | case 20: |
---|
230 | weightStr="gauss20wt" |
---|
231 | zStr="gauss20z" |
---|
232 | if (WaveExists($weightStr) == 0) |
---|
233 | Make/O/D/N=(nord) $weightStr,$zStr |
---|
234 | Make20GaussPoints($weightStr,$zStr) |
---|
235 | endif |
---|
236 | break |
---|
237 | default: |
---|
238 | Abort "Smear_2DModel_PP called with invalid nord value" |
---|
239 | endswitch |
---|
240 | |
---|
241 | Wave/Z wt = $weightStr |
---|
242 | Wave/Z xi = $zStr |
---|
243 | |
---|
244 | /// keep these waves local |
---|
245 | // Make/O/D/N=1 yPtW |
---|
246 | Make/O/D/N=(nord*nord) fcnRet,xptW,res_tot,yptW,wts |
---|
247 | Make/O/D/N=(nord) phi_pt,qpl_pt,qperp_pt |
---|
248 | |
---|
249 | // now just loop over the points as specified |
---|
250 | num=numpnts(s.xw[0]) |
---|
251 | |
---|
252 | answer=0 |
---|
253 | |
---|
254 | Variable spl,spp,apl,app,bpl,bpp |
---|
255 | Variable phi_prime,q_prime |
---|
256 | |
---|
257 | Variable t1=StopMSTimer(-2) |
---|
258 | |
---|
259 | //loop over q-values |
---|
260 | for(ii=0;ii<num;ii+=1) |
---|
261 | |
---|
262 | // if(mod(ii, 1000 ) == 0) |
---|
263 | // Print "ii= ",ii |
---|
264 | // endif |
---|
265 | |
---|
266 | qx = s.xw[0][ii] |
---|
267 | qy = s.xw[1][ii] |
---|
268 | qz = s.qz[ii] |
---|
269 | qval = sqrt(qx^2+qy^2+qz^2) |
---|
270 | spl = s.sQpl[ii] |
---|
271 | spp = s.sQpp[ii] |
---|
272 | fs = s.fs[ii] |
---|
273 | |
---|
274 | normFactor = 2*pi*spl*spp |
---|
275 | |
---|
276 | phi = -1*FindTheta(qx,qy) //Findtheta is an exact duplicate of FindPhi() * -1 |
---|
277 | |
---|
278 | apl = -numStdDev*spl + qval //parallel = q integration limits |
---|
279 | bpl = numStdDev*spl + qval |
---|
280 | /// app = -numStdDev*spp + phi //perpendicular = phi integration limits (WRONG) |
---|
281 | /// bpp = numStdDev*spp + phi |
---|
282 | app = -numStdDev*spp + 0 //q_perp = 0 |
---|
283 | bpp = numStdDev*spp + 0 |
---|
284 | |
---|
285 | //make sure the limits are reasonable. |
---|
286 | if(apl < 0) |
---|
287 | apl = 0 |
---|
288 | endif |
---|
289 | // do I need to specially handle limits when phi ~ 0? |
---|
290 | |
---|
291 | |
---|
292 | // sumOut = 0 |
---|
293 | for(jj=0;jj<nord;jj+=1) // call phi the "outer' |
---|
294 | // phi_pt[jj] = (xi[jj]*(bpp-app)+app+bpp)/2 |
---|
295 | qperp_pt[jj] = (xi[jj]*(bpp-app)+app+bpp)/2 //this is now q_perp |
---|
296 | |
---|
297 | // sumIn=0 |
---|
298 | for(kk=0;kk<nord;kk+=1) //at phi, integrate over Qpl |
---|
299 | |
---|
300 | qpl_pt[kk] = (xi[kk]*(bpl-apl)+apl+bpl)/2 |
---|
301 | |
---|
302 | /// FindQxQy(qpl_pt[kk],phi_pt[jj],qx_pt,qy_pt) //find the corresponding QxQy to the Q,phi |
---|
303 | |
---|
304 | // find QxQy given Qpl and Qperp on the grid |
---|
305 | // |
---|
306 | q_prime = sqrt(qpl_pt[kk]^2+qperp_pt[jj]^2) |
---|
307 | phi_prime = phi + qperp_pt[jj]/qpl_pt[kk] |
---|
308 | FindQxQy(q_prime,phi_prime,qx_pt,qy_pt) |
---|
309 | |
---|
310 | yPtw[nord*jj+kk] = qy_pt //phi is the same in this loop, but qy is not |
---|
311 | xPtW[nord*jj+kk] = qx_pt //qx is different here too, as we're varying Qpl |
---|
312 | |
---|
313 | res_tot[nord*jj+kk] = exp(-0.5*( (qpl_pt[kk]-qval)^2/spl/spl + (qperp_pt[jj])^2/spp/spp ) ) |
---|
314 | // res_tot[nord*jj+kk] = exp(-0.5*( (qpl_pt[kk]-qval)^2/spl/spl + (phi_pt[jj]-phi)^2/spp/spp ) ) |
---|
315 | res_tot[nord*jj+kk] /= normFactor |
---|
316 | // res_tot[kk] *= fs |
---|
317 | |
---|
318 | //weighting |
---|
319 | wts[nord*jj+kk] = wt[jj]*wt[kk] |
---|
320 | endfor |
---|
321 | |
---|
322 | endfor |
---|
323 | |
---|
324 | fcn(s.coefW,fcnRet,xptw,yptw) //calculate nord*nord pts at a time |
---|
325 | |
---|
326 | fcnRet *= wts*res_tot |
---|
327 | // |
---|
328 | answer = (bpl-apl)/2.0*sum(fcnRet) // get the sum, normalize to parallel direction |
---|
329 | answer *= (bpp-app)/2.0 // and normalize to perpendicular direction |
---|
330 | |
---|
331 | s.zw[ii] = answer |
---|
332 | endfor |
---|
333 | |
---|
334 | Variable elap = (StopMSTimer(-2) - t1)/1e6 |
---|
335 | Print "elapsed time = ",elap |
---|
336 | |
---|
337 | return(0) |
---|
338 | end |
---|
339 | |
---|
340 | |
---|
341 | |
---|
342 | |
---|
343 | |
---|
344 | |
---|
345 | //phi is defined from +x axis, proceeding CCW around [0,2Pi] |
---|
346 | //rotate the resolution function by theta, = -phi |
---|
347 | // |
---|
348 | // this is only different by (-1) from FindPhi |
---|
349 | // I'd just call FindPhi, but it's awkward to include |
---|
350 | // |
---|
351 | Threadsafe Function FindTheta(vx,vy) |
---|
352 | variable vx,vy |
---|
353 | |
---|
354 | |
---|
355 | return(-1 * FindPhi(vx,vy)) |
---|
356 | end |
---|
357 | |
---|
358 | |
---|
359 | |
---|
360 | |
---|
361 | ////////////////////////////////////////////////////////// |
---|
362 | ////////////////////////////////////////////////////////// |
---|
363 | /// utility functions to test the calculation of the resolution function and |
---|
364 | // plotting of it for visualization |
---|
365 | // |
---|
366 | // -- these need to have the reduction package included for it to compile |
---|
367 | // |
---|
368 | // this works with the raw 2D data, calculating resolution in place |
---|
369 | // type is "RAW" |
---|
370 | // xx,yy are in detector coordinates |
---|
371 | // |
---|
372 | // call as: |
---|
373 | // PlotResolution_atPixel("RAW",pcsr(A),qcsr(A)) |
---|
374 | // |
---|
375 | // then: |
---|
376 | // Display;AppendImage res;AppendMatrixContour res;ModifyContour res labels=0,autoLevels={*,*,3} |
---|
377 | // |
---|
378 | |
---|
379 | Function PlotResolution_atPixel(type,xx,yy) |
---|
380 | String type |
---|
381 | Variable xx,yy |
---|
382 | |
---|
383 | ///from QxQyExport |
---|
384 | String destStr="",typeStr="" |
---|
385 | Variable step=1,refnum |
---|
386 | destStr = "root:Packages:NIST:"+type |
---|
387 | |
---|
388 | //must select the linear_data to export |
---|
389 | NVAR isLog = $(destStr+":gIsLogScale") |
---|
390 | if(isLog==1) |
---|
391 | typeStr = ":linear_data" |
---|
392 | else |
---|
393 | typeStr = ":data" |
---|
394 | endif |
---|
395 | |
---|
396 | NVAR pixelsX = root:myGlobals:gNPixelsX |
---|
397 | NVAR pixelsY = root:myGlobals:gNPixelsY |
---|
398 | |
---|
399 | Wave data=$(destStr+typeStr) |
---|
400 | WAVE intw=$(destStr + ":integersRead") |
---|
401 | WAVE rw=$(destStr + ":realsRead") |
---|
402 | WAVE/T textw=$(destStr + ":textRead") |
---|
403 | |
---|
404 | // Duplicate/O data,qx_val,qy_val,z_val,qval,qz_val,phi,r_dist |
---|
405 | Variable qx_val,qy_val,z_val,qval,qz_val,phi,r_dist |
---|
406 | |
---|
407 | Variable xctr,yctr,sdd,lambda,pixSize |
---|
408 | xctr = rw[16] |
---|
409 | yctr = rw[17] |
---|
410 | sdd = rw[18] |
---|
411 | lambda = rw[26] |
---|
412 | pixSize = rw[13]/10 //convert mm to cm (x and y are the same size pixels) |
---|
413 | |
---|
414 | qx_val = CalcQx(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) //+1 converts to detector coordinate system |
---|
415 | qy_val = CalcQy(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) |
---|
416 | |
---|
417 | Variable L2 = rw[18] |
---|
418 | Variable BS = rw[21] |
---|
419 | Variable S1 = rw[23] |
---|
420 | Variable S2 = rw[24] |
---|
421 | Variable L1 = rw[25] |
---|
422 | Variable lambdaWidth = rw[27] |
---|
423 | Variable usingLenses = rw[28] //new 2007 |
---|
424 | |
---|
425 | Variable vz_1 = 3.956e5 //velocity [cm/s] of 1 A neutron |
---|
426 | Variable g = 981.0 //gravity acceleration [cm/s^2] |
---|
427 | Variable m_h = 252.8 // m/h [=] s/cm^2 |
---|
428 | |
---|
429 | Variable acc,ssd,lambda0,yg_d,qstar |
---|
430 | |
---|
431 | G = 981. //! ACCELERATION OF GRAVITY, CM/SEC^2 |
---|
432 | acc = vz_1 // 3.956E5 //! CONVERT WAVELENGTH TO VELOCITY CM/SEC |
---|
433 | SDD = L2 *100 //1317 |
---|
434 | SSD = L1 *100 //1627 //cm |
---|
435 | lambda0 = lambda // 15 |
---|
436 | YG_d = -0.5*G*SDD*(SSD+SDD)*(LAMBDA0/acc)^2 |
---|
437 | Print "DISTANCE BEAM FALLS DUE TO GRAVITY (CM) = ",YG_d |
---|
438 | Print "Gravity q* = ",-2*pi/lambda0*2*yg_d/sdd |
---|
439 | qstar = -2*pi/lambda0*2*yg_d/sdd |
---|
440 | |
---|
441 | |
---|
442 | // the gravity center is not the resolution center |
---|
443 | // gravity center = beam center |
---|
444 | // resolution center = offset y = dy + (2)*yg_d |
---|
445 | |
---|
446 | |
---|
447 | qval = CalcQval(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) |
---|
448 | qz_val = CalcQz(xx+1,yy+1,rw[16],rw[17],rw[18],rw[26],rw[13]/10) |
---|
449 | phi = FindPhi( pixSize*((xx+1)-xctr) , pixSize*((yy+1)-yctr)+(2)*yg_d) //(dx,dy+yg_d) |
---|
450 | r_dist = sqrt( (pixSize*((xx+1)-xctr))^2 + (pixSize*((yy+1)-yctr)+(2)*yg_d)^2 ) //radial distance from ctr to pt |
---|
451 | |
---|
452 | Print pixSize*((yy+1)-yctr),pixSize*((yy+1)-yctr)+(2)*yg_d |
---|
453 | |
---|
454 | // Redimension/N=(pixelsX*pixelsY) qz_val,qval,phi,r_dist |
---|
455 | //everything in 1D now |
---|
456 | // Duplicate/O qval SigmaQX,SigmaQY,fsubS |
---|
457 | Variable SigmaQX,SigmaQY,fsubS |
---|
458 | |
---|
459 | //Two parameters DDET and APOFF are instrument dependent. Determine |
---|
460 | //these from the instrument name in the header. |
---|
461 | //From conversation with JB on 01.06.99 these are the current good values |
---|
462 | Variable DDet |
---|
463 | NVAR apOff = root:myGlobals:apOff //in cm |
---|
464 | DDet = rw[10]/10 // header value (X) is in mm, want cm here |
---|
465 | |
---|
466 | Variable ret1,ret2,ret3,del_r |
---|
467 | del_r = rw[10] |
---|
468 | |
---|
469 | get2DResolution(qval,phi,lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,del_r,usingLenses,r_dist,ret1,ret2,ret3) |
---|
470 | SigmaQX = ret1 |
---|
471 | SigmaQY = ret2 |
---|
472 | fsubs = ret3 |
---|
473 | ///// |
---|
474 | |
---|
475 | // Variable theta,phi,qx,qy,sx,sy,a,b,c,val,ii,jj,num=15,x0,y0 |
---|
476 | Variable theta,qx,qy,sx,sy,a,b,c,val,ii,jj,num=15,x0,y0,maxSig,nStdDev=3,normFactor |
---|
477 | Variable qx_ret,qy_ret |
---|
478 | |
---|
479 | // theta = FindPhi(qx_val,qy_val) |
---|
480 | // need to rotate properly - theta is defined a starting from +y axis, moving CW |
---|
481 | // we define phi starting from +x and moving CCW |
---|
482 | theta = -phi //seems to give the right behavior... |
---|
483 | |
---|
484 | |
---|
485 | Print qx_val,qy_val,qval |
---|
486 | Print "phi, theta",phi,theta |
---|
487 | |
---|
488 | // FindQxQy(qval,phi,qx_ret,qy_ret) |
---|
489 | |
---|
490 | sx = SigmaQx |
---|
491 | sy = sigmaQy |
---|
492 | x0 = qx_val |
---|
493 | y0 = qy_val |
---|
494 | |
---|
495 | a = cos(theta)^2/(2*sx*sx) + sin(theta)^2/(2*sy*sy) |
---|
496 | b = -1*sin(2*theta)/(4*sx*sx) + sin(2*theta)/(4*sy*sy) |
---|
497 | c = sin(theta)^2/(2*sx*sx) + cos(theta)^2/(2*sy*sy) |
---|
498 | |
---|
499 | normFactor = pi/sqrt(a*c-b*b) |
---|
500 | |
---|
501 | Make/O/D/N=(num,num) res |
---|
502 | // so the resolution function 'looks right' on a 2D plot - otherwise it always looks like a circle |
---|
503 | maxSig = max(sx,sy) |
---|
504 | Setscale/I x -nStdDev*maxSig+x0,nStdDev*maxSig+x0,res |
---|
505 | Setscale/I y -nStdDev*maxSig+y0,nStdDev*maxSig+y0,res |
---|
506 | ///// Setscale/I x -nStdDev*sx+x0,nStdDev*sx+x0,res |
---|
507 | ///// Setscale/I y -nStdDev*sy+y0,nStdDev*sy+y0,res |
---|
508 | |
---|
509 | Variable xPt,yPt,delx,dely,offx,offy |
---|
510 | delx = DimDelta(res,0) |
---|
511 | dely = DimDelta(res,1) |
---|
512 | offx = DimOffset(res,0) |
---|
513 | offy = DimOffset(res,1) |
---|
514 | |
---|
515 | Print "sx,sy = ",sx,sy |
---|
516 | for(ii=0;ii<num;ii+=1) |
---|
517 | xPt = offx + ii*delx |
---|
518 | for(jj=0;jj<num;jj+=1) |
---|
519 | yPt = offy + jj*dely |
---|
520 | res[ii][jj] = exp(-1*(a*(xPt-x0)^2 + 2*b*(xPt-x0)*(yPt-y0) + c*(yPt-y0)^2)) |
---|
521 | endfor |
---|
522 | endfor |
---|
523 | res /= normFactor |
---|
524 | |
---|
525 | //Print sum(res,-inf,inf)*delx*dely |
---|
526 | if(WaveExists($"coef")==0) |
---|
527 | Make/O/D/N=6 coef |
---|
528 | endif |
---|
529 | Wave coef=coef |
---|
530 | coef[0] = 1 |
---|
531 | coef[1] = qx_val |
---|
532 | coef[2] = qy_val |
---|
533 | coef[3] = sx |
---|
534 | coef[4] = sy |
---|
535 | coef[5] = theta |
---|
536 | |
---|
537 | // Variable t1=StopMSTimer(-2) |
---|
538 | |
---|
539 | // |
---|
540 | do2dIntegrationGauss(-nStdDev*maxSig+x0,nStdDev*maxSig+x0,-nStdDev*maxSig+y0,nStdDev*maxSig+y0) |
---|
541 | // |
---|
542 | |
---|
543 | // Variable elap = (StopMSTimer(-2) - t1)/1e6 |
---|
544 | // Print "elapsed time = ",elap |
---|
545 | // Print "time for 16384 = (minutes)",16384*elap/60 |
---|
546 | return(0) |
---|
547 | End |
---|
548 | |
---|
549 | |
---|
550 | // this is called each time to integrate the gaussian |
---|
551 | Function do2dIntegrationGauss(xMin,xMax,yMin,yMax) |
---|
552 | Variable xMin,xMax,yMin,yMax |
---|
553 | |
---|
554 | Variable/G globalXmin=xMin |
---|
555 | Variable/G globalXmax=xMax |
---|
556 | Variable/G globalY |
---|
557 | |
---|
558 | Variable result=Integrate1d(Gauss2DFuncOuter,yMin,yMax,2,5) |
---|
559 | KillVariables/z globalXmax,globalXmin,globalY |
---|
560 | print "integration of 2D = ",result |
---|
561 | End |
---|
562 | |
---|
563 | Function Gauss2DFuncOuter(inY) |
---|
564 | Variable inY |
---|
565 | |
---|
566 | NVAR globalXmin,globalXmax,globalY |
---|
567 | globalY=inY |
---|
568 | |
---|
569 | return integrate1D(Gauss2DFuncInner,globalXmin,globalXmax,2,5) |
---|
570 | End |
---|
571 | |
---|
572 | Function Gauss2DFuncInner(inX) |
---|
573 | Variable inX |
---|
574 | |
---|
575 | NVAR globalY |
---|
576 | Wave coef=coef |
---|
577 | |
---|
578 | return Gauss2D_theta(coef,inX,GlobalY) |
---|
579 | End |
---|
580 | |
---|
581 | Function Gauss2D_theta(w,x,y) |
---|
582 | Wave w |
---|
583 | Variable x,y |
---|
584 | |
---|
585 | Variable val,a,b,c |
---|
586 | Variable scale,x0,y0,sx,sy,theta,normFactor |
---|
587 | |
---|
588 | scale = w[0] |
---|
589 | x0 = w[1] |
---|
590 | y0 = w[2] |
---|
591 | sx = w[3] |
---|
592 | sy = w[4] |
---|
593 | theta = w[5] |
---|
594 | |
---|
595 | a = cos(theta)^2/(2*sx*sx) + sin(theta)^2/(2*sy*sy) |
---|
596 | b = -1*sin(2*theta)/(4*sx*sx) + sin(2*theta)/(4*sy*sy) |
---|
597 | c = sin(theta)^2/(2*sx*sx) + cos(theta)^2/(2*sy*sy) |
---|
598 | |
---|
599 | val = exp(-1*(a*(x-x0)^2 + 2*b*(x-x0)*(y-y0) + c*(y-y0)^2)) |
---|
600 | |
---|
601 | normFactor = pi/sqrt(a*c-b*b) |
---|
602 | |
---|
603 | return(scale*val/normFactor) |
---|
604 | end |
---|
605 | |
---|
606 | //////////////////////////////////////////////////////////// |
---|
607 | //////////////////////////////////////////////////////////// |
---|