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 | |
---|
422 | // special case of integral limits that are identical |
---|
423 | if(upLim == loLim) |
---|
424 | return( fcn(w,x, loLim)) |
---|
425 | endif |
---|
426 | |
---|
427 | // local variables |
---|
428 | Variable ii,va,vb,summ,yyy,zi |
---|
429 | Variable answer,dum |
---|
430 | String weightStr,zStr |
---|
431 | |
---|
432 | weightStr = "gauss"+num2iStr(nord)+"wt" |
---|
433 | zStr = "gauss"+num2istr(nord)+"z" |
---|
434 | |
---|
435 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
436 | Make/D/N=(nord+1) $weightStr,$zStr |
---|
437 | Wave wt = $weightStr |
---|
438 | Wave xx = $zStr // wave references to pass |
---|
439 | Make_N_GaussPoints(wt,xx) //generates the gauss points and removes the extra point |
---|
440 | else |
---|
441 | if(exists(weightStr) > 1) |
---|
442 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
443 | endif |
---|
444 | Wave wt = $weightStr |
---|
445 | Wave xx = $zStr // create the wave references |
---|
446 | endif |
---|
447 | |
---|
448 | //limits of integration are input to function |
---|
449 | va = loLim |
---|
450 | vb = upLim |
---|
451 | // Using 5 Gauss points |
---|
452 | // remember to index from 0,size-1 |
---|
453 | |
---|
454 | summ = 0.0 // initialize integral |
---|
455 | ii=0 // loop counter |
---|
456 | do |
---|
457 | // calculate Gauss points on integration interval (q-value for evaluation) |
---|
458 | zi = ( xx[ii]*(vb-va) + vb + va )/2.0 |
---|
459 | //calculate partial sum for the passed-in model function |
---|
460 | yyy = wt[ii] * fcn(w,x,zi) |
---|
461 | summ += yyy //add to the running total of the quadrature |
---|
462 | ii+=1 |
---|
463 | while (ii<nord) // end of loop over quadrature points |
---|
464 | |
---|
465 | // calculate value of integral to return |
---|
466 | answer = (vb-va)/2.0*summ |
---|
467 | |
---|
468 | Return (answer) |
---|
469 | End |
---|
470 | |
---|
471 | |
---|
472 | |
---|
473 | //////////// |
---|
474 | Function IntegrateFn5(fcn,loLim,upLim,w,x) |
---|
475 | FUNCREF GenericQuadrature_proto fcn |
---|
476 | Variable loLim,upLim //limits of integration |
---|
477 | Wave w //coefficients of function fcn(w,x) |
---|
478 | Variable x //x-value (q) for the calculation |
---|
479 | |
---|
480 | |
---|
481 | // special case of integral limits that are identical |
---|
482 | if(upLim == loLim) |
---|
483 | return( fcn(w,x, loLim)) |
---|
484 | endif |
---|
485 | |
---|
486 | // local variables |
---|
487 | Variable nord,ii,va,vb,summ,yyy,zi |
---|
488 | Variable answer,dum |
---|
489 | String weightStr,zStr |
---|
490 | |
---|
491 | weightStr = "gauss5wt" |
---|
492 | zStr = "gauss5z" |
---|
493 | nord = 5 |
---|
494 | |
---|
495 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
496 | Make/D/N=5 $weightStr,$zStr |
---|
497 | Wave w5 = $weightStr |
---|
498 | Wave z5 = $zStr // wave references to pass |
---|
499 | Make5GaussPoints(w5,z5) |
---|
500 | else |
---|
501 | if(exists(weightStr) > 1) |
---|
502 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
503 | endif |
---|
504 | Wave w5 = $weightStr |
---|
505 | Wave z5 = $zStr // create the wave references |
---|
506 | endif |
---|
507 | |
---|
508 | //limits of integration are input to function |
---|
509 | va = loLim |
---|
510 | vb = upLim |
---|
511 | // Using 5 Gauss points |
---|
512 | // remember to index from 0,size-1 |
---|
513 | |
---|
514 | summ = 0.0 // initialize integral |
---|
515 | ii=0 // loop counter |
---|
516 | do |
---|
517 | // calculate Gauss points on integration interval (q-value for evaluation) |
---|
518 | zi = ( z5[ii]*(vb-va) + vb + va )/2.0 |
---|
519 | //calculate partial sum for the passed-in model function |
---|
520 | yyy = w5[ii] * fcn(w,x,zi) |
---|
521 | summ += yyy //add to the running total of the quadrature |
---|
522 | ii+=1 |
---|
523 | while (ii<nord) // end of loop over quadrature points |
---|
524 | |
---|
525 | // calculate value of integral to return |
---|
526 | answer = (vb-va)/2.0*summ |
---|
527 | |
---|
528 | Return (answer) |
---|
529 | End |
---|
530 | /////////////////////////////////////////////////////////////// |
---|
531 | |
---|
532 | //////////// |
---|
533 | Function IntegrateFn10(fcn,loLim,upLim,w,x) |
---|
534 | FUNCREF GenericQuadrature_proto fcn |
---|
535 | Variable loLim,upLim //limits of integration |
---|
536 | Wave w //coefficients of function fcn(w,x) |
---|
537 | Variable x //x-value (q) for the calculation |
---|
538 | |
---|
539 | // special case of integral limits that are identical |
---|
540 | if(upLim == loLim) |
---|
541 | return( fcn(w,x, loLim)) |
---|
542 | endif |
---|
543 | |
---|
544 | // local variables |
---|
545 | Variable nord,ii,va,vb,summ,yyy,zi |
---|
546 | Variable answer,dum |
---|
547 | String weightStr,zStr |
---|
548 | |
---|
549 | weightStr = "gauss10wt" |
---|
550 | zStr = "gauss10z" |
---|
551 | nord = 10 |
---|
552 | |
---|
553 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
554 | Make/D/N=10 $weightStr,$zStr |
---|
555 | Wave w10 = $weightStr |
---|
556 | Wave z10 = $zStr // wave references to pass |
---|
557 | Make10GaussPoints(w10,z10) |
---|
558 | else |
---|
559 | if(exists(weightStr) > 1) |
---|
560 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
561 | endif |
---|
562 | Wave w10 = $weightStr |
---|
563 | Wave z10 = $zStr // create the wave references |
---|
564 | endif |
---|
565 | |
---|
566 | //limits of integration are input to function |
---|
567 | va = loLim |
---|
568 | vb = upLim |
---|
569 | // Using 10 Gauss points |
---|
570 | // remember to index from 0,size-1 |
---|
571 | |
---|
572 | summ = 0.0 // initialize integral |
---|
573 | ii=0 // loop counter |
---|
574 | do |
---|
575 | // calculate Gauss points on integration interval (q-value for evaluation) |
---|
576 | zi = ( z10[ii]*(vb-va) + vb + va )/2.0 |
---|
577 | //calculate partial sum for the passed-in model function |
---|
578 | yyy = w10[ii] * fcn(w,x,zi) |
---|
579 | summ += yyy //add to the running total of the quadrature |
---|
580 | ii+=1 |
---|
581 | while (ii<nord) // end of loop over quadrature points |
---|
582 | |
---|
583 | // calculate value of integral to return |
---|
584 | answer = (vb-va)/2.0*summ |
---|
585 | |
---|
586 | Return (answer) |
---|
587 | End |
---|
588 | /////////////////////////////////////////////////////////////// |
---|
589 | |
---|
590 | //////////// |
---|
591 | Function IntegrateFn20(fcn,loLim,upLim,w,x) |
---|
592 | FUNCREF GenericQuadrature_proto fcn |
---|
593 | Variable loLim,upLim //limits of integration |
---|
594 | Wave w //coefficients of function fcn(w,x) |
---|
595 | Variable x //x-value (q) for the calculation |
---|
596 | |
---|
597 | // special case of integral limits that are identical |
---|
598 | if(upLim == loLim) |
---|
599 | return( fcn(w,x, loLim)) |
---|
600 | endif |
---|
601 | |
---|
602 | // local variables |
---|
603 | Variable nord,ii,va,vb,summ,yyy,zi |
---|
604 | Variable answer,dum |
---|
605 | String weightStr,zStr |
---|
606 | |
---|
607 | weightStr = "gauss20wt" |
---|
608 | zStr = "gauss20z" |
---|
609 | nord = 20 |
---|
610 | |
---|
611 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
612 | Make/D/N=20 $weightStr,$zStr |
---|
613 | Wave w20 = $weightStr |
---|
614 | Wave z20 = $zStr // wave references to pass |
---|
615 | Make20GaussPoints(w20,z20) |
---|
616 | else |
---|
617 | if(exists(weightStr) > 1) |
---|
618 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
619 | endif |
---|
620 | Wave w20 = $weightStr |
---|
621 | Wave z20 = $zStr // create the wave references |
---|
622 | endif |
---|
623 | |
---|
624 | //limits of integration are input to function |
---|
625 | va = loLim |
---|
626 | vb = upLim |
---|
627 | // Using 20 Gauss points |
---|
628 | // remember to index from 0,size-1 |
---|
629 | |
---|
630 | summ = 0.0 // initialize integral |
---|
631 | ii=0 // loop counter |
---|
632 | do |
---|
633 | // calculate Gauss points on integration interval (q-value for evaluation) |
---|
634 | zi = ( z20[ii]*(vb-va) + vb + va )/2.0 |
---|
635 | //calculate partial sum for the passed-in model function |
---|
636 | yyy = w20[ii] * fcn(w,x,zi) |
---|
637 | summ += yyy //add to the running total of the quadrature |
---|
638 | ii+=1 |
---|
639 | while (ii<nord) // end of loop over quadrature points |
---|
640 | |
---|
641 | // calculate value of integral to return |
---|
642 | answer = (vb-va)/2.0*summ |
---|
643 | |
---|
644 | Return (answer) |
---|
645 | End |
---|
646 | /////////////////////////////////////////////////////////////// |
---|
647 | |
---|
648 | Function IntegrateFn76(fcn,loLim,upLim,w,x) |
---|
649 | FUNCREF GenericQuadrature_proto fcn |
---|
650 | Variable loLim,upLim //limits of integration |
---|
651 | Wave w //coefficients of function fcn(w,x) |
---|
652 | Variable x //x-value (q) for the calculation |
---|
653 | |
---|
654 | //**** The coefficient wave is passed into this function and straight through to the unsmeared model function |
---|
655 | // special case of integral limits that are identical |
---|
656 | if(upLim == loLim) |
---|
657 | return( fcn(w,x, loLim)) |
---|
658 | endif |
---|
659 | |
---|
660 | // local variables |
---|
661 | Variable nord,ii,va,vb,summ,yyy,zi |
---|
662 | Variable answer,dum |
---|
663 | String weightStr,zStr |
---|
664 | |
---|
665 | weightStr = "gauss76wt" |
---|
666 | zStr = "gauss76z" |
---|
667 | nord = 76 |
---|
668 | |
---|
669 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
670 | Make/D/N=76 $weightStr,$zStr |
---|
671 | Wave w76 = $weightStr |
---|
672 | Wave z76 = $zStr // wave references to pass |
---|
673 | Make76GaussPoints(w76,z76) |
---|
674 | else |
---|
675 | if(exists(weightStr) > 1) |
---|
676 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
677 | endif |
---|
678 | Wave w76 = $weightStr |
---|
679 | Wave z76 = $zStr // create the wave references |
---|
680 | endif |
---|
681 | |
---|
682 | //limits of integration are input to function |
---|
683 | va = loLim |
---|
684 | vb = upLim |
---|
685 | // Using 76 Gauss points |
---|
686 | // remember to index from 0,size-1 |
---|
687 | |
---|
688 | summ = 0.0 // initialize integral |
---|
689 | ii=0 // loop counter |
---|
690 | do |
---|
691 | // calculate Gauss points on integration interval (q-value for evaluation) |
---|
692 | zi = ( z76[ii]*(vb-va) + vb + va )/2.0 |
---|
693 | //calculate partial sum for the passed-in model function |
---|
694 | yyy = w76[ii] * fcn(w,x,zi) |
---|
695 | summ += yyy //add to the running total of the quadrature |
---|
696 | ii+=1 |
---|
697 | while (ii<nord) // end of loop over quadrature points |
---|
698 | |
---|
699 | // calculate value of integral to return |
---|
700 | answer = (vb-va)/2.0*summ |
---|
701 | |
---|
702 | Return (answer) |
---|
703 | End |
---|
704 | /////////////////////////////////////////////////////////////// |
---|
705 | |
---|
706 | //////Resolution Smearing Utilities |
---|
707 | |
---|
708 | // To check for the existence of all waves needed for smearing |
---|
709 | // returns 1 if any waves are missing, 0 if all is OK |
---|
710 | Function ResolutionWavesMissing() |
---|
711 | |
---|
712 | SVAR/Z sq = gSig_Q |
---|
713 | SVAR/Z qb = gQ_bar |
---|
714 | SVAR/Z sh = gShadow |
---|
715 | SVAR/Z gQ = gQVals |
---|
716 | |
---|
717 | Variable err=0 |
---|
718 | if(!SVAR_Exists(sq) || !SVAR_Exists(qb) || !SVAR_Exists(sh) || !SVAR_Exists(gQ)) |
---|
719 | DoAlert 0,"Some 6-column QSIG waves are missing. Re-load experimental data with LoadOneDData macro" |
---|
720 | return(1) |
---|
721 | endif |
---|
722 | |
---|
723 | if(WaveExists($sq) == 0) //wave ref does not exist |
---|
724 | err=1 |
---|
725 | endif |
---|
726 | if(WaveExists($qb) == 0) //wave ref does not exist |
---|
727 | err=1 |
---|
728 | endif |
---|
729 | if(WaveExists($sh) == 0) //wave ref does not exist |
---|
730 | err=1 |
---|
731 | endif |
---|
732 | if(WaveExists($gQ) == 0) //wave ref does not exist |
---|
733 | err=1 |
---|
734 | endif |
---|
735 | |
---|
736 | if(err) |
---|
737 | DoAlert 0,"Some 6-column QSIG waves are missing. Re-load experimental data with 'Load SANS or USANS Data' macro" |
---|
738 | endif |
---|
739 | |
---|
740 | return(err) |
---|
741 | end |
---|
742 | |
---|
743 | //////Resolution Smearing Utilities |
---|
744 | |
---|
745 | // To check for the existence of all waves needed for smearing |
---|
746 | // returns 1 if any waves are missing, 0 if all is OK |
---|
747 | // str passed in is the data folder containing the data |
---|
748 | // |
---|
749 | // 19 JUN07 using new data folder structure for loading |
---|
750 | // and resolution matrix |
---|
751 | Function ResolutionWavesMissingDF(str) |
---|
752 | String str |
---|
753 | |
---|
754 | String DF="root:"+str+":" |
---|
755 | |
---|
756 | WAVE/Z res = $(DF+str+"_res") |
---|
757 | |
---|
758 | if(!WaveExists(res)) |
---|
759 | DoAlert 0,"The resolution matrix is missing. Re-load experimental data with the 'Load SANS or USANS Data' macro" |
---|
760 | return(1) |
---|
761 | endif |
---|
762 | |
---|
763 | return(0) |
---|
764 | end |
---|
765 | |
---|
766 | /////////////////////////////////////////////////////////////// |
---|
767 | |
---|
768 | // "backwards" wrapped to reduce redundant code |
---|
769 | // there are only 4 choices of N (5,10,20,76) for smearing |
---|
770 | // |
---|
771 | // |
---|
772 | // 4 MAR 2011 |
---|
773 | // Note: In John's paper, he integrated the Gaussian to +/- 3 sigma and then renormalized |
---|
774 | // to an integral of 1. This "truncated" gaussian was a somewhat better approximation |
---|
775 | // to the triangular resolution function. Here, I integrate to +/- 3 sigma and |
---|
776 | // now correctly renormalize the integral to 1. Hence the smeared calculation in the past was 0.27% low. |
---|
777 | // Confimation of the integral is easily seen by smearing a constant value. |
---|
778 | // |
---|
779 | // Using 5 quadrature points is not recommended, as it doesn't normalize properly using .9973 |
---|
780 | // -- instead, it normalizes to 1.0084, |
---|
781 | // |
---|
782 | Function Smear_Model_N(fcn,w,x,resW,wi,zi,nord) |
---|
783 | FUNCREF SANSModelAAO_proto fcn |
---|
784 | Wave w //coefficients of function fcn(w,x) |
---|
785 | Variable x //x-value (q) for the calculation THIS IS PASSED IN AS A WAVE |
---|
786 | Wave resW // Nx4 or NxN matrix of resolution |
---|
787 | Wave wi //weight wave |
---|
788 | Wave zi //abscissa wave |
---|
789 | Variable nord //order of integration |
---|
790 | |
---|
791 | NVAR dQv = root:Packages:NIST:USANS_dQv |
---|
792 | |
---|
793 | // local variables |
---|
794 | Variable ii,va,vb |
---|
795 | Variable answer,i_shad,i_qbar,i_sigq,normalize=1 |
---|
796 | |
---|
797 | // current x point is the q-value for evaluation |
---|
798 | // |
---|
799 | // * for the input x, the resolution function waves are interpolated to get the correct values for |
---|
800 | // sigq, qbar and shad - since the model x-spacing may not be the same as |
---|
801 | // the experimental QSIG data. This is always the case when curve fitting, since fit_wave is |
---|
802 | // Igor-defined as 200 points and has its own (linear) q-(x)-scaling which will be quite different |
---|
803 | // from experimental data. |
---|
804 | // **note** if the (x) passed in is the experimental q-values, these values are |
---|
805 | // returned from the interpolation (as expected) |
---|
806 | |
---|
807 | Make/O/D/N=(DimSize(resW, 0)) sigQ,qbar,shad,qvals |
---|
808 | sigq = resW[p][0] //std dev of resolution fn |
---|
809 | qbar = resW[p][1] //mean q-value |
---|
810 | shad = resW[p][2] //beamstop shadow factor |
---|
811 | qvals = resW[p][3] //q-values where R(q) is known |
---|
812 | |
---|
813 | i_shad = interp(x,qvals,shad) |
---|
814 | i_qbar = interp(x,qvals,qbar) |
---|
815 | i_sigq = interp(x,qvals,sigq) |
---|
816 | |
---|
817 | // set up the integration |
---|
818 | // number of Gauss Quadrature points |
---|
819 | |
---|
820 | if (isSANSResolution(i_sigq)) |
---|
821 | |
---|
822 | // end points of integration |
---|
823 | // limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero |
---|
824 | // +/- 3 sigq catches 99.73% of distrubution |
---|
825 | // change limits (and spacing of zi) at each evaluation based on R() |
---|
826 | //integration from va to vb |
---|
827 | |
---|
828 | // for +/- 3 sigma ONLY |
---|
829 | if(nord == 5) |
---|
830 | normalize = 1.0057 //empirical correction, N=5 shouldn't be any different |
---|
831 | else |
---|
832 | normalize = 0.9973 |
---|
833 | endif |
---|
834 | |
---|
835 | va = -3*i_sigq + i_qbar |
---|
836 | if (va<0) |
---|
837 | va=0 //to avoid numerical error when va<0 (-ve q-value) |
---|
838 | // Print "truncated Gaussian at nominal q = ",x |
---|
839 | endif |
---|
840 | vb = 3*i_sigq + i_qbar |
---|
841 | |
---|
842 | |
---|
843 | // Using 20 Gauss points |
---|
844 | ii=0 // loop counter |
---|
845 | // do the calculation as a single pass w/AAO function |
---|
846 | Make/O/D/N=(nord) Resoln,yyy,xGauss |
---|
847 | do |
---|
848 | // calculate Gauss points on integration interval (q-value for evaluation) |
---|
849 | xGauss[ii] = ( zi[ii]*(vb-va) + vb + va )/2.0 |
---|
850 | // calculate resolution function at input q-value (use the interpolated values and zi) |
---|
851 | Resoln[ii] = i_shad/sqrt(2*pi*i_sigq*i_sigq) |
---|
852 | Resoln[ii] *= exp((-1*(xGauss[ii] - i_qbar)^2)/(2*i_sigq*i_sigq)) |
---|
853 | ii+=1 |
---|
854 | while (ii<nord) // end of loop over quadrature points |
---|
855 | |
---|
856 | fcn(w,yyy,xGauss) //yyy is the return value as a wave |
---|
857 | |
---|
858 | yyy *= wi *Resoln //multiply function by resolution and weights |
---|
859 | // calculate value of integral to return |
---|
860 | answer = (vb-va)/2.0*sum(yyy) |
---|
861 | // all scaling, background addition... etc. is done in the model calculation |
---|
862 | |
---|
863 | // renormalize to 1 |
---|
864 | answer /= normalize |
---|
865 | else |
---|
866 | //smear with the USANS routine |
---|
867 | // Make global string and local variables |
---|
868 | // now data folder aware, necessary for GlobalFit = FULL path to wave |
---|
869 | String/G gTrap_coefStr = GetWavesDataFolder(w, 2 ) |
---|
870 | Variable maxiter=20, tol=1e-4 |
---|
871 | |
---|
872 | // set up limits for the integration |
---|
873 | va=0 |
---|
874 | vb=abs(dQv) |
---|
875 | |
---|
876 | Variable/G gEvalQval = x |
---|
877 | |
---|
878 | // call qtrap to do actual work |
---|
879 | answer = qtrap_USANS(fcn,va,vb,tol,maxiter) |
---|
880 | answer /= vb |
---|
881 | |
---|
882 | endif |
---|
883 | |
---|
884 | //killing these waves is cleaner, but MUCH SLOWER |
---|
885 | // Killwaves/Z Resoln,yyy,xGauss |
---|
886 | // Killwaves/Z sigQ,qbar,shad,qvals |
---|
887 | Return (answer) |
---|
888 | |
---|
889 | End |
---|
890 | |
---|
891 | //resolution smearing, using only 5 Gauss points |
---|
892 | // |
---|
893 | // |
---|
894 | Function Smear_Model_5(fcn,w,x,answer,resW) |
---|
895 | FUNCREF SANSModelAAO_proto fcn |
---|
896 | Wave w //coefficients of function fcn(w,x) |
---|
897 | Wave x //x-value (q) for the calculation |
---|
898 | Wave answer // ywave for calculation result |
---|
899 | Wave resW // Nx4 or NxN matrix of resolution |
---|
900 | NVAR useTrap = root:Packages:NIST:USANSUseTrap |
---|
901 | |
---|
902 | String weightStr,zStr |
---|
903 | Variable nord=5 |
---|
904 | |
---|
905 | if (dimsize(resW,1) > 4 && useTrap !=1) |
---|
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 = "gauss5wt" |
---|
918 | zStr = "gauss5z" |
---|
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 | Make5GaussPoints(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 | Smear_Model_N_AAO(fcn,w,x,resW,weightW,abscissW,nord,answer) |
---|
936 | |
---|
937 | Return (0) |
---|
938 | endif |
---|
939 | |
---|
940 | End |
---|
941 | |
---|
942 | //resolution smearing, using only 10 Gauss points |
---|
943 | // |
---|
944 | // |
---|
945 | Function Smear_Model_10(fcn,w,x,answer,resW) |
---|
946 | FUNCREF SANSModelAAO_proto fcn |
---|
947 | Wave w //coefficients of function fcn(w,x) |
---|
948 | Wave x //x-value (q) for the calculation |
---|
949 | Wave answer // ywave for calculation result |
---|
950 | Wave resW // Nx4 or NxN matrix of resolution |
---|
951 | |
---|
952 | String weightStr,zStr |
---|
953 | Variable nord=10 |
---|
954 | |
---|
955 | if (dimsize(resW,1) > 4) |
---|
956 | if(dimsize(resW,1) != dimsize(answer,0) ) |
---|
957 | Abort "ResW and answer are different dimensions - (res,ans)"+num2str(dimsize(resW,1))+","+num2str(dimsize(answer,0)) |
---|
958 | endif |
---|
959 | //USANS Weighting matrix is present. |
---|
960 | fcn(w,answer,x) |
---|
961 | |
---|
962 | MatrixOP/O answer = resW x answer |
---|
963 | //Duplicate/O answer,tmpMat |
---|
964 | //MatrixOP/O answer = resW x tmpMat |
---|
965 | Return(0) |
---|
966 | else |
---|
967 | weightStr = "gauss10wt" |
---|
968 | zStr = "gauss10z" |
---|
969 | |
---|
970 | // if wt,z waves don't exist, create them (only check for weight, should really check for both) |
---|
971 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
972 | Make/D/N=(nord) $weightStr,$zStr |
---|
973 | Wave weightW = $weightStr |
---|
974 | Wave abscissW = $zStr // wave references to pass |
---|
975 | Make10GaussPoints(weightW,abscissW) |
---|
976 | else |
---|
977 | if(exists(weightStr) > 1) |
---|
978 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
979 | endif |
---|
980 | Wave weightW = $weightStr |
---|
981 | Wave abscissW = $zStr // create the wave references |
---|
982 | endif |
---|
983 | |
---|
984 | // answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord) |
---|
985 | Smear_Model_N_AAO(fcn,w,x,resW,weightW,abscissW,nord,answer) |
---|
986 | |
---|
987 | Return (0) |
---|
988 | endif |
---|
989 | |
---|
990 | End |
---|
991 | |
---|
992 | // |
---|
993 | //Smear_Model_20(SphereForm,s.coefW,s.yW,s.xW,s.resW) |
---|
994 | // |
---|
995 | // Wave sigq //std dev of resolution fn |
---|
996 | // Wave qbar //mean q-value |
---|
997 | // Wave shad //beamstop shadow factor |
---|
998 | // Wave qvals //q-values where R(q) is known |
---|
999 | // |
---|
1000 | Function Smear_Model_20(fcn,w,x,answer,resW) |
---|
1001 | FUNCREF SANSModelAAO_proto fcn |
---|
1002 | Wave w //coefficients of function fcn(w,x) |
---|
1003 | Wave x //x-value (q) for the calculation |
---|
1004 | Wave answer // ywave for calculation result |
---|
1005 | Wave resW // Nx4 or NxN matrix of resolution |
---|
1006 | NVAR useTrap = root:Packages:NIST:USANSUseTrap |
---|
1007 | |
---|
1008 | String weightStr,zStr |
---|
1009 | Variable nord=20 |
---|
1010 | |
---|
1011 | |
---|
1012 | if (dimsize(resW,1) > 4 && useTrap != 1) |
---|
1013 | if(dimsize(resW,1) != dimsize(answer,0) ) |
---|
1014 | Abort "ResW and answer are different dimensions - (res,ans)"+num2str(dimsize(resW,1))+","+num2str(dimsize(answer,0)) |
---|
1015 | endif |
---|
1016 | //USANS Weighting matrix is present. |
---|
1017 | fcn(w,answer,x) |
---|
1018 | |
---|
1019 | MatrixOP/O answer = resW x answer |
---|
1020 | //Duplicate/O answer,tmpMat |
---|
1021 | //MatrixOP/O answer = resW x tmpMat |
---|
1022 | Return(0) |
---|
1023 | else |
---|
1024 | weightStr = "gauss20wt" |
---|
1025 | zStr = "gauss20z" |
---|
1026 | |
---|
1027 | // if wt,z waves don't exist, create them (only check for weight, should really check for both) |
---|
1028 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
1029 | Make/D/N=(nord) $weightStr,$zStr |
---|
1030 | Wave weightW = $weightStr |
---|
1031 | Wave abscissW = $zStr // wave references to pass |
---|
1032 | Make20GaussPoints(weightW,abscissW) |
---|
1033 | else |
---|
1034 | if(exists(weightStr) > 1) |
---|
1035 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
1036 | endif |
---|
1037 | Wave weightW = $weightStr |
---|
1038 | Wave abscissW = $zStr // create the wave references |
---|
1039 | endif |
---|
1040 | |
---|
1041 | // answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord) |
---|
1042 | Smear_Model_N_AAO(fcn,w,x,resW,weightW,abscissW,nord,answer) |
---|
1043 | |
---|
1044 | Return (0) |
---|
1045 | endif |
---|
1046 | |
---|
1047 | End |
---|
1048 | /////////////////////////////////////////////////////////////// |
---|
1049 | Function Smear_Model_76(fcn,w,x,answer,resW) |
---|
1050 | FUNCREF SANSModelAAO_proto fcn |
---|
1051 | Wave w //coefficients of function fcn(w,x) |
---|
1052 | Wave x //x-value (q) for the calculation |
---|
1053 | Wave answer // ywave for calculation result |
---|
1054 | Wave resW // Nx4 or NxN matrix of resolution |
---|
1055 | NVAR useTrap = root:Packages:NIST:USANSUseTrap |
---|
1056 | |
---|
1057 | String weightStr,zStr |
---|
1058 | Variable nord=76 |
---|
1059 | |
---|
1060 | if (dimsize(resW,1) > 4 && useTrap != 1) |
---|
1061 | if(dimsize(resW,1) != dimsize(answer,0) ) |
---|
1062 | Abort "ResW and answer are different dimensions - (res,ans)"+num2str(dimsize(resW,1))+","+num2str(dimsize(answer,0)) |
---|
1063 | endif |
---|
1064 | //USANS Weighting matrix is present. |
---|
1065 | fcn(w,answer,x) |
---|
1066 | |
---|
1067 | MatrixOP/O answer = resW x answer |
---|
1068 | //Duplicate/O answer,tmpMat |
---|
1069 | //MatrixOP/O answer = resW x tmpMat |
---|
1070 | Return(0) |
---|
1071 | else |
---|
1072 | weightStr = "gauss76wt" |
---|
1073 | zStr = "gauss76z" |
---|
1074 | |
---|
1075 | // if wt,z waves don't exist, create them (only check for weight, should really check for both) |
---|
1076 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
1077 | Make/D/N=(nord) $weightStr,$zStr |
---|
1078 | Wave weightW = $weightStr |
---|
1079 | Wave abscissW = $zStr // wave references to pass |
---|
1080 | Make76GaussPoints(weightW,abscissW) |
---|
1081 | else |
---|
1082 | if(exists(weightStr) > 1) |
---|
1083 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
1084 | endif |
---|
1085 | Wave weightW = $weightStr |
---|
1086 | Wave abscissW = $zStr // create the wave references |
---|
1087 | endif |
---|
1088 | |
---|
1089 | // answer = Smear_Model_N(fcn,w,x,resW,weightW,abscissW,nord) |
---|
1090 | Smear_Model_N_AAO(fcn,w,x,resW,weightW,abscissW,nord,answer) |
---|
1091 | |
---|
1092 | Return (0) |
---|
1093 | endif |
---|
1094 | |
---|
1095 | End |
---|
1096 | |
---|
1097 | |
---|
1098 | /////////////////////////////////////////////////////////////// |
---|
1099 | |
---|
1100 | //typically, the first point (or any point) of sigQ is passed |
---|
1101 | // if negative, it's USANS data... |
---|
1102 | Function isSANSResolution(val) |
---|
1103 | Variable val |
---|
1104 | if(val>=0) |
---|
1105 | return(1) //true, SANS data |
---|
1106 | else |
---|
1107 | return(0) //false, USANS |
---|
1108 | endif |
---|
1109 | End |
---|
1110 | |
---|
1111 | Function GenericQuadrature_proto(w,x,dum) |
---|
1112 | Wave w |
---|
1113 | Variable x,dum |
---|
1114 | |
---|
1115 | Print "in GenericQuadrature_proto function" |
---|
1116 | return(1) |
---|
1117 | end |
---|
1118 | |
---|
1119 | // prototype function for smearing routines, Smear_Model_N |
---|
1120 | // and trapzd_USANS() and qtrap_USANS() |
---|
1121 | // it intentionally does nothing |
---|
1122 | Function SANSModel_proto(w,x) |
---|
1123 | Wave w |
---|
1124 | Variable x |
---|
1125 | |
---|
1126 | Print "in SANSModel_proto function" |
---|
1127 | return(1) |
---|
1128 | end |
---|
1129 | |
---|
1130 | // prototype function for smearing routines, Smear_Model_N |
---|
1131 | // and trapzd_USANS() and qtrap_USANS() |
---|
1132 | // it intentionally does nothing |
---|
1133 | Function SANSModelAAO_proto(w,yw,xw) |
---|
1134 | Wave w,yw,xw |
---|
1135 | |
---|
1136 | Print "in SANSModelAAO_proto function" |
---|
1137 | return(1) |
---|
1138 | end |
---|
1139 | |
---|
1140 | // prototype function for fit wrapper |
---|
1141 | // it intentionally does nothing |
---|
1142 | Function SANSModelSTRUCT_proto(s) |
---|
1143 | Struct ResSmearAAOStruct &s |
---|
1144 | |
---|
1145 | Print "in SANSModelSTRUCT_proto function" |
---|
1146 | return(1) |
---|
1147 | end |
---|
1148 | |
---|
1149 | // prototype function for 2D smearing routine |
---|
1150 | ThreadSafe Function SANS_2D_ModelAAO_proto(w,zw,xw,yw) |
---|
1151 | Wave w,zw,xw,yw |
---|
1152 | |
---|
1153 | Print "in SANSModelAAO_proto function" |
---|
1154 | return(1) |
---|
1155 | end |
---|
1156 | |
---|
1157 | // prototype function for fit wrapper using 2D smearing |
---|
1158 | // not used (yet) |
---|
1159 | Function SANS_2D_ModelSTRUCT_proto(s) |
---|
1160 | Struct ResSmear_2D_AAOStruct &s |
---|
1161 | |
---|
1162 | Print "in SANSModelSTRUCT_proto function" |
---|
1163 | return(1) |
---|
1164 | end |
---|
1165 | |
---|
1166 | //Numerical Recipes routine to calculate the nn(th) stage |
---|
1167 | //refinement of a trapezoid integration |
---|
1168 | // |
---|
1169 | //must be called sequentially from nn=1...n from qtrap() |
---|
1170 | // to cumulatively refine the integration value |
---|
1171 | // |
---|
1172 | // in the conversion: |
---|
1173 | // -- s was replaced with sVal and declared global (rather than static) |
---|
1174 | // so that the nn-1 value would be available during the nn(th) call |
---|
1175 | // |
---|
1176 | // -- the specific coefficient wave for func() is passed in as a |
---|
1177 | // global string (then converted to a wave reference), since |
---|
1178 | // func() will eventually call sphereForm() |
---|
1179 | // |
---|
1180 | Function trapzd_USANS(fcn,aa,bb,nn) |
---|
1181 | FUNCREF SANSModelAAO_proto fcn |
---|
1182 | Variable aa,bb,nn |
---|
1183 | |
---|
1184 | Variable xx,tnm,summ,del |
---|
1185 | Variable it,jj,arg1,arg2 |
---|
1186 | NVAR sVal=sVal //calling function must initialize this global |
---|
1187 | NVAR qval = gEvalQval |
---|
1188 | SVAR cwStr = gTrap_CoefStr //pass in the coefficient wave (string) |
---|
1189 | Wave cw=$cwStr |
---|
1190 | Variable temp=0 |
---|
1191 | Make/D/O/N=2 tmp_xw,tmp_yw |
---|
1192 | if(nn==1) |
---|
1193 | tmp_xw[0] = sqrt(qval^2 + aa^2) |
---|
1194 | tmp_xw[1] = sqrt(qval^2 + bb^2) |
---|
1195 | fcn(cw,tmp_yw,tmp_xw) |
---|
1196 | temp = 0.5*(bb-aa)*(tmp_yw[0] + tmp_yw[1]) |
---|
1197 | sval = temp |
---|
1198 | return(sVal) |
---|
1199 | else |
---|
1200 | it=1 |
---|
1201 | it= 2^(nn-2) //done in NR with a bit shift <<= |
---|
1202 | tnm = it |
---|
1203 | del = (bb - aa)/tnm //this is the spacing of points to add |
---|
1204 | xx = aa+0.5*del |
---|
1205 | summ=0 |
---|
1206 | for(jj=1;jj<=it;jj+=1) |
---|
1207 | tmp_xw = sqrt(qval^2 + xx^2) |
---|
1208 | fcn(cw,tmp_yw,tmp_xw) //not the most efficient... but replaced by the matrix method |
---|
1209 | summ += tmp_yw[0] |
---|
1210 | xx += del |
---|
1211 | endfor |
---|
1212 | sval = 0.5*(sval+(bb-aa)*summ/tnm) //replaces sval with its refined value |
---|
1213 | return (sval) |
---|
1214 | endif |
---|
1215 | //KillWaves/Z tmp_xw,tmp_yw |
---|
1216 | End |
---|
1217 | |
---|
1218 | ////Numerical Recipes routine to calculate the nn(th) stage |
---|
1219 | ////refinement of a trapezoid integration |
---|
1220 | //// |
---|
1221 | ////must be called sequentially from nn=1...n from qtrap() |
---|
1222 | //// to cumulatively refine the integration value |
---|
1223 | //// |
---|
1224 | //// in the conversion: |
---|
1225 | //// -- s was replaced with sVal and declared global (rather than static) |
---|
1226 | //// so that the nn-1 value would be available during the nn(th) call |
---|
1227 | //// |
---|
1228 | //// -- the specific coefficient wave for func() is passed in as a |
---|
1229 | //// global string (then converted to a wave reference), since |
---|
1230 | //// func() will eventually call sphereForm() |
---|
1231 | //// |
---|
1232 | //Function trapzd_USANS_point(fcn,aa,bb,nn) |
---|
1233 | // FUNCREF SANSModel_proto fcn |
---|
1234 | // Variable aa,bb,nn |
---|
1235 | // |
---|
1236 | // Variable xx,tnm,summ,del |
---|
1237 | // Variable it,jj,arg1,arg2 |
---|
1238 | // NVAR sVal=sVal //calling function must initialize this global |
---|
1239 | // NVAR qval = gEvalQval |
---|
1240 | // SVAR cwStr = gTrap_CoefStr //pass in the coefficient wave (string) |
---|
1241 | // Wave cw=$cwStr |
---|
1242 | // Variable temp=0 |
---|
1243 | // if(nn==1) |
---|
1244 | // arg1 = qval^2 + aa^2 |
---|
1245 | // arg2 = qval^2 + bb^2 |
---|
1246 | // temp = 0.5*(bb-aa)*(fcn(cw,sqrt(arg1)) + fcn(cw,sqrt(arg2))) |
---|
1247 | // sval = temp |
---|
1248 | // return(sVal) |
---|
1249 | // else |
---|
1250 | // it=1 |
---|
1251 | // it= 2^(nn-2) //done in NR with a bit shift <<= |
---|
1252 | // tnm = it |
---|
1253 | // del = (bb - aa)/tnm //this is the spacing of points to add |
---|
1254 | // xx = aa+0.5*del |
---|
1255 | // summ=0 |
---|
1256 | // for(jj=1;jj<=it;jj+=1) |
---|
1257 | // arg1 = qval^2 + xx^2 |
---|
1258 | // summ += fcn(cw,sqrt(arg1)) |
---|
1259 | // xx += del |
---|
1260 | // endfor |
---|
1261 | // sval = 0.5*(sval+(bb-aa)*summ/tnm) //replaces sval with its refined value |
---|
1262 | // return (sval) |
---|
1263 | // endif |
---|
1264 | // |
---|
1265 | //End |
---|
1266 | |
---|
1267 | // Numerical Recipes routine to calculate the integral of a |
---|
1268 | // specified function, trapezoid rule is used to a user-specified |
---|
1269 | // level of refinement using sequential calls to trapzd() |
---|
1270 | // |
---|
1271 | // in NR, eps and maxIt were global, pass them in here... |
---|
1272 | // eps typically 1e-5 |
---|
1273 | // maxit typically 20 |
---|
1274 | Function qtrap_USANS(fcn,aa,bb,eps,maxIt) |
---|
1275 | FUNCREF SANSModelAAO_proto fcn |
---|
1276 | Variable aa,bb,eps,maxit |
---|
1277 | |
---|
1278 | Variable/G sVal=0 //create and initialize what trapzd will return |
---|
1279 | Variable jj,ss,olds |
---|
1280 | |
---|
1281 | olds = -1e30 //any number that is not the avg. of endpoints of the funciton |
---|
1282 | for(jj=1;jj<=maxit;jj+=1) //call trapzd() repeatedly until convergence |
---|
1283 | ss = trapzd_USANS(fcn,aa,bb,jj) |
---|
1284 | if( abs(ss-olds) < eps*abs(olds) ) // good enough, stop now |
---|
1285 | return ss |
---|
1286 | endif |
---|
1287 | if( (ss == 0.0) && (olds == 0.0) && (jj>6) ) //no progress? |
---|
1288 | return ss |
---|
1289 | endif |
---|
1290 | olds = ss |
---|
1291 | endfor |
---|
1292 | |
---|
1293 | Print "Maxit exceeded in qtrap. If you're here, there was an error in qtrap" |
---|
1294 | return(ss) //should never get here if function is well-behaved |
---|
1295 | |
---|
1296 | End |
---|
1297 | |
---|
1298 | //// Numerical Recipes routine to calculate the integral of a |
---|
1299 | //// specified function, trapezoid rule is used to a user-specified |
---|
1300 | //// level of refinement using sequential calls to trapzd() |
---|
1301 | //// |
---|
1302 | //// in NR, eps and maxIt were global, pass them in here... |
---|
1303 | //// eps typically 1e-5 |
---|
1304 | //// maxit typically 20 |
---|
1305 | //Function qtrap_USANS_point(fcn,aa,bb,eps,maxIt) |
---|
1306 | // FUNCREF SANSModel_proto fcn |
---|
1307 | // Variable aa,bb,eps,maxit |
---|
1308 | // |
---|
1309 | // Variable/G sVal=0 //create and initialize what trapzd will return |
---|
1310 | // Variable jj,ss,olds |
---|
1311 | // |
---|
1312 | // olds = -1e30 //any number that is not the avg. of endpoints of the funciton |
---|
1313 | // for(jj=1;jj<=maxit;jj+=1) //call trapzd() repeatedly until convergence |
---|
1314 | // ss = trapzd_USANS_point(fcn,aa,bb,jj) |
---|
1315 | // if( abs(ss-olds) < eps*abs(olds) ) // good enough, stop now |
---|
1316 | // return ss |
---|
1317 | // endif |
---|
1318 | // if( (ss == 0.0) && (olds == 0.0) && (jj>6) ) //no progress? |
---|
1319 | // return ss |
---|
1320 | // endif |
---|
1321 | // olds = ss |
---|
1322 | // endfor |
---|
1323 | // |
---|
1324 | // Print "Maxit exceeded in qtrap. If you're here, there was an error in qtrap" |
---|
1325 | // return(ss) //should never get here if function is well-behaved |
---|
1326 | // |
---|
1327 | //End |
---|
1328 | |
---|
1329 | Proc RRW() |
---|
1330 | ResetResolutionWaves() |
---|
1331 | End |
---|
1332 | |
---|
1333 | //utility procedures that are currently untied to any actions, although useful... |
---|
1334 | Proc ResetResolutionWaves(str) |
---|
1335 | String Str |
---|
1336 | Prompt str,"Pick the intensity wave with the resolution you want",popup,WaveList("*_i",";","") |
---|
1337 | |
---|
1338 | |
---|
1339 | Abort "This function is not data floder aware and does nothing..." |
---|
1340 | |
---|
1341 | String/G root:gQvals = str[0,strlen(str)-3]+"_q" |
---|
1342 | String/G root:gSig_Q = str[0,strlen(str)-3]+"sq" |
---|
1343 | String/G root:gQ_bar = str[0,strlen(str)-3]+"qb" |
---|
1344 | String/G root:gShadow = str[0,strlen(str)-3]+"fs" |
---|
1345 | |
---|
1346 | //touch everything to make sure that the dependencies are |
---|
1347 | //properly updated - especially the $gQvals reference in the |
---|
1348 | // dependency assignment |
---|
1349 | fKillDependencies("Smear*") |
---|
1350 | |
---|
1351 | //replace the q-values and intensity (so they're the right length) |
---|
1352 | fResetSmearedModels("Smear*",root:gQvals) |
---|
1353 | |
---|
1354 | fRestoreDependencies("Smear*") |
---|
1355 | End |
---|
1356 | |
---|
1357 | // pass "*" as the matchString to do ALL dependent waves |
---|
1358 | // or "abc*" to get just the matching waves |
---|
1359 | // |
---|
1360 | Function fKillDependencies(matchStr) |
---|
1361 | String matchStr |
---|
1362 | |
---|
1363 | String str=WaveList(matchStr, ";", "" ),item,formula |
---|
1364 | Variable ii |
---|
1365 | |
---|
1366 | for(ii=0;ii<ItemsInList(str ,";");ii+=1) |
---|
1367 | item = StringFromList(ii, str ,";") |
---|
1368 | formula = GetFormula($item) |
---|
1369 | if(cmpstr("",formula)!=0) |
---|
1370 | Printf "wave %s had the formula %s removed\r",item,formula |
---|
1371 | Note $item, "FORMULA:"+formula |
---|
1372 | SetFormula $item, "" //clears the formula |
---|
1373 | endif |
---|
1374 | endfor |
---|
1375 | return(0) |
---|
1376 | end |
---|
1377 | |
---|
1378 | // pass "*" as the matchString to do ALL dependent waves |
---|
1379 | // or "abc*" to get just the matching waves |
---|
1380 | // |
---|
1381 | Function fRestoreDependencies(matchStr) |
---|
1382 | String matchStr |
---|
1383 | |
---|
1384 | String str=WaveList(matchStr, ";", "" ),item,formula |
---|
1385 | Variable ii |
---|
1386 | |
---|
1387 | for(ii=0;ii<ItemsInList(str ,";");ii+=1) |
---|
1388 | item = StringFromList(ii, str ,";") |
---|
1389 | formula = StringByKey("FORMULA", note($item),":",";") |
---|
1390 | if(cmpstr("",formula)!=0) |
---|
1391 | Printf "wave %s had the formula %s restored\r",item,formula |
---|
1392 | Note/K $item |
---|
1393 | SetFormula $item, formula //restores the formula |
---|
1394 | endif |
---|
1395 | endfor |
---|
1396 | return(0) |
---|
1397 | end |
---|
1398 | |
---|
1399 | Function fResetSmearedModels(matchStr,qStr) |
---|
1400 | String matchStr,qStr |
---|
1401 | |
---|
1402 | Duplicate/O $qStr root:smeared_qvals |
---|
1403 | |
---|
1404 | String str=WaveList(matchStr, ";", "" ),item,formula |
---|
1405 | Variable ii |
---|
1406 | |
---|
1407 | for(ii=0;ii<ItemsInList(str ,";");ii+=1) |
---|
1408 | item = StringFromList(ii, str ,";") |
---|
1409 | formula = StringByKey("FORMULA", note($item),":",";") |
---|
1410 | if(cmpstr("",formula)!=0) |
---|
1411 | Printf "wave %s has been duplicated to gQvals\r",item |
---|
1412 | Duplicate/O $qStr $item |
---|
1413 | Note $item, "FORMULA:"+formula //be sure to keep the formula note |
---|
1414 | Print " and the formula is",formula |
---|
1415 | endif |
---|
1416 | endfor |
---|
1417 | return(0) |
---|
1418 | end |
---|
1419 | |
---|
1420 | |
---|
1421 | //// |
---|
1422 | //// moved from RawWindowHook to here - where the Q calculations are available to |
---|
1423 | // reduction and analysis |
---|
1424 | // |
---|
1425 | |
---|
1426 | //phi is defined from +x axis, proceeding CCW around [0,2Pi] |
---|
1427 | Threadsafe Function FindPhi(vx,vy) |
---|
1428 | variable vx,vy |
---|
1429 | |
---|
1430 | variable phi |
---|
1431 | |
---|
1432 | phi = atan(vy/vx) //returns a value from -pi/2 to pi/2 |
---|
1433 | |
---|
1434 | // special cases |
---|
1435 | if(vx==0 && vy > 0) |
---|
1436 | return(pi/2) |
---|
1437 | endif |
---|
1438 | if(vx==0 && vy < 0) |
---|
1439 | return(3*pi/2) |
---|
1440 | endif |
---|
1441 | if(vx >= 0 && vy == 0) |
---|
1442 | return(0) |
---|
1443 | endif |
---|
1444 | if(vx < 0 && vy == 0) |
---|
1445 | return(pi) |
---|
1446 | endif |
---|
1447 | |
---|
1448 | |
---|
1449 | if(vx > 0 && vy > 0) |
---|
1450 | return(phi) |
---|
1451 | endif |
---|
1452 | if(vx < 0 && vy > 0) |
---|
1453 | return(phi + pi) |
---|
1454 | endif |
---|
1455 | if(vx < 0 && vy < 0) |
---|
1456 | return(phi + pi) |
---|
1457 | endif |
---|
1458 | if( vx > 0 && vy < 0) |
---|
1459 | return(phi + 2*pi) |
---|
1460 | endif |
---|
1461 | |
---|
1462 | return(phi) |
---|
1463 | end |
---|
1464 | |
---|
1465 | |
---|
1466 | //function to calculate the overall q-value, given all of the necesary trig inputs |
---|
1467 | //NOTE: detector locations passed in are pixels = 0.5cm real space on the detector |
---|
1468 | //and are in detector coordinates (1,128) rather than axis values |
---|
1469 | //the pixel locations need not be integers, reals are ok inputs |
---|
1470 | //sdd is in meters |
---|
1471 | //wavelength is in Angstroms |
---|
1472 | // |
---|
1473 | //returned magnitude of Q is in 1/Angstroms |
---|
1474 | // |
---|
1475 | Function CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) |
---|
1476 | Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize |
---|
1477 | |
---|
1478 | Variable dx,dy,qval,two_theta,dist |
---|
1479 | |
---|
1480 | Variable pixSizeX=pixSize |
---|
1481 | Variable pixSizeY=pixSize |
---|
1482 | |
---|
1483 | sdd *=100 //convert to cm |
---|
1484 | dx = (xaxval - xctr)*pixSizeX //delta x in cm |
---|
1485 | dy = (yaxval - yctr)*pixSizeY //delta y in cm |
---|
1486 | dist = sqrt(dx^2 + dy^2) |
---|
1487 | |
---|
1488 | two_theta = atan(dist/sdd) |
---|
1489 | |
---|
1490 | qval = 4*Pi/lam*sin(two_theta/2) |
---|
1491 | |
---|
1492 | return qval |
---|
1493 | End |
---|
1494 | |
---|
1495 | //calculates just the q-value in the x-direction on the detector |
---|
1496 | //input/output is the same as CalcQval() |
---|
1497 | //ALL inputs are in detector coordinates |
---|
1498 | // |
---|
1499 | //NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector |
---|
1500 | //sdd is in meters |
---|
1501 | //wavelength is in Angstroms |
---|
1502 | // |
---|
1503 | // repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst) |
---|
1504 | // now properly accounts for qz |
---|
1505 | // |
---|
1506 | Function CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) |
---|
1507 | Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize |
---|
1508 | |
---|
1509 | Variable qx,qval,phi,dx,dy,dist,two_theta |
---|
1510 | |
---|
1511 | qval = CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) |
---|
1512 | |
---|
1513 | sdd *=100 //convert to cm |
---|
1514 | dx = (xaxval - xctr)*pixSize //delta x in cm |
---|
1515 | dy = (yaxval - yctr)*pixSize //delta y in cm |
---|
1516 | phi = FindPhi(dx,dy) |
---|
1517 | |
---|
1518 | //get scattering angle to project onto flat detector => Qr = qval*cos(theta) |
---|
1519 | dist = sqrt(dx^2 + dy^2) |
---|
1520 | two_theta = atan(dist/sdd) |
---|
1521 | |
---|
1522 | qx = qval*cos(two_theta/2)*cos(phi) |
---|
1523 | |
---|
1524 | return qx |
---|
1525 | End |
---|
1526 | |
---|
1527 | //calculates just the q-value in the y-direction on the detector |
---|
1528 | //input/output is the same as CalcQval() |
---|
1529 | //ALL inputs are in detector coordinates |
---|
1530 | //NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector |
---|
1531 | //sdd is in meters |
---|
1532 | //wavelength is in Angstroms |
---|
1533 | // |
---|
1534 | // repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst) |
---|
1535 | // now properly accounts for qz |
---|
1536 | // |
---|
1537 | Function CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) |
---|
1538 | Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize |
---|
1539 | |
---|
1540 | Variable dy,qval,dx,phi,qy,dist,two_theta |
---|
1541 | |
---|
1542 | qval = CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) |
---|
1543 | |
---|
1544 | sdd *=100 //convert to cm |
---|
1545 | dx = (xaxval - xctr)*pixSize //delta x in cm |
---|
1546 | dy = (yaxval - yctr)*pixSize //delta y in cm |
---|
1547 | phi = FindPhi(dx,dy) |
---|
1548 | |
---|
1549 | //get scattering angle to project onto flat detector => Qr = qval*cos(theta) |
---|
1550 | dist = sqrt(dx^2 + dy^2) |
---|
1551 | two_theta = atan(dist/sdd) |
---|
1552 | |
---|
1553 | qy = qval*cos(two_theta/2)*sin(phi) |
---|
1554 | |
---|
1555 | return qy |
---|
1556 | End |
---|
1557 | |
---|
1558 | //calculates just the z-component of the q-vector, not measured on the detector |
---|
1559 | //input/output is the same as CalcQval() |
---|
1560 | //ALL inputs are in detector coordinates |
---|
1561 | //NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector |
---|
1562 | //sdd is in meters |
---|
1563 | //wavelength is in Angstroms |
---|
1564 | // |
---|
1565 | // not actually used, but here for completeness if anyone asks |
---|
1566 | // |
---|
1567 | Function CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) |
---|
1568 | Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSize |
---|
1569 | |
---|
1570 | Variable dy,qval,dx,phi,qz,dist,two_theta |
---|
1571 | |
---|
1572 | qval = CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSize) |
---|
1573 | |
---|
1574 | sdd *=100 //convert to cm |
---|
1575 | |
---|
1576 | //get scattering angle to project onto flat detector => Qr = qval*cos(theta) |
---|
1577 | dx = (xaxval - xctr)*pixSize //delta x in cm |
---|
1578 | dy = (yaxval - yctr)*pixSize //delta y in cm |
---|
1579 | dist = sqrt(dx^2 + dy^2) |
---|
1580 | two_theta = atan(dist/sdd) |
---|
1581 | |
---|
1582 | qz = qval*sin(two_theta/2) |
---|
1583 | |
---|
1584 | return qz |
---|
1585 | End |
---|
1586 | |
---|
1587 | //for command-line testing, replace the function declaration |
---|
1588 | //Function FindQxQy(qq,phi) |
---|
1589 | // Variable qq,phi |
---|
1590 | // Variable qx,qy |
---|
1591 | // |
---|
1592 | // |
---|
1593 | ThreadSafe Function FindQxQy(qq,phi,qx,qy) |
---|
1594 | Variable qq,phi,&qx,&qy |
---|
1595 | |
---|
1596 | qx = sqrt(qq^2/(1+tan(phi)*tan(phi))) |
---|
1597 | qy = qx*tan(phi) |
---|
1598 | |
---|
1599 | if(phi >= 0 && phi <= pi/2) |
---|
1600 | qx = abs(qx) |
---|
1601 | qy = abs(qy) |
---|
1602 | endif |
---|
1603 | |
---|
1604 | if(phi > pi/2 && phi <= pi) |
---|
1605 | qx = -abs(qx) |
---|
1606 | qy = abs(qy) |
---|
1607 | endif |
---|
1608 | |
---|
1609 | if(phi > pi && phi <= pi*3/2) |
---|
1610 | qx = -abs(qx) |
---|
1611 | qy = -abs(qy) |
---|
1612 | endif |
---|
1613 | |
---|
1614 | if(phi > pi*3/2 && phi < 2*pi) |
---|
1615 | qx = abs(qx) |
---|
1616 | qy = -abs(qy) |
---|
1617 | endif |
---|
1618 | |
---|
1619 | |
---|
1620 | // Print "recalculated qx,qy,q = ",qx,qy,sqrt(qx*qx+qy*qy) |
---|
1621 | |
---|
1622 | return(0) |
---|
1623 | end |
---|
1624 | |
---|
1625 | |
---|
1626 | // 7 MAR 2011 SRK |
---|
1627 | // |
---|
1628 | // calculate the resolution smearing AAO |
---|
1629 | // |
---|
1630 | // - many of the form factor calculations are threaded, so they benefit |
---|
1631 | // from being passed large numbers of q-values at once, rather than suffering the |
---|
1632 | // overhead penalty of setting up threads. |
---|
1633 | // |
---|
1634 | // In general, single integral functions benefit from this, multiple integrals not so much. |
---|
1635 | // As an example, a fit using SmearedCylinderForm took 4.3s passing nord (=20) q-values |
---|
1636 | // at a time, but only 1.1s by passing all (Nq*nord) q-values at once. For Cyl_polyRad, |
---|
1637 | // the difference was not so large, 16.2s vs. 11.9s. This is due to CylPolyRad being a |
---|
1638 | // double integral and slow enough of a calculation that passing even 20 points at once |
---|
1639 | // provides some speedup. |
---|
1640 | // |
---|
1641 | //// APRIL 2011 *** this function is now cursor aware. The whole input x-wave is interpolated |
---|
1642 | // |
---|
1643 | // |
---|
1644 | // 4 MAR 2011 SRK |
---|
1645 | // Note: In John's paper, he integrated the Gaussian to +/- 3 sigma and then renormalized |
---|
1646 | // to an integral of 1. This "truncated" gaussian was a somewhat better approximation |
---|
1647 | // to the triangular resolution function. Here, I integrate to +/- 3 sigma and |
---|
1648 | // do not renormalize the integral to 1. Hence the smeared calculation is 0.27% low. |
---|
1649 | // This is easily seen by smearing a constant value. |
---|
1650 | // |
---|
1651 | // Using 5 quadrature points is not recommended, as it doesn't normalize properly using .9973 |
---|
1652 | // -- instead, it normalizes to 1.0084. |
---|
1653 | // |
---|
1654 | Function Smear_Model_N_AAO(fcn,w,x,resW,wi,zi,nord,sm_ans) |
---|
1655 | FUNCREF SANSModelAAO_proto fcn |
---|
1656 | Wave w //coefficients of function fcn(w,x) |
---|
1657 | WAVE x //x-value (q) for the calculation THIS IS PASSED IN AS A WAVE |
---|
1658 | Wave resW // Nx4 or NxN matrix of resolution |
---|
1659 | Wave wi //weight wave |
---|
1660 | Wave zi //abscissa wave |
---|
1661 | Variable nord //order of integration |
---|
1662 | Wave sm_ans // wave returned with the smeared model |
---|
1663 | |
---|
1664 | NVAR dQv = root:Packages:NIST:USANS_dQv |
---|
1665 | |
---|
1666 | // local variables |
---|
1667 | Variable ii,jj |
---|
1668 | Variable normalize=1 |
---|
1669 | Variable nTot,num,block_sum |
---|
1670 | |
---|
1671 | |
---|
1672 | // current x point is the q-value for evaluation |
---|
1673 | // |
---|
1674 | // * for the input x, the resolution function waves are interpolated to get the correct values for |
---|
1675 | // sigq, qbar and shad - since the model x-spacing may not be the same as |
---|
1676 | // the experimental QSIG data. This is always the case when curve fitting, since fit_wave is |
---|
1677 | // Igor-defined as 200 points and has its own (linear) q-(x)-scaling which will be quite different |
---|
1678 | // from experimental data. |
---|
1679 | // **note** if the (x) passed in is the experimental q-values, these values are |
---|
1680 | // returned from the interpolation (as expected) |
---|
1681 | |
---|
1682 | Make/O/D/N=(numpnts(x)) sigQ,qbar,shad,qvals,va,vb |
---|
1683 | Make/O/D/N=(DimSize(resW, 0)) tmpsigQ,tmpqbar,tmpshad,tmpqvals |
---|
1684 | tmpsigq = resW[p][0] //std dev of resolution fn |
---|
1685 | tmpqbar = resW[p][1] //mean q-value |
---|
1686 | tmpshad = resW[p][2] //beamstop shadow factor |
---|
1687 | tmpqvals = resW[p][3] //q-values where R(q) is known |
---|
1688 | |
---|
1689 | //interpolate the whole input x-wave to make sure that the resolution and input x are in sync if cursors are used |
---|
1690 | shad = interp(x,tmpqvals,tmpshad) |
---|
1691 | qbar = interp(x,tmpqvals,tmpqbar) |
---|
1692 | sigq = interp(x,tmpqvals,tmpsigq) |
---|
1693 | |
---|
1694 | // if USANS data, handle separately |
---|
1695 | // -- but this would only ever be used if the calculation was forced to use trapezoid integration |
---|
1696 | if ( ! isSANSResolution(sigq[0]) ) |
---|
1697 | //smear with the USANS routine |
---|
1698 | // Make global string and local variables |
---|
1699 | // now data folder aware, necessary for GlobalFit = FULL path to wave |
---|
1700 | String/G gTrap_coefStr = GetWavesDataFolder(w, 2 ) |
---|
1701 | Variable maxiter=20, tol=1e-4,uva,uvb |
---|
1702 | |
---|
1703 | num=numpnts(x) |
---|
1704 | // set up limits for the integration |
---|
1705 | uva=0 |
---|
1706 | uvb=abs(dQv) |
---|
1707 | |
---|
1708 | //loop over the q-values |
---|
1709 | for(jj=0;jj<num;jj+=1) |
---|
1710 | Variable/G gEvalQval = x[jj] |
---|
1711 | |
---|
1712 | // call qtrap to do actual work |
---|
1713 | sm_ans[jj] = qtrap_USANS(fcn,uva,uvb,tol,maxiter) |
---|
1714 | sm_ans[jj] /= (uvb - uva) |
---|
1715 | endfor |
---|
1716 | |
---|
1717 | return(0) |
---|
1718 | endif |
---|
1719 | |
---|
1720 | |
---|
1721 | // now the idea is to calculate a long vector of all of the zi's (Nq * nord) |
---|
1722 | // and pass these AAO to the AAO function, to make the most use of the threading |
---|
1723 | // passing repeated short lengths of q to the function can actually be slower |
---|
1724 | // due to the overhead. |
---|
1725 | |
---|
1726 | num = numpnts(x) |
---|
1727 | nTot = nord*num |
---|
1728 | |
---|
1729 | Make/O/D/N=(nTot) Resoln,yyy,xGauss,wts |
---|
1730 | |
---|
1731 | //loop over q |
---|
1732 | for(jj=0;jj<num;jj+=1) |
---|
1733 | |
---|
1734 | //for each q, set up the integration range |
---|
1735 | // end points of integration limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero |
---|
1736 | // +/- 3 sigq catches 99.73% of distrubution |
---|
1737 | // change limits (and spacing of zi) at each evaluation based on R() |
---|
1738 | //integration from va to vb |
---|
1739 | |
---|
1740 | va[jj] = -3*sigq[jj] + qbar[jj] |
---|
1741 | if (va[jj]<0) |
---|
1742 | va[jj]=0 //to avoid numerical error when va<0 (-ve q-value) |
---|
1743 | // Print "truncated Gaussian at nominal q = ",x |
---|
1744 | endif |
---|
1745 | vb[jj] = 3*sigq[jj] + qbar[jj] |
---|
1746 | |
---|
1747 | // loop over the Gauss points |
---|
1748 | for(ii=0;ii<nord;ii+=1) |
---|
1749 | // calculate Gauss points on integration interval (q-value for evaluation) |
---|
1750 | xGauss[nord*jj+ii] = ( zi[ii]*(vb[jj]-va[jj]) + vb[jj] + va[jj] )/2.0 |
---|
1751 | // calculate resolution function at input q-value (use the interpolated values and zi) |
---|
1752 | Resoln[nord*jj+ii] = shad[jj]/sqrt(2*pi*sigq[jj]*sigq[jj]) |
---|
1753 | Resoln[nord*jj+ii] *= exp((-1*(xGauss[nord*jj+ii] - qbar[jj])^2)/(2*sigq[jj]*sigq[jj])) |
---|
1754 | // Resoln[nord*jj+ii] *= exp((-1*(xGauss[nord*jj+ii] - qvals[jj])^2)/(2*sigq[jj]*sigq[jj])) //WRONG, but just for testing |
---|
1755 | // carry a copy of the weights |
---|
1756 | wts[nord*jj+ii] = wi[ii] |
---|
1757 | endfor // end of loop over quadrature points |
---|
1758 | |
---|
1759 | endfor //loop over q |
---|
1760 | |
---|
1761 | //calculate AAO |
---|
1762 | yyy = 0 |
---|
1763 | fcn(w,yyy,xGauss) //yyy is the return value as a wave |
---|
1764 | |
---|
1765 | //multiply by weights |
---|
1766 | yyy *= wts*Resoln //multiply function by resolution and weights |
---|
1767 | |
---|
1768 | //sum up blockwise to get the final answer |
---|
1769 | for(jj=0;jj<num;jj+=1) |
---|
1770 | block_sum = 0 |
---|
1771 | for(ii=0;ii<nord;ii+=1) |
---|
1772 | block_sum += yyy[nord*jj+ii] |
---|
1773 | endfor |
---|
1774 | sm_ans[jj] = (vb[jj]-va[jj])/2.0*block_sum |
---|
1775 | endfor |
---|
1776 | |
---|
1777 | |
---|
1778 | // then normalize for +/- 3 sigma ONLY |
---|
1779 | if(nord == 5) |
---|
1780 | normalize = 1.0057 //empirical correction, N=5 shouldn't be any different |
---|
1781 | else |
---|
1782 | normalize = 0.9973 |
---|
1783 | endif |
---|
1784 | |
---|
1785 | sm_ans /= normalize |
---|
1786 | |
---|
1787 | KillWaves/Z tmpsigQ,tmpqbar,tmpshad,tmpqvals |
---|
1788 | |
---|
1789 | return(0) |
---|
1790 | |
---|
1791 | End |
---|
1792 | |
---|