1 | #pragma rtGlobals=1 // Use modern global access method. |
---|
2 | #pragma version=4.00 |
---|
3 | #pragma IgorVersion=6.1 |
---|
4 | |
---|
5 | // GaussUtils.ipf |
---|
6 | // 22dec97 SRK |
---|
7 | |
---|
8 | //// NEW |
---|
9 | // |
---|
10 | // added two utility functions to do numerical |
---|
11 | // integration of an arbitrary function |
---|
12 | // Gaussian quadrature of 20 or 76 points is used |
---|
13 | // |
---|
14 | // fcn must be defined as in the prototype function |
---|
15 | // which is similar to the normal fitting function definition |
---|
16 | // |
---|
17 | // 09 SEP 03 SRK |
---|
18 | // |
---|
19 | // |
---|
20 | // 04DEC03 - added routines for 5 and 10 Gauss points |
---|
21 | // |
---|
22 | // 13 DEC 04 BSG/SRK - Modified Smear_Model_(5)(20) to allow |
---|
23 | // smearing of USANS data. dQv is passed in as negative value |
---|
24 | // in all resolution columns. the dQv value is read from the SigQ columnn |
---|
25 | // (the 4th column) - and all values are read, so fill the whole column |
---|
26 | // |
---|
27 | // Adaptive trapezoidal integration is used, 76 Gauss pts |
---|
28 | // did not give sufficient accuracy. |
---|
29 | // |
---|
30 | //////////////////////////////////////////////// |
---|
31 | // 23 JUL 2007 SRK |
---|
32 | // major overhaul to use AAO functions in the smearing calculations |
---|
33 | // |
---|
34 | // - re-definition of function prototype in Smear_Model_N |
---|
35 | // - re-definition in trapezoidal routines too |
---|
36 | // - calls from model functions are somewhat different, so this will generate a lot of errors |
---|
37 | // until all can be changed |
---|
38 | // - now is looking for a resolution matrix that contains the 3 resolution waves plus Q |
---|
39 | // = mat[num_Qvals][0] = sigQ |
---|
40 | // = mat[num_Qvals][1] = Qbar |
---|
41 | // = mat[num_Qvals][2] = fShad |
---|
42 | // = mat[num_Qvals][3] = qvals |
---|
43 | // |
---|
44 | // -- does not yet use the matrix calculation for USANS |
---|
45 | // -- SO THERE IS NO SWITCH YET IN Smear_Model_N() |
---|
46 | // |
---|
47 | // |
---|
48 | |
---|
49 | //maybe not the optimal set of parameters for the STRUCT |
---|
50 | // |
---|
51 | // resW is /N=(np,4) if SANS data |
---|
52 | // resW is /N=(np,np) if USANS data |
---|
53 | // |
---|
54 | // info may be useful later as a "KEY=value;" string to carry additional information |
---|
55 | Structure ResSmearAAOStruct |
---|
56 | Wave coefW |
---|
57 | Wave yW |
---|
58 | Wave xW |
---|
59 | Wave resW |
---|
60 | String info |
---|
61 | EndStructure |
---|
62 | |
---|
63 | |
---|
64 | //// tentative pass at 2D resolution smearing |
---|
65 | //// |
---|
66 | //Structure ResSmear_2D_AAOStruct |
---|
67 | // Wave coefW |
---|
68 | // Wave zw //answer |
---|
69 | // Wave qy // q-value |
---|
70 | // Wave qx |
---|
71 | // Wave qz |
---|
72 | // Wave sQpl //resolution parallel to Q |
---|
73 | // Wave sQpp //resolution perpendicular to Q |
---|
74 | // Wave fs |
---|
75 | // String info |
---|
76 | //EndStructure |
---|
77 | |
---|
78 | // reformat the structure ?? WM Fit compatible |
---|
79 | // -- 2D resolution smearing |
---|
80 | // |
---|
81 | Structure ResSmear_2D_AAOStruct |
---|
82 | Wave coefW |
---|
83 | Wave zw //answer |
---|
84 | Wave xw[2] // qx-value is [0], qy is xw[1] |
---|
85 | Wave qz |
---|
86 | Wave sQpl //resolution parallel to Q |
---|
87 | Wave sQpp //resolution perpendicular to Q |
---|
88 | Wave fs |
---|
89 | String info |
---|
90 | EndStructure |
---|
91 | |
---|
92 | |
---|
93 | |
---|
94 | Function Make5GaussPoints(w5,z5) |
---|
95 | Wave w5,z5 |
---|
96 | |
---|
97 | // printf "in make Gauss Pts\r" |
---|
98 | |
---|
99 | z5[0] = -.906179845938664 |
---|
100 | z5[1] = -.538469310105683 |
---|
101 | z5[2] = -.0000000000000 |
---|
102 | z5[3] = .538469310105683 |
---|
103 | z5[4] = .906179845938664 |
---|
104 | |
---|
105 | w5[0] = .236926885056189 |
---|
106 | w5[1] = .478628670499366 |
---|
107 | w5[2] = .56888888888889 |
---|
108 | w5[3] = .478628670499366 |
---|
109 | w5[4] = .236926885056189 |
---|
110 | |
---|
111 | // printf "w[0],z[0] = %g %g\r", w5[0],z5[0] |
---|
112 | End |
---|
113 | |
---|
114 | Function Make10GaussPoints(w10,z10) |
---|
115 | Wave w10,z10 |
---|
116 | |
---|
117 | // printf "in make Gauss Pts\r" |
---|
118 | z10[0] = -.973906528517172 |
---|
119 | z10[1] = -.865063366688985 |
---|
120 | z10[2] = -.679409568299024 |
---|
121 | z10[3] = -.433395394129247 |
---|
122 | z10[4] = -.148874338981631 |
---|
123 | z10[5] = .148874338981631 |
---|
124 | z10[6] = .433395394129247 |
---|
125 | z10[7] = .679409568299024 |
---|
126 | z10[8] = .865063366688985 |
---|
127 | z10[9] = .973906528517172 |
---|
128 | |
---|
129 | w10[0] = .066671344308688 |
---|
130 | w10[1] = 0.149451349150581 |
---|
131 | w10[2] = 0.219086362515982 |
---|
132 | w10[3] = .269266719309996 |
---|
133 | w10[4] = 0.295524224714753 |
---|
134 | w10[5] = 0.295524224714753 |
---|
135 | w10[6] = .269266719309996 |
---|
136 | w10[7] = 0.219086362515982 |
---|
137 | w10[8] = 0.149451349150581 |
---|
138 | w10[9] = .066671344308688 |
---|
139 | |
---|
140 | // printf "w[0],z[0] = %g %g\r", w10[0],z10[0] |
---|
141 | End |
---|
142 | |
---|
143 | Function Make20GaussPoints(w20,z20) |
---|
144 | Wave w20,z20 |
---|
145 | |
---|
146 | // printf "in make Gauss Pts\r" |
---|
147 | |
---|
148 | z20[0] = -.993128599185095 |
---|
149 | z20[1] = -.963971927277914 |
---|
150 | z20[2] = -.912234428251326 |
---|
151 | z20[3] = -.839116971822219 |
---|
152 | z20[4] = -.746331906460151 |
---|
153 | z20[5] = -.636053680726515 |
---|
154 | z20[6] = -.510867001950827 |
---|
155 | z20[7] = -.37370608871542 |
---|
156 | z20[8] = -.227785851141645 |
---|
157 | z20[9] = -.076526521133497 |
---|
158 | z20[10] = .0765265211334973 |
---|
159 | z20[11] = .227785851141645 |
---|
160 | z20[12] = .37370608871542 |
---|
161 | z20[13] = .510867001950827 |
---|
162 | z20[14] = .636053680726515 |
---|
163 | z20[15] = .746331906460151 |
---|
164 | z20[16] = .839116971822219 |
---|
165 | z20[17] = .912234428251326 |
---|
166 | z20[18] = .963971927277914 |
---|
167 | z20[19] = .993128599185095 |
---|
168 | |
---|
169 | w20[0] = .0176140071391521 |
---|
170 | w20[1] = .0406014298003869 |
---|
171 | w20[2] = .0626720483341091 |
---|
172 | w20[3] = .0832767415767047 |
---|
173 | w20[4] = .10193011981724 |
---|
174 | w20[5] = .118194531961518 |
---|
175 | w20[6] = .131688638449177 |
---|
176 | w20[7] = .142096109318382 |
---|
177 | w20[8] = .149172986472604 |
---|
178 | w20[9] = .152753387130726 |
---|
179 | w20[10] = .152753387130726 |
---|
180 | w20[11] = .149172986472604 |
---|
181 | w20[12] = .142096109318382 |
---|
182 | w20[13] = .131688638449177 |
---|
183 | w20[14] = .118194531961518 |
---|
184 | w20[15] = .10193011981724 |
---|
185 | w20[16] = .0832767415767047 |
---|
186 | w20[17] = .0626720483341091 |
---|
187 | w20[18] = .0406014298003869 |
---|
188 | w20[19] = .0176140071391521 |
---|
189 | // printf "w[0],z[0] = %g %g\r", w20[0],z20[0] |
---|
190 | End |
---|
191 | |
---|
192 | Function Make76GaussPoints(w76,z76) |
---|
193 | Wave w76,z76 |
---|
194 | |
---|
195 | // printf "in make Gauss Pts\r" |
---|
196 | |
---|
197 | z76[0] = .999505948362153*(-1.0) |
---|
198 | z76[75] = -.999505948362153*(-1.0) |
---|
199 | z76[1] = .997397786355355*(-1.0) |
---|
200 | z76[74] = -.997397786355355*(-1.0) |
---|
201 | z76[2] = .993608772723527*(-1.0) |
---|
202 | z76[73] = -.993608772723527*(-1.0) |
---|
203 | z76[3] = .988144453359837*(-1.0) |
---|
204 | z76[72] = -.988144453359837*(-1.0) |
---|
205 | z76[4] = .981013938975656*(-1.0) |
---|
206 | z76[71] = -.981013938975656*(-1.0) |
---|
207 | z76[5] = .972229228520377*(-1.0) |
---|
208 | z76[70] = -.972229228520377*(-1.0) |
---|
209 | z76[6] = .961805126758768*(-1.0) |
---|
210 | z76[69] = -.961805126758768*(-1.0) |
---|
211 | z76[7] = .949759207710896*(-1.0) |
---|
212 | z76[68] = -.949759207710896*(-1.0) |
---|
213 | z76[8] = .936111781934811*(-1.0) |
---|
214 | z76[67] = -.936111781934811*(-1.0) |
---|
215 | z76[9] = .92088586125215*(-1.0) |
---|
216 | z76[66] = -.92088586125215*(-1.0) |
---|
217 | z76[10] = .904107119545567*(-1.0) |
---|
218 | z76[65] = -.904107119545567*(-1.0) |
---|
219 | z76[11] = .885803849292083*(-1.0) |
---|
220 | z76[64] = -.885803849292083*(-1.0) |
---|
221 | z76[12] = .866006913771982*(-1.0) |
---|
222 | z76[63] = -.866006913771982*(-1.0) |
---|
223 | z76[13] = .844749694983342*(-1.0) |
---|
224 | z76[62] = -.844749694983342*(-1.0) |
---|
225 | z76[14] = .822068037328975*(-1.0) |
---|
226 | z76[61] = -.822068037328975*(-1.0) |
---|
227 | z76[15] = .7980001871612*(-1.0) |
---|
228 | z76[60] = -.7980001871612*(-1.0) |
---|
229 | z76[16] = .77258672828181*(-1.0) |
---|
230 | z76[59] = -.77258672828181*(-1.0) |
---|
231 | z76[17] = .74587051350361*(-1.0) |
---|
232 | z76[58] = -.74587051350361*(-1.0) |
---|
233 | z76[18] = .717896592387704*(-1.0) |
---|
234 | z76[57] = -.717896592387704*(-1.0) |
---|
235 | z76[19] = .688712135277641*(-1.0) |
---|
236 | z76[56] = -.688712135277641*(-1.0) |
---|
237 | z76[20] = .658366353758143*(-1.0) |
---|
238 | z76[55] = -.658366353758143*(-1.0) |
---|
239 | z76[21] = .626910417672267*(-1.0) |
---|
240 | z76[54] = -.626910417672267*(-1.0) |
---|
241 | z76[22] = .594397368836793*(-1.0) |
---|
242 | z76[53] = -.594397368836793*(-1.0) |
---|
243 | z76[23] = .560882031601237*(-1.0) |
---|
244 | z76[52] = -.560882031601237*(-1.0) |
---|
245 | z76[24] = .526420920401243*(-1.0) |
---|
246 | z76[51] = -.526420920401243*(-1.0) |
---|
247 | z76[25] = .491072144462194*(-1.0) |
---|
248 | z76[50] = -.491072144462194*(-1.0) |
---|
249 | z76[26] = .454895309813726*(-1.0) |
---|
250 | z76[49] = -.454895309813726*(-1.0) |
---|
251 | z76[27] = .417951418780327*(-1.0) |
---|
252 | z76[48] = -.417951418780327*(-1.0) |
---|
253 | z76[28] = .380302767117504*(-1.0) |
---|
254 | z76[47] = -.380302767117504*(-1.0) |
---|
255 | z76[29] = .342012838966962*(-1.0) |
---|
256 | z76[46] = -.342012838966962*(-1.0) |
---|
257 | z76[30] = .303146199807908*(-1.0) |
---|
258 | z76[45] = -.303146199807908*(-1.0) |
---|
259 | z76[31] = .263768387584994*(-1.0) |
---|
260 | z76[44] = -.263768387584994*(-1.0) |
---|
261 | z76[32] = .223945802196474*(-1.0) |
---|
262 | z76[43] = -.223945802196474*(-1.0) |
---|
263 | z76[33] = .183745593528914*(-1.0) |
---|
264 | z76[42] = -.183745593528914*(-1.0) |
---|
265 | z76[34] = .143235548227268*(-1.0) |
---|
266 | z76[41] = -.143235548227268*(-1.0) |
---|
267 | z76[35] = .102483975391227*(-1.0) |
---|
268 | z76[40] = -.102483975391227*(-1.0) |
---|
269 | z76[36] = .0615595913906112*(-1.0) |
---|
270 | z76[39] = -.0615595913906112*(-1.0) |
---|
271 | z76[37] = .0205314039939986*(-1.0) |
---|
272 | z76[38] = -.0205314039939986*(-1.0) |
---|
273 | |
---|
274 | w76[0] = .00126779163408536 |
---|
275 | w76[75] = .00126779163408536 |
---|
276 | w76[1] = .00294910295364247 |
---|
277 | w76[74] = .00294910295364247 |
---|
278 | w76[2] = .00462793522803742 |
---|
279 | w76[73] = .00462793522803742 |
---|
280 | w76[3] = .00629918049732845 |
---|
281 | w76[72] = .00629918049732845 |
---|
282 | w76[4] = .00795984747723973 |
---|
283 | w76[71] = .00795984747723973 |
---|
284 | w76[5] = .00960710541471375 |
---|
285 | w76[70] = .00960710541471375 |
---|
286 | w76[6] = .0112381685696677 |
---|
287 | w76[69] = .0112381685696677 |
---|
288 | w76[7] = .0128502838475101 |
---|
289 | w76[68] = .0128502838475101 |
---|
290 | w76[8] = .0144407317482767 |
---|
291 | w76[67] = .0144407317482767 |
---|
292 | w76[9] = .0160068299122486 |
---|
293 | w76[66] = .0160068299122486 |
---|
294 | w76[10] = .0175459372914742 |
---|
295 | w76[65] = .0175459372914742 |
---|
296 | w76[11] = .0190554584671906 |
---|
297 | w76[64] = .0190554584671906 |
---|
298 | w76[12] = .020532847967908 |
---|
299 | w76[63] = .020532847967908 |
---|
300 | w76[13] = .0219756145344162 |
---|
301 | w76[62] = .0219756145344162 |
---|
302 | w76[14] = .0233813253070112 |
---|
303 | w76[61] = .0233813253070112 |
---|
304 | w76[15] = .0247476099206597 |
---|
305 | w76[60] = .0247476099206597 |
---|
306 | w76[16] = .026072164497986 |
---|
307 | w76[59] = .026072164497986 |
---|
308 | w76[17] = .0273527555318275 |
---|
309 | w76[58] = .0273527555318275 |
---|
310 | w76[18] = .028587223650054 |
---|
311 | w76[57] = .028587223650054 |
---|
312 | w76[19] = .029773487255905 |
---|
313 | w76[56] = .029773487255905 |
---|
314 | w76[20] = .0309095460374916 |
---|
315 | w76[55] = .0309095460374916 |
---|
316 | w76[21] = .0319934843404216 |
---|
317 | w76[54] = .0319934843404216 |
---|
318 | w76[22] = .0330234743977917 |
---|
319 | w76[53] = .0330234743977917 |
---|
320 | w76[23] = .0339977794120564 |
---|
321 | w76[52] = .0339977794120564 |
---|
322 | w76[24] = .0349147564835508 |
---|
323 | w76[51] = .0349147564835508 |
---|
324 | w76[25] = .0357728593807139 |
---|
325 | w76[50] = .0357728593807139 |
---|
326 | w76[26] = .0365706411473296 |
---|
327 | w76[49] = .0365706411473296 |
---|
328 | w76[27] = .0373067565423816 |
---|
329 | w76[48] = .0373067565423816 |
---|
330 | w76[28] = .0379799643084053 |
---|
331 | w76[47] = .0379799643084053 |
---|
332 | w76[29] = .0385891292645067 |
---|
333 | w76[46] = .0385891292645067 |
---|
334 | w76[30] = .0391332242205184 |
---|
335 | w76[45] = .0391332242205184 |
---|
336 | w76[31] = .0396113317090621 |
---|
337 | w76[44] = .0396113317090621 |
---|
338 | w76[32] = .0400226455325968 |
---|
339 | w76[43] = .0400226455325968 |
---|
340 | w76[33] = .040366472122844 |
---|
341 | w76[42] = .040366472122844 |
---|
342 | w76[34] = .0406422317102947 |
---|
343 | w76[41] = .0406422317102947 |
---|
344 | w76[35] = .0408494593018285 |
---|
345 | w76[40] = .0408494593018285 |
---|
346 | w76[36] = .040987805464794 |
---|
347 | w76[39] = .040987805464794 |
---|
348 | w76[37] = .0410570369162294 |
---|
349 | w76[38] = .0410570369162294 |
---|
350 | |
---|
351 | End //Make76GaussPoints() |
---|
352 | |
---|
353 | // !!!!! reduces the length of wt and zi by one !!!!! |
---|
354 | // |
---|
355 | Function Make_N_GaussPoints(wt,zi) |
---|
356 | Wave wt,zi |
---|
357 | |
---|
358 | Variable num |
---|
359 | num = numpnts(wt) - 1 |
---|
360 | |
---|
361 | gauleg(-1,1,zi,wt,num) |
---|
362 | |
---|
363 | DeletePoints 0,1,wt,zi |
---|
364 | |
---|
365 | return(0) |
---|
366 | End |
---|
367 | |
---|
368 | /// gauleg subroutine from NR to calculate weights and abscissae for |
---|
369 | // Gauss-Legendre quadrature |
---|
370 | // |
---|
371 | // |
---|
372 | // arrays are indexed from 1 |
---|
373 | // |
---|
374 | Function gauleg( x1, x2, x, w, n) |
---|
375 | Variable x1, x2 |
---|
376 | Wave x, w |
---|
377 | Variable n |
---|
378 | |
---|
379 | variable m,j,i |
---|
380 | variable z1,z,xm,xl,pp,p3,p2,p1 |
---|
381 | Variable eps = 3e-11 |
---|
382 | |
---|
383 | m=(n+1)/2 |
---|
384 | xm=0.5*(x2+x1) |
---|
385 | xl=0.5*(x2-x1) |
---|
386 | for (i=1;i<=m;i+=1) |
---|
387 | z=cos(pi*(i-0.25)/(n+0.5)) |
---|
388 | do |
---|
389 | p1=1.0 |
---|
390 | p2=0.0 |
---|
391 | for (j=1;j<=n;j+=1) |
---|
392 | p3=p2 |
---|
393 | p2=p1 |
---|
394 | p1=((2.0*j-1.0)*z*p2-(j-1.0)*p3)/j |
---|
395 | endfor |
---|
396 | pp=n*(z*p1-p2)/(z*z-1.0) |
---|
397 | z1=z |
---|
398 | z=z1-p1/pp |
---|
399 | while (abs(z-z1) > EPS) |
---|
400 | x[i]=xm-xl*z |
---|
401 | x[n+1-i]=xm+xl*z |
---|
402 | w[i]=2.0*xl/((1.0-z*z)*pp*pp) |
---|
403 | w[n+1-i]=w[i] |
---|
404 | Endfor |
---|
405 | End |
---|
406 | |
---|
407 | |
---|
408 | /// uses a user-supplied number of Gauss points, and generates them on-the-fly as needed |
---|
409 | // using a Numerical Recipes routine |
---|
410 | // |
---|
411 | // - note that this has an extra input parameter, nord |
---|
412 | // |
---|
413 | //////////// |
---|
414 | Function IntegrateFn_N(fcn,loLim,upLim,w,x,nord) |
---|
415 | FUNCREF GenericQuadrature_proto fcn |
---|
416 | Variable loLim,upLim //limits of integration |
---|
417 | Wave w //coefficients of function fcn(w,x) |
---|
418 | Variable x //x-value (q) for the calculation |
---|
419 | Variable nord //number of quadrature points to used |
---|
420 | |
---|
421 | // local variables |
---|
422 | Variable ii,va,vb,summ,yyy,zi |
---|
423 | Variable answer,dum |
---|
424 | String weightStr,zStr |
---|
425 | |
---|
426 | weightStr = "gauss"+num2iStr(nord)+"wt" |
---|
427 | zStr = "gauss"+num2istr(nord)+"z" |
---|
428 | |
---|
429 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
430 | Make/D/N=(nord+1) $weightStr,$zStr |
---|
431 | Wave wt = $weightStr |
---|
432 | Wave xx = $zStr // wave references to pass |
---|
433 | Make_N_GaussPoints(wt,xx) //generates the gauss points and removes the extra point |
---|
434 | else |
---|
435 | if(exists(weightStr) > 1) |
---|
436 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
437 | endif |
---|
438 | Wave wt = $weightStr |
---|
439 | Wave xx = $zStr // create the wave references |
---|
440 | endif |
---|
441 | |
---|
442 | //limits of integration are input to function |
---|
443 | va = loLim |
---|
444 | vb = upLim |
---|
445 | // Using 5 Gauss points |
---|
446 | // remember to index from 0,size-1 |
---|
447 | |
---|
448 | summ = 0.0 // initialize integral |
---|
449 | ii=0 // loop counter |
---|
450 | do |
---|
451 | // calculate Gauss points on integration interval (q-value for evaluation) |
---|
452 | zi = ( xx[ii]*(vb-va) + vb + va )/2.0 |
---|
453 | //calculate partial sum for the passed-in model function |
---|
454 | yyy = wt[ii] * fcn(w,x,zi) |
---|
455 | summ += yyy //add to the running total of the quadrature |
---|
456 | ii+=1 |
---|
457 | while (ii<nord) // end of loop over quadrature points |
---|
458 | |
---|
459 | // calculate value of integral to return |
---|
460 | answer = (vb-va)/2.0*summ |
---|
461 | |
---|
462 | Return (answer) |
---|
463 | End |
---|
464 | |
---|
465 | |
---|
466 | |
---|
467 | //////////// |
---|
468 | Function IntegrateFn5(fcn,loLim,upLim,w,x) |
---|
469 | FUNCREF GenericQuadrature_proto fcn |
---|
470 | Variable loLim,upLim //limits of integration |
---|
471 | Wave w //coefficients of function fcn(w,x) |
---|
472 | Variable x //x-value (q) for the calculation |
---|
473 | |
---|
474 | // local variables |
---|
475 | Variable nord,ii,va,vb,summ,yyy,zi |
---|
476 | Variable answer,dum |
---|
477 | String weightStr,zStr |
---|
478 | |
---|
479 | weightStr = "gauss5wt" |
---|
480 | zStr = "gauss5z" |
---|
481 | nord = 5 |
---|
482 | |
---|
483 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
484 | Make/D/N=5 $weightStr,$zStr |
---|
485 | Wave w5 = $weightStr |
---|
486 | Wave z5 = $zStr // wave references to pass |
---|
487 | Make5GaussPoints(w5,z5) |
---|
488 | else |
---|
489 | if(exists(weightStr) > 1) |
---|
490 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
491 | endif |
---|
492 | Wave w5 = $weightStr |
---|
493 | Wave z5 = $zStr // create the wave references |
---|
494 | endif |
---|
495 | |
---|
496 | //limits of integration are input to function |
---|
497 | va = loLim |
---|
498 | vb = upLim |
---|
499 | // Using 5 Gauss points |
---|
500 | // remember to index from 0,size-1 |
---|
501 | |
---|
502 | summ = 0.0 // initialize integral |
---|
503 | ii=0 // loop counter |
---|
504 | do |
---|
505 | // calculate Gauss points on integration interval (q-value for evaluation) |
---|
506 | zi = ( z5[ii]*(vb-va) + vb + va )/2.0 |
---|
507 | //calculate partial sum for the passed-in model function |
---|
508 | yyy = w5[ii] * fcn(w,x,zi) |
---|
509 | summ += yyy //add to the running total of the quadrature |
---|
510 | ii+=1 |
---|
511 | while (ii<nord) // end of loop over quadrature points |
---|
512 | |
---|
513 | // calculate value of integral to return |
---|
514 | answer = (vb-va)/2.0*summ |
---|
515 | |
---|
516 | Return (answer) |
---|
517 | End |
---|
518 | /////////////////////////////////////////////////////////////// |
---|
519 | |
---|
520 | //////////// |
---|
521 | Function IntegrateFn10(fcn,loLim,upLim,w,x) |
---|
522 | FUNCREF GenericQuadrature_proto fcn |
---|
523 | Variable loLim,upLim //limits of integration |
---|
524 | Wave w //coefficients of function fcn(w,x) |
---|
525 | Variable x //x-value (q) for the calculation |
---|
526 | |
---|
527 | // local variables |
---|
528 | Variable nord,ii,va,vb,summ,yyy,zi |
---|
529 | Variable answer,dum |
---|
530 | String weightStr,zStr |
---|
531 | |
---|
532 | weightStr = "gauss10wt" |
---|
533 | zStr = "gauss10z" |
---|
534 | nord = 10 |
---|
535 | |
---|
536 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
537 | Make/D/N=10 $weightStr,$zStr |
---|
538 | Wave w10 = $weightStr |
---|
539 | Wave z10 = $zStr // wave references to pass |
---|
540 | Make10GaussPoints(w10,z10) |
---|
541 | else |
---|
542 | if(exists(weightStr) > 1) |
---|
543 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
544 | endif |
---|
545 | Wave w10 = $weightStr |
---|
546 | Wave z10 = $zStr // create the wave references |
---|
547 | endif |
---|
548 | |
---|
549 | //limits of integration are input to function |
---|
550 | va = loLim |
---|
551 | vb = upLim |
---|
552 | // Using 10 Gauss points |
---|
553 | // remember to index from 0,size-1 |
---|
554 | |
---|
555 | summ = 0.0 // initialize integral |
---|
556 | ii=0 // loop counter |
---|
557 | do |
---|
558 | // calculate Gauss points on integration interval (q-value for evaluation) |
---|
559 | zi = ( z10[ii]*(vb-va) + vb + va )/2.0 |
---|
560 | //calculate partial sum for the passed-in model function |
---|
561 | yyy = w10[ii] * fcn(w,x,zi) |
---|
562 | summ += yyy //add to the running total of the quadrature |
---|
563 | ii+=1 |
---|
564 | while (ii<nord) // end of loop over quadrature points |
---|
565 | |
---|
566 | // calculate value of integral to return |
---|
567 | answer = (vb-va)/2.0*summ |
---|
568 | |
---|
569 | Return (answer) |
---|
570 | End |
---|
571 | /////////////////////////////////////////////////////////////// |
---|
572 | |
---|
573 | //////////// |
---|
574 | Function IntegrateFn20(fcn,loLim,upLim,w,x) |
---|
575 | FUNCREF GenericQuadrature_proto fcn |
---|
576 | Variable loLim,upLim //limits of integration |
---|
577 | Wave w //coefficients of function fcn(w,x) |
---|
578 | Variable x //x-value (q) for the calculation |
---|
579 | |
---|
580 | // local variables |
---|
581 | Variable nord,ii,va,vb,summ,yyy,zi |
---|
582 | Variable answer,dum |
---|
583 | String weightStr,zStr |
---|
584 | |
---|
585 | weightStr = "gauss20wt" |
---|
586 | zStr = "gauss20z" |
---|
587 | nord = 20 |
---|
588 | |
---|
589 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
590 | Make/D/N=20 $weightStr,$zStr |
---|
591 | Wave w20 = $weightStr |
---|
592 | Wave z20 = $zStr // wave references to pass |
---|
593 | Make20GaussPoints(w20,z20) |
---|
594 | else |
---|
595 | if(exists(weightStr) > 1) |
---|
596 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
597 | endif |
---|
598 | Wave w20 = $weightStr |
---|
599 | Wave z20 = $zStr // create the wave references |
---|
600 | endif |
---|
601 | |
---|
602 | //limits of integration are input to function |
---|
603 | va = loLim |
---|
604 | vb = upLim |
---|
605 | // Using 20 Gauss points |
---|
606 | // remember to index from 0,size-1 |
---|
607 | |
---|
608 | summ = 0.0 // initialize integral |
---|
609 | ii=0 // loop counter |
---|
610 | do |
---|
611 | // calculate Gauss points on integration interval (q-value for evaluation) |
---|
612 | zi = ( z20[ii]*(vb-va) + vb + va )/2.0 |
---|
613 | //calculate partial sum for the passed-in model function |
---|
614 | yyy = w20[ii] * fcn(w,x,zi) |
---|
615 | summ += yyy //add to the running total of the quadrature |
---|
616 | ii+=1 |
---|
617 | while (ii<nord) // end of loop over quadrature points |
---|
618 | |
---|
619 | // calculate value of integral to return |
---|
620 | answer = (vb-va)/2.0*summ |
---|
621 | |
---|
622 | Return (answer) |
---|
623 | End |
---|
624 | /////////////////////////////////////////////////////////////// |
---|
625 | |
---|
626 | Function IntegrateFn76(fcn,loLim,upLim,w,x) |
---|
627 | FUNCREF GenericQuadrature_proto fcn |
---|
628 | Variable loLim,upLim //limits of integration |
---|
629 | Wave w //coefficients of function fcn(w,x) |
---|
630 | Variable x //x-value (q) for the calculation |
---|
631 | |
---|
632 | //**** The coefficient wave is passed into this function and straight through to the unsmeared model function |
---|
633 | |
---|
634 | // local variables |
---|
635 | Variable nord,ii,va,vb,summ,yyy,zi |
---|
636 | Variable answer,dum |
---|
637 | String weightStr,zStr |
---|
638 | |
---|
639 | weightStr = "gauss76wt" |
---|
640 | zStr = "gauss76z" |
---|
641 | nord = 76 |
---|
642 | |
---|
643 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
644 | Make/D/N=76 $weightStr,$zStr |
---|
645 | Wave w76 = $weightStr |
---|
646 | Wave z76 = $zStr // wave references to pass |
---|
647 | Make76GaussPoints(w76,z76) |
---|
648 | else |
---|
649 | if(exists(weightStr) > 1) |
---|
650 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
651 | endif |
---|
652 | Wave w76 = $weightStr |
---|
653 | Wave z76 = $zStr // create the wave references |
---|
654 | endif |
---|
655 | |
---|
656 | //limits of integration are input to function |
---|
657 | va = loLim |
---|
658 | vb = upLim |
---|
659 | // Using 76 Gauss points |
---|
660 | // remember to index from 0,size-1 |
---|
661 | |
---|
662 | summ = 0.0 // initialize integral |
---|
663 | ii=0 // loop counter |
---|
664 | do |
---|
665 | // calculate Gauss points on integration interval (q-value for evaluation) |
---|
666 | zi = ( z76[ii]*(vb-va) + vb + va )/2.0 |
---|
667 | //calculate partial sum for the passed-in model function |
---|
668 | yyy = w76[ii] * fcn(w,x,zi) |
---|
669 | summ += yyy //add to the running total of the quadrature |
---|
670 | ii+=1 |
---|
671 | while (ii<nord) // end of loop over quadrature points |
---|
672 | |
---|
673 | // calculate value of integral to return |
---|
674 | answer = (vb-va)/2.0*summ |
---|
675 | |
---|
676 | Return (answer) |
---|
677 | End |
---|
678 | /////////////////////////////////////////////////////////////// |
---|
679 | |
---|
680 | //////Resolution Smearing Utilities |
---|
681 | |
---|
682 | // To check for the existence of all waves needed for smearing |
---|
683 | // returns 1 if any waves are missing, 0 if all is OK |
---|
684 | Function ResolutionWavesMissing() |
---|
685 | |
---|
686 | SVAR/Z sq = gSig_Q |
---|
687 | SVAR/Z qb = gQ_bar |
---|
688 | SVAR/Z sh = gShadow |
---|
689 | SVAR/Z gQ = gQVals |
---|
690 | |
---|
691 | Variable err=0 |
---|
692 | if(!SVAR_Exists(sq) || !SVAR_Exists(qb) || !SVAR_Exists(sh) || !SVAR_Exists(gQ)) |
---|
693 | DoAlert 0,"Some 6-column QSIG waves are missing. Re-load experimental data with LoadOneDData macro" |
---|
694 | return(1) |
---|
695 | endif |
---|
696 | |
---|
697 | if(WaveExists($sq) == 0) //wave ref does not exist |
---|
698 | err=1 |
---|
699 | endif |
---|
700 | if(WaveExists($qb) == 0) //wave ref does not exist |
---|
701 | err=1 |
---|
702 | endif |
---|
703 | if(WaveExists($sh) == 0) //wave ref does not exist |
---|
704 | err=1 |
---|
705 | endif |
---|
706 | if(WaveExists($gQ) == 0) //wave ref does not exist |
---|
707 | err=1 |
---|
708 | endif |
---|
709 | |
---|
710 | if(err) |
---|
711 | DoAlert 0,"Some 6-column QSIG waves are missing. Re-load experimental data with 'Load SANS or USANS Data' macro" |
---|
712 | endif |
---|
713 | |
---|
714 | return(err) |
---|
715 | end |
---|
716 | |
---|
717 | //////Resolution Smearing Utilities |
---|
718 | |
---|
719 | // To check for the existence of all waves needed for smearing |
---|
720 | // returns 1 if any waves are missing, 0 if all is OK |
---|
721 | // str passed in is the data folder containing the data |
---|
722 | // |
---|
723 | // 19 JUN07 using new data folder structure for loading |
---|
724 | // and resolution matrix |
---|
725 | Function ResolutionWavesMissingDF(str) |
---|
726 | String str |
---|
727 | |
---|
728 | String DF="root:"+str+":" |
---|
729 | |
---|
730 | WAVE/Z res = $(DF+str+"_res") |
---|
731 | |
---|
732 | if(!WaveExists(res)) |
---|
733 | DoAlert 0,"The resolution matrix is missing. Re-load experimental data with the 'Load SANS or USANS Data' macro" |
---|
734 | return(1) |
---|
735 | endif |
---|
736 | |
---|
737 | return(0) |
---|
738 | end |
---|
739 | |
---|
740 | /////////////////////////////////////////////////////////////// |
---|
741 | |
---|
742 | // "backwards" wrapped to reduce redundant code |
---|
743 | // there are only 4 choices of N (5,10,20,76) for smearing |
---|
744 | Function Smear_Model_N(fcn,w,x,resW,wi,zi,nord) |
---|
745 | FUNCREF SANSModelAAO_proto fcn |
---|
746 | Wave w //coefficients of function fcn(w,x) |
---|
747 | Variable x //x-value (q) for the calculation THIS IS PASSED IN AS A WAVE |
---|
748 | Wave resW // Nx4 or NxN matrix of resolution |
---|
749 | Wave wi //weight wave |
---|
750 | Wave zi //abscissa wave |
---|
751 | Variable nord //order of integration |
---|
752 | |
---|
753 | NVAR dQv = root:Packages:NIST:USANS_dQv |
---|
754 | |
---|
755 | // local variables |
---|
756 | Variable ii,va,vb |
---|
757 | Variable answer,i_shad,i_qbar,i_sigq |
---|
758 | |
---|
759 | // current x point is the q-value for evaluation |
---|
760 | // |
---|
761 | // * for the input x, the resolution function waves are interpolated to get the correct values for |
---|
762 | // sigq, qbar and shad - since the model x-spacing may not be the same as |
---|
763 | // the experimental QSIG data. This is always the case when curve fitting, since fit_wave is |
---|
764 | // Igor-defined as 200 points and has its own (linear) q-(x)-scaling which will be quite different |
---|
765 | // from experimental data. |
---|
766 | // **note** if the (x) passed in is the experimental q-values, these values are |
---|
767 | // returned from the interpolation (as expected) |
---|
768 | |
---|
769 | Make/O/D/N=(DimSize(resW, 0)) sigQ,qbar,shad,qvals |
---|
770 | sigq = resW[p][0] //std dev of resolution fn |
---|
771 | qbar = resW[p][1] //mean q-value |
---|
772 | shad = resW[p][2] //beamstop shadow factor |
---|
773 | qvals = resW[p][3] //q-values where R(q) is known |
---|
774 | |
---|
775 | i_shad = interp(x,qvals,shad) |
---|
776 | i_qbar = interp(x,qvals,qbar) |
---|
777 | i_sigq = interp(x,qvals,sigq) |
---|
778 | |
---|
779 | // set up the integration |
---|
780 | // number of Gauss Quadrature points |
---|
781 | |
---|
782 | if (isSANSResolution(i_sigq)) |
---|
783 | |
---|
784 | // end points of integration |
---|
785 | // limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero |
---|
786 | // +/- 3 sigq catches 99.73% of distrubution |
---|
787 | // change limits (and spacing of zi) at each evaluation based on R() |
---|
788 | //integration from va to vb |
---|
789 | |
---|
790 | va = -3*i_sigq + i_qbar |
---|
791 | if (va<0) |
---|
792 | va=0 //to avoid numerical error when va<0 (-ve q-value) |
---|
793 | // Print "truncated Gaussian at nominal q = ",x |
---|
794 | endif |
---|
795 | vb = 3*i_sigq + i_qbar |
---|
796 | |
---|
797 | // Using 20 Gauss points |
---|
798 | ii=0 // loop counter |
---|
799 | // do the calculation as a single pass w/AAO function |
---|
800 | Make/O/D/N=(nord) Resoln,yyy,xGauss |
---|
801 | do |
---|
802 | // calculate Gauss points on integration interval (q-value for evaluation) |
---|
803 | xGauss[ii] = ( zi[ii]*(vb-va) + vb + va )/2.0 |
---|
804 | // calculate resolution function at input q-value (use the interpolated values and zi) |
---|
805 | Resoln[ii] = i_shad/sqrt(2*pi*i_sigq*i_sigq) |
---|
806 | Resoln[ii] *= exp((-1*(xGauss[ii] - i_qbar)^2)/(2*i_sigq*i_sigq)) |
---|
807 | ii+=1 |
---|
808 | while (ii<nord) // end of loop over quadrature points |
---|
809 | |
---|
810 | fcn(w,yyy,xGauss) //yyy is the return value as a wave |
---|
811 | |
---|
812 | yyy *= wi *Resoln //multiply function by resolution and weights |
---|
813 | // calculate value of integral to return |
---|
814 | answer = (vb-va)/2.0*sum(yyy) |
---|
815 | // all scaling, background addition... etc. is done in the model calculation |
---|
816 | |
---|
817 | else |
---|
818 | //smear with the USANS routine |
---|
819 | // Make global string and local variables |
---|
820 | // now data folder aware, necessary for GlobalFit = FULL path to wave |
---|
821 | String/G gTrap_coefStr = GetWavesDataFolder(w, 2 ) |
---|
822 | Variable maxiter=20, tol=1e-4 |
---|
823 | |
---|
824 | // set up limits for the integration |
---|
825 | va=0 |
---|
826 | vb=abs(dQv) |
---|
827 | |
---|
828 | Variable/G gEvalQval = x |
---|
829 | |
---|
830 | // call qtrap to do actual work |
---|
831 | answer = qtrap_USANS(fcn,va,vb,tol,maxiter) |
---|
832 | answer /= vb |
---|
833 | |
---|
834 | endif |
---|
835 | |
---|
836 | //killing these waves is cleaner, but MUCH SLOWER |
---|
837 | // Killwaves/Z Resoln,yyy,xGauss |
---|
838 | // Killwaves/Z sigQ,qbar,shad,qvals |
---|
839 | Return (answer) |
---|
840 | |
---|
841 | End |
---|
842 | |
---|
843 | //resolution smearing, using only 5 Gauss points |
---|
844 | // |
---|
845 | // |
---|
846 | Function Smear_Model_5(fcn,w,x,answer,resW) |
---|
847 | FUNCREF SANSModelAAO_proto fcn |
---|
848 | Wave w //coefficients of function fcn(w,x) |
---|
849 | Wave x //x-value (q) for the calculation |
---|
850 | Wave answer // ywave for calculation result |
---|
851 | Wave resW // Nx4 or NxN matrix of resolution |
---|
852 | NVAR useTrap = root:Packages:NIST:USANSUseTrap |
---|
853 | |
---|
854 | String weightStr,zStr |
---|
855 | Variable nord=5 |
---|
856 | |
---|
857 | if (dimsize(resW,1) > 4 && useTrap !=1) |
---|
858 | if(dimsize(resW,1) != dimsize(answer,0) ) |
---|
859 | Abort "ResW and answer are different dimensions - (res,ans)"+num2str(dimsize(resW,1))+","+num2str(dimsize(answer,0)) |
---|
860 | endif |
---|
861 | //USANS Weighting matrix is present. |
---|
862 | fcn(w,answer,x) |
---|
863 | |
---|
864 | MatrixOP/O answer = resW x answer |
---|
865 | //Duplicate/O answer,tmpMat |
---|
866 | //MatrixOP/O answer = resW x tmpMat |
---|
867 | Return(0) |
---|
868 | else |
---|
869 | weightStr = "gauss5wt" |
---|
870 | zStr = "gauss5z" |
---|
871 | |
---|
872 | // if wt,z waves don't exist, create them (only check for weight, should really check for both) |
---|
873 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
874 | Make/D/N=(nord) $weightStr,$zStr |
---|
875 | Wave weightW = $weightStr |
---|
876 | Wave abscissW = $zStr // wave references to pass |
---|
877 | Make5GaussPoints(weightW,abscissW) |
---|
878 | else |
---|
879 | if(exists(weightStr) > 1) |
---|
880 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
881 | endif |
---|
882 | Wave weightW = $weightStr |
---|
883 | Wave abscissW = $zStr // create the wave references |
---|
884 | endif |
---|
885 | |
---|
886 | answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord) |
---|
887 | Return (0) |
---|
888 | endif |
---|
889 | |
---|
890 | End |
---|
891 | |
---|
892 | //resolution smearing, using only 10 Gauss points |
---|
893 | // |
---|
894 | // |
---|
895 | Function Smear_Model_10(fcn,w,x,answer,resW) |
---|
896 | FUNCREF SANSModelAAO_proto fcn |
---|
897 | Wave w //coefficients of function fcn(w,x) |
---|
898 | Wave x //x-value (q) for the calculation |
---|
899 | Wave answer // ywave for calculation result |
---|
900 | Wave resW // Nx4 or NxN matrix of resolution |
---|
901 | |
---|
902 | String weightStr,zStr |
---|
903 | Variable nord=10 |
---|
904 | |
---|
905 | if (dimsize(resW,1) > 4) |
---|
906 | if(dimsize(resW,1) != dimsize(answer,0) ) |
---|
907 | Abort "ResW and answer are different dimensions - (res,ans)"+num2str(dimsize(resW,1))+","+num2str(dimsize(answer,0)) |
---|
908 | endif |
---|
909 | //USANS Weighting matrix is present. |
---|
910 | fcn(w,answer,x) |
---|
911 | |
---|
912 | MatrixOP/O answer = resW x answer |
---|
913 | //Duplicate/O answer,tmpMat |
---|
914 | //MatrixOP/O answer = resW x tmpMat |
---|
915 | Return(0) |
---|
916 | else |
---|
917 | weightStr = "gauss10wt" |
---|
918 | zStr = "gauss10z" |
---|
919 | |
---|
920 | // if wt,z waves don't exist, create them (only check for weight, should really check for both) |
---|
921 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
922 | Make/D/N=(nord) $weightStr,$zStr |
---|
923 | Wave weightW = $weightStr |
---|
924 | Wave abscissW = $zStr // wave references to pass |
---|
925 | Make10GaussPoints(weightW,abscissW) |
---|
926 | else |
---|
927 | if(exists(weightStr) > 1) |
---|
928 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
929 | endif |
---|
930 | Wave weightW = $weightStr |
---|
931 | Wave abscissW = $zStr // create the wave references |
---|
932 | endif |
---|
933 | |
---|
934 | answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord) |
---|
935 | Return (0) |
---|
936 | endif |
---|
937 | |
---|
938 | End |
---|
939 | |
---|
940 | // |
---|
941 | //Smear_Model_20(SphereForm,s.coefW,s.yW,s.xW,s.resW) |
---|
942 | // |
---|
943 | // Wave sigq //std dev of resolution fn |
---|
944 | // Wave qbar //mean q-value |
---|
945 | // Wave shad //beamstop shadow factor |
---|
946 | // Wave qvals //q-values where R(q) is known |
---|
947 | // |
---|
948 | Function Smear_Model_20(fcn,w,x,answer,resW) |
---|
949 | FUNCREF SANSModelAAO_proto fcn |
---|
950 | Wave w //coefficients of function fcn(w,x) |
---|
951 | Wave x //x-value (q) for the calculation |
---|
952 | Wave answer // ywave for calculation result |
---|
953 | Wave resW // Nx4 or NxN matrix of resolution |
---|
954 | NVAR useTrap = root:Packages:NIST:USANSUseTrap |
---|
955 | |
---|
956 | String weightStr,zStr |
---|
957 | Variable nord=20 |
---|
958 | |
---|
959 | if (dimsize(resW,1) > 4 && useTrap != 1) |
---|
960 | if(dimsize(resW,1) != dimsize(answer,0) ) |
---|
961 | Abort "ResW and answer are different dimensions - (res,ans)"+num2str(dimsize(resW,1))+","+num2str(dimsize(answer,0)) |
---|
962 | endif |
---|
963 | //USANS Weighting matrix is present. |
---|
964 | fcn(w,answer,x) |
---|
965 | |
---|
966 | MatrixOP/O answer = resW x answer |
---|
967 | //Duplicate/O answer,tmpMat |
---|
968 | //MatrixOP/O answer = resW x tmpMat |
---|
969 | Return(0) |
---|
970 | else |
---|
971 | weightStr = "gauss20wt" |
---|
972 | zStr = "gauss20z" |
---|
973 | |
---|
974 | // if wt,z waves don't exist, create them (only check for weight, should really check for both) |
---|
975 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
976 | Make/D/N=(nord) $weightStr,$zStr |
---|
977 | Wave weightW = $weightStr |
---|
978 | Wave abscissW = $zStr // wave references to pass |
---|
979 | Make20GaussPoints(weightW,abscissW) |
---|
980 | else |
---|
981 | if(exists(weightStr) > 1) |
---|
982 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
983 | endif |
---|
984 | Wave weightW = $weightStr |
---|
985 | Wave abscissW = $zStr // create the wave references |
---|
986 | endif |
---|
987 | |
---|
988 | answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord) |
---|
989 | Return (0) |
---|
990 | endif |
---|
991 | |
---|
992 | End |
---|
993 | /////////////////////////////////////////////////////////////// |
---|
994 | Function Smear_Model_76(fcn,w,x,answer,resW) |
---|
995 | FUNCREF SANSModelAAO_proto fcn |
---|
996 | Wave w //coefficients of function fcn(w,x) |
---|
997 | Wave x //x-value (q) for the calculation |
---|
998 | Wave answer // ywave for calculation result |
---|
999 | Wave resW // Nx4 or NxN matrix of resolution |
---|
1000 | NVAR useTrap = root:Packages:NIST:USANSUseTrap |
---|
1001 | |
---|
1002 | String weightStr,zStr |
---|
1003 | Variable nord=76 |
---|
1004 | |
---|
1005 | if (dimsize(resW,1) > 4 && useTrap != 1) |
---|
1006 | if(dimsize(resW,1) != dimsize(answer,0) ) |
---|
1007 | Abort "ResW and answer are different dimensions - (res,ans)"+num2str(dimsize(resW,1))+","+num2str(dimsize(answer,0)) |
---|
1008 | endif |
---|
1009 | //USANS Weighting matrix is present. |
---|
1010 | fcn(w,answer,x) |
---|
1011 | |
---|
1012 | MatrixOP/O answer = resW x answer |
---|
1013 | //Duplicate/O answer,tmpMat |
---|
1014 | //MatrixOP/O answer = resW x tmpMat |
---|
1015 | Return(0) |
---|
1016 | else |
---|
1017 | weightStr = "gauss76wt" |
---|
1018 | zStr = "gauss76z" |
---|
1019 | |
---|
1020 | // if wt,z waves don't exist, create them (only check for weight, should really check for both) |
---|
1021 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
1022 | Make/D/N=(nord) $weightStr,$zStr |
---|
1023 | Wave weightW = $weightStr |
---|
1024 | Wave abscissW = $zStr // wave references to pass |
---|
1025 | Make76GaussPoints(weightW,abscissW) |
---|
1026 | else |
---|
1027 | if(exists(weightStr) > 1) |
---|
1028 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
1029 | endif |
---|
1030 | Wave weightW = $weightStr |
---|
1031 | Wave abscissW = $zStr // create the wave references |
---|
1032 | endif |
---|
1033 | |
---|
1034 | answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord) |
---|
1035 | Return (0) |
---|
1036 | endif |
---|
1037 | |
---|
1038 | End |
---|
1039 | |
---|
1040 | |
---|
1041 | /////////////////////////////////////////////////////////////// |
---|
1042 | |
---|
1043 | //typically, the first point (or any point) of sigQ is passed |
---|
1044 | // if negative, it's USANS data... |
---|
1045 | Function isSANSResolution(val) |
---|
1046 | Variable val |
---|
1047 | if(val>=0) |
---|
1048 | return(1) //true, SANS data |
---|
1049 | else |
---|
1050 | return(0) //false, USANS |
---|
1051 | endif |
---|
1052 | End |
---|
1053 | |
---|
1054 | Function GenericQuadrature_proto(w,x,dum) |
---|
1055 | Wave w |
---|
1056 | Variable x,dum |
---|
1057 | |
---|
1058 | Print "in GenericQuadrature_proto function" |
---|
1059 | return(1) |
---|
1060 | end |
---|
1061 | |
---|
1062 | // prototype function for smearing routines, Smear_Model_N |
---|
1063 | // and trapzd_USANS() and qtrap_USANS() |
---|
1064 | // it intentionally does nothing |
---|
1065 | Function SANSModel_proto(w,x) |
---|
1066 | Wave w |
---|
1067 | Variable x |
---|
1068 | |
---|
1069 | Print "in SANSModel_proto function" |
---|
1070 | return(1) |
---|
1071 | end |
---|
1072 | |
---|
1073 | // prototype function for smearing routines, Smear_Model_N |
---|
1074 | // and trapzd_USANS() and qtrap_USANS() |
---|
1075 | // it intentionally does nothing |
---|
1076 | Function SANSModelAAO_proto(w,yw,xw) |
---|
1077 | Wave w,yw,xw |
---|
1078 | |
---|
1079 | Print "in SANSModelAAO_proto function" |
---|
1080 | return(1) |
---|
1081 | end |
---|
1082 | |
---|
1083 | // prototype function for fit wrapper |
---|
1084 | // it intentionally does nothing |
---|
1085 | Function SANSModelSTRUCT_proto(s) |
---|
1086 | Struct ResSmearAAOStruct &s |
---|
1087 | |
---|
1088 | Print "in SANSModelSTRUCT_proto function" |
---|
1089 | return(1) |
---|
1090 | end |
---|
1091 | |
---|
1092 | // prototype function for 2D smearing routine |
---|
1093 | ThreadSafe Function SANS_2D_ModelAAO_proto(w,zw,xw,yw) |
---|
1094 | Wave w,zw,xw,yw |
---|
1095 | |
---|
1096 | Print "in SANSModelAAO_proto function" |
---|
1097 | return(1) |
---|
1098 | end |
---|
1099 | |
---|
1100 | // prototype function for fit wrapper using 2D smearing |
---|
1101 | // not used (yet) |
---|
1102 | Function SANS_2D_ModelSTRUCT_proto(s) |
---|
1103 | Struct ResSmear_2D_AAOStruct &s |
---|
1104 | |
---|
1105 | Print "in SANSModelSTRUCT_proto function" |
---|
1106 | return(1) |
---|
1107 | end |
---|
1108 | |
---|
1109 | //Numerical Recipes routine to calculate the nn(th) stage |
---|
1110 | //refinement of a trapezoid integration |
---|
1111 | // |
---|
1112 | //must be called sequentially from nn=1...n from qtrap() |
---|
1113 | // to cumulatively refine the integration value |
---|
1114 | // |
---|
1115 | // in the conversion: |
---|
1116 | // -- s was replaced with sVal and declared global (rather than static) |
---|
1117 | // so that the nn-1 value would be available during the nn(th) call |
---|
1118 | // |
---|
1119 | // -- the specific coefficient wave for func() is passed in as a |
---|
1120 | // global string (then converted to a wave reference), since |
---|
1121 | // func() will eventually call sphereForm() |
---|
1122 | // |
---|
1123 | Function trapzd_USANS(fcn,aa,bb,nn) |
---|
1124 | FUNCREF SANSModelAAO_proto fcn |
---|
1125 | Variable aa,bb,nn |
---|
1126 | |
---|
1127 | Variable xx,tnm,summ,del |
---|
1128 | Variable it,jj,arg1,arg2 |
---|
1129 | NVAR sVal=sVal //calling function must initialize this global |
---|
1130 | NVAR qval = gEvalQval |
---|
1131 | SVAR cwStr = gTrap_CoefStr //pass in the coefficient wave (string) |
---|
1132 | Wave cw=$cwStr |
---|
1133 | Variable temp=0 |
---|
1134 | Make/D/O/N=2 tmp_xw,tmp_yw |
---|
1135 | if(nn==1) |
---|
1136 | tmp_xw[0] = sqrt(qval^2 + aa^2) |
---|
1137 | tmp_xw[1] = sqrt(qval^2 + bb^2) |
---|
1138 | fcn(cw,tmp_yw,tmp_xw) |
---|
1139 | temp = 0.5*(bb-aa)*(tmp_yw[0] + tmp_yw[1]) |
---|
1140 | sval = temp |
---|
1141 | return(sVal) |
---|
1142 | else |
---|
1143 | it=1 |
---|
1144 | it= 2^(nn-2) //done in NR with a bit shift <<= |
---|
1145 | tnm = it |
---|
1146 | del = (bb - aa)/tnm //this is the spacing of points to add |
---|
1147 | xx = aa+0.5*del |
---|
1148 | summ=0 |
---|
1149 | for(jj=1;jj<=it;jj+=1) |
---|
1150 | tmp_xw = sqrt(qval^2 + xx^2) |
---|
1151 | fcn(cw,tmp_yw,tmp_xw) //not the most efficient... but replaced by the matrix method |
---|
1152 | summ += tmp_yw[0] |
---|
1153 | xx += del |
---|
1154 | endfor |
---|
1155 | sval = 0.5*(sval+(bb-aa)*summ/tnm) //replaces sval with its refined value |
---|
1156 | return (sval) |
---|
1157 | endif |
---|
1158 | //KillWaves/Z tmp_xw,tmp_yw |
---|
1159 | End |
---|
1160 | |
---|
1161 | ////Numerical Recipes routine to calculate the nn(th) stage |
---|
1162 | ////refinement of a trapezoid integration |
---|
1163 | //// |
---|
1164 | ////must be called sequentially from nn=1...n from qtrap() |
---|
1165 | //// to cumulatively refine the integration value |
---|
1166 | //// |
---|
1167 | //// in the conversion: |
---|
1168 | //// -- s was replaced with sVal and declared global (rather than static) |
---|
1169 | //// so that the nn-1 value would be available during the nn(th) call |
---|
1170 | //// |
---|
1171 | //// -- the specific coefficient wave for func() is passed in as a |
---|
1172 | //// global string (then converted to a wave reference), since |
---|
1173 | //// func() will eventually call sphereForm() |
---|
1174 | //// |
---|
1175 | //Function trapzd_USANS_point(fcn,aa,bb,nn) |
---|
1176 | // FUNCREF SANSModel_proto fcn |
---|
1177 | // Variable aa,bb,nn |
---|
1178 | // |
---|
1179 | // Variable xx,tnm,summ,del |
---|
1180 | // Variable it,jj,arg1,arg2 |
---|
1181 | // NVAR sVal=sVal //calling function must initialize this global |
---|
1182 | // NVAR qval = gEvalQval |
---|
1183 | // SVAR cwStr = gTrap_CoefStr //pass in the coefficient wave (string) |
---|
1184 | // Wave cw=$cwStr |
---|
1185 | // Variable temp=0 |
---|
1186 | // if(nn==1) |
---|
1187 | // arg1 = qval^2 + aa^2 |
---|
1188 | // arg2 = qval^2 + bb^2 |
---|
1189 | // temp = 0.5*(bb-aa)*(fcn(cw,sqrt(arg1)) + fcn(cw,sqrt(arg2))) |
---|
1190 | // sval = temp |
---|
1191 | // return(sVal) |
---|
1192 | // else |
---|
1193 | // it=1 |
---|
1194 | // it= 2^(nn-2) //done in NR with a bit shift <<= |
---|
1195 | // tnm = it |
---|
1196 | // del = (bb - aa)/tnm //this is the spacing of points to add |
---|
1197 | // xx = aa+0.5*del |
---|
1198 | // summ=0 |
---|
1199 | // for(jj=1;jj<=it;jj+=1) |
---|
1200 | // arg1 = qval^2 + xx^2 |
---|
1201 | // summ += fcn(cw,sqrt(arg1)) |
---|
1202 | // xx += del |
---|
1203 | // endfor |
---|
1204 | // sval = 0.5*(sval+(bb-aa)*summ/tnm) //replaces sval with its refined value |
---|
1205 | // return (sval) |
---|
1206 | // endif |
---|
1207 | // |
---|
1208 | //End |
---|
1209 | |
---|
1210 | // Numerical Recipes routine to calculate the integral of a |
---|
1211 | // specified function, trapezoid rule is used to a user-specified |
---|
1212 | // level of refinement using sequential calls to trapzd() |
---|
1213 | // |
---|
1214 | // in NR, eps and maxIt were global, pass them in here... |
---|
1215 | // eps typically 1e-5 |
---|
1216 | // maxit typically 20 |
---|
1217 | Function qtrap_USANS(fcn,aa,bb,eps,maxIt) |
---|
1218 | FUNCREF SANSModelAAO_proto fcn |
---|
1219 | Variable aa,bb,eps,maxit |
---|
1220 | |
---|
1221 | Variable/G sVal=0 //create and initialize what trapzd will return |
---|
1222 | Variable jj,ss,olds |
---|
1223 | |
---|
1224 | olds = -1e30 //any number that is not the avg. of endpoints of the funciton |
---|
1225 | for(jj=1;jj<=maxit;jj+=1) //call trapzd() repeatedly until convergence |
---|
1226 | ss = trapzd_USANS(fcn,aa,bb,jj) |
---|
1227 | if( abs(ss-olds) < eps*abs(olds) ) // good enough, stop now |
---|
1228 | return ss |
---|
1229 | endif |
---|
1230 | if( (ss == 0.0) && (olds == 0.0) && (jj>6) ) //no progress? |
---|
1231 | return ss |
---|
1232 | endif |
---|
1233 | olds = ss |
---|
1234 | endfor |
---|
1235 | |
---|
1236 | Print "Maxit exceeded in qtrap. If you're here, there was an error in qtrap" |
---|
1237 | return(ss) //should never get here if function is well-behaved |
---|
1238 | |
---|
1239 | End |
---|
1240 | |
---|
1241 | //// Numerical Recipes routine to calculate the integral of a |
---|
1242 | //// specified function, trapezoid rule is used to a user-specified |
---|
1243 | //// level of refinement using sequential calls to trapzd() |
---|
1244 | //// |
---|
1245 | //// in NR, eps and maxIt were global, pass them in here... |
---|
1246 | //// eps typically 1e-5 |
---|
1247 | //// maxit typically 20 |
---|
1248 | //Function qtrap_USANS_point(fcn,aa,bb,eps,maxIt) |
---|
1249 | // FUNCREF SANSModel_proto fcn |
---|
1250 | // Variable aa,bb,eps,maxit |
---|
1251 | // |
---|
1252 | // Variable/G sVal=0 //create and initialize what trapzd will return |
---|
1253 | // Variable jj,ss,olds |
---|
1254 | // |
---|
1255 | // olds = -1e30 //any number that is not the avg. of endpoints of the funciton |
---|
1256 | // for(jj=1;jj<=maxit;jj+=1) //call trapzd() repeatedly until convergence |
---|
1257 | // ss = trapzd_USANS_point(fcn,aa,bb,jj) |
---|
1258 | // if( abs(ss-olds) < eps*abs(olds) ) // good enough, stop now |
---|
1259 | // return ss |
---|
1260 | // endif |
---|
1261 | // if( (ss == 0.0) && (olds == 0.0) && (jj>6) ) //no progress? |
---|
1262 | // return ss |
---|
1263 | // endif |
---|
1264 | // olds = ss |
---|
1265 | // endfor |
---|
1266 | // |
---|
1267 | // Print "Maxit exceeded in qtrap. If you're here, there was an error in qtrap" |
---|
1268 | // return(ss) //should never get here if function is well-behaved |
---|
1269 | // |
---|
1270 | //End |
---|
1271 | |
---|
1272 | Proc RRW() |
---|
1273 | ResetResolutionWaves() |
---|
1274 | End |
---|
1275 | |
---|
1276 | //utility procedures that are currently untied to any actions, although useful... |
---|
1277 | Proc ResetResolutionWaves(str) |
---|
1278 | String Str |
---|
1279 | Prompt str,"Pick the intensity wave with the resolution you want",popup,WaveList("*_i",";","") |
---|
1280 | |
---|
1281 | |
---|
1282 | Abort "This function is not data floder aware and does nothing..." |
---|
1283 | |
---|
1284 | String/G root:gQvals = str[0,strlen(str)-3]+"_q" |
---|
1285 | String/G root:gSig_Q = str[0,strlen(str)-3]+"sq" |
---|
1286 | String/G root:gQ_bar = str[0,strlen(str)-3]+"qb" |
---|
1287 | String/G root:gShadow = str[0,strlen(str)-3]+"fs" |
---|
1288 | |
---|
1289 | //touch everything to make sure that the dependencies are |
---|
1290 | //properly updated - especially the $gQvals reference in the |
---|
1291 | // dependency assignment |
---|
1292 | fKillDependencies("Smear*") |
---|
1293 | |
---|
1294 | //replace the q-values and intensity (so they're the right length) |
---|
1295 | fResetSmearedModels("Smear*",root:gQvals) |
---|
1296 | |
---|
1297 | fRestoreDependencies("Smear*") |
---|
1298 | End |
---|
1299 | |
---|
1300 | // pass "*" as the matchString to do ALL dependent waves |
---|
1301 | // or "abc*" to get just the matching waves |
---|
1302 | // |
---|
1303 | Function fKillDependencies(matchStr) |
---|
1304 | String matchStr |
---|
1305 | |
---|
1306 | String str=WaveList(matchStr, ";", "" ),item,formula |
---|
1307 | Variable ii |
---|
1308 | |
---|
1309 | for(ii=0;ii<ItemsInList(str ,";");ii+=1) |
---|
1310 | item = StringFromList(ii, str ,";") |
---|
1311 | formula = GetFormula($item) |
---|
1312 | if(cmpstr("",formula)!=0) |
---|
1313 | Printf "wave %s had the formula %s removed\r",item,formula |
---|
1314 | Note $item, "FORMULA:"+formula |
---|
1315 | SetFormula $item, "" //clears the formula |
---|
1316 | endif |
---|
1317 | endfor |
---|
1318 | return(0) |
---|
1319 | end |
---|
1320 | |
---|
1321 | // pass "*" as the matchString to do ALL dependent waves |
---|
1322 | // or "abc*" to get just the matching waves |
---|
1323 | // |
---|
1324 | Function fRestoreDependencies(matchStr) |
---|
1325 | String matchStr |
---|
1326 | |
---|
1327 | String str=WaveList(matchStr, ";", "" ),item,formula |
---|
1328 | Variable ii |
---|
1329 | |
---|
1330 | for(ii=0;ii<ItemsInList(str ,";");ii+=1) |
---|
1331 | item = StringFromList(ii, str ,";") |
---|
1332 | formula = StringByKey("FORMULA", note($item),":",";") |
---|
1333 | if(cmpstr("",formula)!=0) |
---|
1334 | Printf "wave %s had the formula %s restored\r",item,formula |
---|
1335 | Note/K $item |
---|
1336 | SetFormula $item, formula //restores the formula |
---|
1337 | endif |
---|
1338 | endfor |
---|
1339 | return(0) |
---|
1340 | end |
---|
1341 | |
---|
1342 | Function fResetSmearedModels(matchStr,qStr) |
---|
1343 | String matchStr,qStr |
---|
1344 | |
---|
1345 | Duplicate/O $qStr root:smeared_qvals |
---|
1346 | |
---|
1347 | String str=WaveList(matchStr, ";", "" ),item,formula |
---|
1348 | Variable ii |
---|
1349 | |
---|
1350 | for(ii=0;ii<ItemsInList(str ,";");ii+=1) |
---|
1351 | item = StringFromList(ii, str ,";") |
---|
1352 | formula = StringByKey("FORMULA", note($item),":",";") |
---|
1353 | if(cmpstr("",formula)!=0) |
---|
1354 | Printf "wave %s has been duplicated to gQvals\r",item |
---|
1355 | Duplicate/O $qStr $item |
---|
1356 | Note $item, "FORMULA:"+formula //be sure to keep the formula note |
---|
1357 | Print " and the formula is",formula |
---|
1358 | endif |
---|
1359 | endfor |
---|
1360 | return(0) |
---|
1361 | end |
---|
1362 | |
---|
1363 | |
---|
1364 | //// |
---|
1365 | //// moved from RawWindowHook to here - where the Q calculations are available to |
---|
1366 | // reduction and analysis |
---|
1367 | // |
---|
1368 | |
---|
1369 | //phi is defined from +x axis, proceeding CCW around [0,2Pi] |
---|
1370 | Threadsafe Function FindPhi(vx,vy) |
---|
1371 | variable vx,vy |
---|
1372 | |
---|
1373 | variable phi |
---|
1374 | |
---|
1375 | phi = atan(vy/vx) //returns a value from -pi/2 to pi/2 |
---|
1376 | |
---|
1377 | // special cases |
---|
1378 | if(vx==0 && vy > 0) |
---|
1379 | return(pi/2) |
---|
1380 | endif |
---|
1381 | if(vx==0 && vy < 0) |
---|
1382 | return(3*pi/2) |
---|
1383 | endif |
---|
1384 | if(vx >= 0 && vy == 0) |
---|
1385 | return(0) |
---|
1386 | endif |
---|
1387 | if(vx < 0 && vy == 0) |
---|
1388 | return(pi) |
---|
1389 | endif |
---|
1390 | |
---|
1391 | |
---|
1392 | if(vx > 0 && vy > 0) |
---|
1393 | return(phi) |
---|
1394 | endif |
---|
1395 | if(vx < 0 && vy > 0) |
---|
1396 | return(phi + pi) |
---|
1397 | endif |
---|
1398 | if(vx < 0 && vy < 0) |
---|
1399 | return(phi + pi) |
---|
1400 | endif |
---|
1401 | if( vx > 0 && vy < 0) |
---|
1402 | return(phi + 2*pi) |
---|
1403 | endif |
---|
1404 | |
---|
1405 | return(phi) |
---|
1406 | end |
---|
1407 | |
---|
1408 | |
---|
1409 | //function to calculate the overall q-value, given all of the necesary trig inputs |
---|
1410 | //NOTE: detector locations passed in are pixels = 0.5cm real space on the detector |
---|
1411 | //and are in detector coordinates (1,128) rather than axis values |
---|
1412 | //the pixel locations need not be integers, reals are ok inputs |
---|
1413 | //sdd is in meters |
---|
1414 | //wavelength is in Angstroms |
---|
1415 | // |
---|
1416 | //returned magnitude of Q is in 1/Angstroms |
---|
1417 | // |
---|
1418 | Function CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) |
---|
1419 | Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize |
---|
1420 | |
---|
1421 | Variable dx,dy,qval,two_theta,dist |
---|
1422 | |
---|
1423 | Variable pixSizeX=pixSize |
---|
1424 | Variable pixSizeY=pixSize |
---|
1425 | |
---|
1426 | sdd *=100 //convert to cm |
---|
1427 | dx = (xaxval - xctr)*pixSizeX //delta x in cm |
---|
1428 | dy = (yaxval - yctr)*pixSizeY //delta y in cm |
---|
1429 | dist = sqrt(dx^2 + dy^2) |
---|
1430 | |
---|
1431 | two_theta = atan(dist/sdd) |
---|
1432 | |
---|
1433 | qval = 4*Pi/lam*sin(two_theta/2) |
---|
1434 | |
---|
1435 | return qval |
---|
1436 | End |
---|
1437 | |
---|
1438 | //calculates just the q-value in the x-direction on the detector |
---|
1439 | //input/output is the same as CalcQval() |
---|
1440 | //ALL inputs are in detector coordinates |
---|
1441 | // |
---|
1442 | //NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector |
---|
1443 | //sdd is in meters |
---|
1444 | //wavelength is in Angstroms |
---|
1445 | // |
---|
1446 | // repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst) |
---|
1447 | // now properly accounts for qz |
---|
1448 | // |
---|
1449 | Function CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) |
---|
1450 | Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize |
---|
1451 | |
---|
1452 | Variable qx,qval,phi,dx,dy,dist,two_theta |
---|
1453 | |
---|
1454 | qval = CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) |
---|
1455 | |
---|
1456 | sdd *=100 //convert to cm |
---|
1457 | dx = (xaxval - xctr)*pixSize //delta x in cm |
---|
1458 | dy = (yaxval - yctr)*pixSize //delta y in cm |
---|
1459 | phi = FindPhi(dx,dy) |
---|
1460 | |
---|
1461 | //get scattering angle to project onto flat detector => Qr = qval*cos(theta) |
---|
1462 | dist = sqrt(dx^2 + dy^2) |
---|
1463 | two_theta = atan(dist/sdd) |
---|
1464 | |
---|
1465 | qx = qval*cos(two_theta/2)*cos(phi) |
---|
1466 | |
---|
1467 | return qx |
---|
1468 | End |
---|
1469 | |
---|
1470 | //calculates just the q-value in the y-direction on the detector |
---|
1471 | //input/output is the same as CalcQval() |
---|
1472 | //ALL inputs are in detector coordinates |
---|
1473 | //NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector |
---|
1474 | //sdd is in meters |
---|
1475 | //wavelength is in Angstroms |
---|
1476 | // |
---|
1477 | // repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst) |
---|
1478 | // now properly accounts for qz |
---|
1479 | // |
---|
1480 | Function CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) |
---|
1481 | Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize |
---|
1482 | |
---|
1483 | Variable dy,qval,dx,phi,qy,dist,two_theta |
---|
1484 | |
---|
1485 | qval = CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) |
---|
1486 | |
---|
1487 | sdd *=100 //convert to cm |
---|
1488 | dx = (xaxval - xctr)*pixSize //delta x in cm |
---|
1489 | dy = (yaxval - yctr)*pixSize //delta y in cm |
---|
1490 | phi = FindPhi(dx,dy) |
---|
1491 | |
---|
1492 | //get scattering angle to project onto flat detector => Qr = qval*cos(theta) |
---|
1493 | dist = sqrt(dx^2 + dy^2) |
---|
1494 | two_theta = atan(dist/sdd) |
---|
1495 | |
---|
1496 | qy = qval*cos(two_theta/2)*sin(phi) |
---|
1497 | |
---|
1498 | return qy |
---|
1499 | End |
---|
1500 | |
---|
1501 | //calculates just the z-component of the q-vector, not measured on the detector |
---|
1502 | //input/output is the same as CalcQval() |
---|
1503 | //ALL inputs are in detector coordinates |
---|
1504 | //NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector |
---|
1505 | //sdd is in meters |
---|
1506 | //wavelength is in Angstroms |
---|
1507 | // |
---|
1508 | // not actually used, but here for completeness if anyone asks |
---|
1509 | // |
---|
1510 | Function CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) |
---|
1511 | Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize |
---|
1512 | |
---|
1513 | Variable dy,qval,dx,phi,qz,dist,two_theta |
---|
1514 | |
---|
1515 | qval = CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) |
---|
1516 | |
---|
1517 | sdd *=100 //convert to cm |
---|
1518 | |
---|
1519 | //get scattering angle to project onto flat detector => Qr = qval*cos(theta) |
---|
1520 | dx = (xaxval - xctr)*pixSize //delta x in cm |
---|
1521 | dy = (yaxval - yctr)*pixSize //delta y in cm |
---|
1522 | dist = sqrt(dx^2 + dy^2) |
---|
1523 | two_theta = atan(dist/sdd) |
---|
1524 | |
---|
1525 | qz = qval*sin(two_theta/2) |
---|
1526 | |
---|
1527 | return qz |
---|
1528 | End |
---|
1529 | |
---|
1530 | //for command-line testing, replace the function declaration |
---|
1531 | //Function FindQxQy(qq,phi) |
---|
1532 | // Variable qq,phi |
---|
1533 | // Variable qx,qy |
---|
1534 | // |
---|
1535 | // |
---|
1536 | ThreadSafe Function FindQxQy(qq,phi,qx,qy) |
---|
1537 | Variable qq,phi,&qx,&qy |
---|
1538 | |
---|
1539 | qx = sqrt(qq^2/(1+tan(phi)*tan(phi))) |
---|
1540 | qy = qx*tan(phi) |
---|
1541 | |
---|
1542 | if(phi >= 0 && phi <= pi/2) |
---|
1543 | qx = abs(qx) |
---|
1544 | qy = abs(qy) |
---|
1545 | endif |
---|
1546 | |
---|
1547 | if(phi > pi/2 && phi <= pi) |
---|
1548 | qx = -abs(qx) |
---|
1549 | qy = abs(qy) |
---|
1550 | endif |
---|
1551 | |
---|
1552 | if(phi > pi && phi <= pi*3/2) |
---|
1553 | qx = -abs(qx) |
---|
1554 | qy = -abs(qy) |
---|
1555 | endif |
---|
1556 | |
---|
1557 | if(phi > pi*3/2 && phi < 2*pi) |
---|
1558 | qx = abs(qx) |
---|
1559 | qy = -abs(qy) |
---|
1560 | endif |
---|
1561 | |
---|
1562 | |
---|
1563 | // Print "recalculated qx,qy,q = ",qx,qy,sqrt(qx*qx+qy*qy) |
---|
1564 | |
---|
1565 | return(0) |
---|
1566 | end |
---|
1567 | |
---|