1 | #pragma rtGlobals=1 // Use modern global access method. |
---|
2 | #pragma version=3.00 |
---|
3 | #pragma IgorVersion=4.0 |
---|
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 | Function Make5GaussPoints(w5,z5) |
---|
32 | Wave w5,z5 |
---|
33 | |
---|
34 | // printf "in make Gauss Pts\r" |
---|
35 | |
---|
36 | z5[0] = -.906179845938664 |
---|
37 | z5[1] = -.538469310105683 |
---|
38 | z5[2] = -.0000000000000 |
---|
39 | z5[3] = .538469310105683 |
---|
40 | z5[4] = .906179845938664 |
---|
41 | |
---|
42 | w5[0] = .236926885056189 |
---|
43 | w5[1] = .478628670499366 |
---|
44 | w5[2] = .56888888888889 |
---|
45 | w5[3] = .478628670499366 |
---|
46 | w5[4] = .236926885056189 |
---|
47 | |
---|
48 | // printf "w[0],z[0] = %g %g\r", w5[0],z5[0] |
---|
49 | End |
---|
50 | |
---|
51 | Function Make10GaussPoints(w10,z10) |
---|
52 | Wave w10,z10 |
---|
53 | |
---|
54 | // printf "in make Gauss Pts\r" |
---|
55 | z10[0] = -.973906528517172 |
---|
56 | z10[1] = -.865063366688985 |
---|
57 | z10[2] = -.679409568299024 |
---|
58 | z10[3] = -.433395394129247 |
---|
59 | z10[4] = -.148874338981631 |
---|
60 | z10[5] = .148874338981631 |
---|
61 | z10[6] = .433395394129247 |
---|
62 | z10[7] = .679409568299024 |
---|
63 | z10[8] = .865063366688985 |
---|
64 | z10[9] = .973906528517172 |
---|
65 | |
---|
66 | w10[0] = .066671344308688 |
---|
67 | w10[1] = 0.149451349150581 |
---|
68 | w10[2] = 0.219086362515982 |
---|
69 | w10[3] = .269266719309996 |
---|
70 | w10[4] = 0.295524224714753 |
---|
71 | w10[5] = 0.295524224714753 |
---|
72 | w10[6] = .269266719309996 |
---|
73 | w10[7] = 0.219086362515982 |
---|
74 | w10[8] = 0.149451349150581 |
---|
75 | w10[9] = .066671344308688 |
---|
76 | |
---|
77 | // printf "w[0],z[0] = %g %g\r", w10[0],z10[0] |
---|
78 | End |
---|
79 | |
---|
80 | Function Make20GaussPoints(w20,z20) |
---|
81 | Wave w20,z20 |
---|
82 | |
---|
83 | // printf "in make Gauss Pts\r" |
---|
84 | |
---|
85 | z20[0] = -.993128599185095 |
---|
86 | z20[1] = -.963971927277914 |
---|
87 | z20[2] = -.912234428251326 |
---|
88 | z20[3] = -.839116971822219 |
---|
89 | z20[4] = -.746331906460151 |
---|
90 | z20[5] = -.636053680726515 |
---|
91 | z20[6] = -.510867001950827 |
---|
92 | z20[7] = -.37370608871542 |
---|
93 | z20[8] = -.227785851141645 |
---|
94 | z20[9] = -.076526521133497 |
---|
95 | z20[10] = .0765265211334973 |
---|
96 | z20[11] = .227785851141645 |
---|
97 | z20[12] = .37370608871542 |
---|
98 | z20[13] = .510867001950827 |
---|
99 | z20[14] = .636053680726515 |
---|
100 | z20[15] = .746331906460151 |
---|
101 | z20[16] = .839116971822219 |
---|
102 | z20[17] = .912234428251326 |
---|
103 | z20[18] = .963971927277914 |
---|
104 | z20[19] = .993128599185095 |
---|
105 | |
---|
106 | w20[0] = .0176140071391521 |
---|
107 | w20[1] = .0406014298003869 |
---|
108 | w20[2] = .0626720483341091 |
---|
109 | w20[3] = .0832767415767047 |
---|
110 | w20[4] = .10193011981724 |
---|
111 | w20[5] = .118194531961518 |
---|
112 | w20[6] = .131688638449177 |
---|
113 | w20[7] = .142096109318382 |
---|
114 | w20[8] = .149172986472604 |
---|
115 | w20[9] = .152753387130726 |
---|
116 | w20[10] = .152753387130726 |
---|
117 | w20[11] = .149172986472604 |
---|
118 | w20[12] = .142096109318382 |
---|
119 | w20[13] = .131688638449177 |
---|
120 | w20[14] = .118194531961518 |
---|
121 | w20[15] = .10193011981724 |
---|
122 | w20[16] = .0832767415767047 |
---|
123 | w20[17] = .0626720483341091 |
---|
124 | w20[18] = .0406014298003869 |
---|
125 | w20[19] = .0176140071391521 |
---|
126 | // printf "w[0],z[0] = %g %g\r", w20[0],z20[0] |
---|
127 | End |
---|
128 | |
---|
129 | |
---|
130 | |
---|
131 | |
---|
132 | Function Make76GaussPoints(w76,z76) |
---|
133 | Wave w76,z76 |
---|
134 | |
---|
135 | // printf "in make Gauss Pts\r" |
---|
136 | |
---|
137 | z76[0] = .999505948362153*(-1.0) |
---|
138 | z76[75] = -.999505948362153*(-1.0) |
---|
139 | z76[1] = .997397786355355*(-1.0) |
---|
140 | z76[74] = -.997397786355355*(-1.0) |
---|
141 | z76[2] = .993608772723527*(-1.0) |
---|
142 | z76[73] = -.993608772723527*(-1.0) |
---|
143 | z76[3] = .988144453359837*(-1.0) |
---|
144 | z76[72] = -.988144453359837*(-1.0) |
---|
145 | z76[4] = .981013938975656*(-1.0) |
---|
146 | z76[71] = -.981013938975656*(-1.0) |
---|
147 | z76[5] = .972229228520377*(-1.0) |
---|
148 | z76[70] = -.972229228520377*(-1.0) |
---|
149 | z76[6] = .961805126758768*(-1.0) |
---|
150 | z76[69] = -.961805126758768*(-1.0) |
---|
151 | z76[7] = .949759207710896*(-1.0) |
---|
152 | z76[68] = -.949759207710896*(-1.0) |
---|
153 | z76[8] = .936111781934811*(-1.0) |
---|
154 | z76[67] = -.936111781934811*(-1.0) |
---|
155 | z76[9] = .92088586125215*(-1.0) |
---|
156 | z76[66] = -.92088586125215*(-1.0) |
---|
157 | z76[10] = .904107119545567*(-1.0) |
---|
158 | z76[65] = -.904107119545567*(-1.0) |
---|
159 | z76[11] = .885803849292083*(-1.0) |
---|
160 | z76[64] = -.885803849292083*(-1.0) |
---|
161 | z76[12] = .866006913771982*(-1.0) |
---|
162 | z76[63] = -.866006913771982*(-1.0) |
---|
163 | z76[13] = .844749694983342*(-1.0) |
---|
164 | z76[62] = -.844749694983342*(-1.0) |
---|
165 | z76[14] = .822068037328975*(-1.0) |
---|
166 | z76[61] = -.822068037328975*(-1.0) |
---|
167 | z76[15] = .7980001871612*(-1.0) |
---|
168 | z76[60] = -.7980001871612*(-1.0) |
---|
169 | z76[16] = .77258672828181*(-1.0) |
---|
170 | z76[59] = -.77258672828181*(-1.0) |
---|
171 | z76[17] = .74587051350361*(-1.0) |
---|
172 | z76[58] = -.74587051350361*(-1.0) |
---|
173 | z76[18] = .717896592387704*(-1.0) |
---|
174 | z76[57] = -.717896592387704*(-1.0) |
---|
175 | z76[19] = .688712135277641*(-1.0) |
---|
176 | z76[56] = -.688712135277641*(-1.0) |
---|
177 | z76[20] = .658366353758143*(-1.0) |
---|
178 | z76[55] = -.658366353758143*(-1.0) |
---|
179 | z76[21] = .626910417672267*(-1.0) |
---|
180 | z76[54] = -.626910417672267*(-1.0) |
---|
181 | z76[22] = .594397368836793*(-1.0) |
---|
182 | z76[53] = -.594397368836793*(-1.0) |
---|
183 | z76[23] = .560882031601237*(-1.0) |
---|
184 | z76[52] = -.560882031601237*(-1.0) |
---|
185 | z76[24] = .526420920401243*(-1.0) |
---|
186 | z76[51] = -.526420920401243*(-1.0) |
---|
187 | z76[25] = .491072144462194*(-1.0) |
---|
188 | z76[50] = -.491072144462194*(-1.0) |
---|
189 | z76[26] = .454895309813726*(-1.0) |
---|
190 | z76[49] = -.454895309813726*(-1.0) |
---|
191 | z76[27] = .417951418780327*(-1.0) |
---|
192 | z76[48] = -.417951418780327*(-1.0) |
---|
193 | z76[28] = .380302767117504*(-1.0) |
---|
194 | z76[47] = -.380302767117504*(-1.0) |
---|
195 | z76[29] = .342012838966962*(-1.0) |
---|
196 | z76[46] = -.342012838966962*(-1.0) |
---|
197 | z76[30] = .303146199807908*(-1.0) |
---|
198 | z76[45] = -.303146199807908*(-1.0) |
---|
199 | z76[31] = .263768387584994*(-1.0) |
---|
200 | z76[44] = -.263768387584994*(-1.0) |
---|
201 | z76[32] = .223945802196474*(-1.0) |
---|
202 | z76[43] = -.223945802196474*(-1.0) |
---|
203 | z76[33] = .183745593528914*(-1.0) |
---|
204 | z76[42] = -.183745593528914*(-1.0) |
---|
205 | z76[34] = .143235548227268*(-1.0) |
---|
206 | z76[41] = -.143235548227268*(-1.0) |
---|
207 | z76[35] = .102483975391227*(-1.0) |
---|
208 | z76[40] = -.102483975391227*(-1.0) |
---|
209 | z76[36] = .0615595913906112*(-1.0) |
---|
210 | z76[39] = -.0615595913906112*(-1.0) |
---|
211 | z76[37] = .0205314039939986*(-1.0) |
---|
212 | z76[38] = -.0205314039939986*(-1.0) |
---|
213 | |
---|
214 | w76[0] = .00126779163408536 |
---|
215 | w76[75] = .00126779163408536 |
---|
216 | w76[1] = .00294910295364247 |
---|
217 | w76[74] = .00294910295364247 |
---|
218 | w76[2] = .00462793522803742 |
---|
219 | w76[73] = .00462793522803742 |
---|
220 | w76[3] = .00629918049732845 |
---|
221 | w76[72] = .00629918049732845 |
---|
222 | w76[4] = .00795984747723973 |
---|
223 | w76[71] = .00795984747723973 |
---|
224 | w76[5] = .00960710541471375 |
---|
225 | w76[70] = .00960710541471375 |
---|
226 | w76[6] = .0112381685696677 |
---|
227 | w76[69] = .0112381685696677 |
---|
228 | w76[7] = .0128502838475101 |
---|
229 | w76[68] = .0128502838475101 |
---|
230 | w76[8] = .0144407317482767 |
---|
231 | w76[67] = .0144407317482767 |
---|
232 | w76[9] = .0160068299122486 |
---|
233 | w76[66] = .0160068299122486 |
---|
234 | w76[10] = .0175459372914742 |
---|
235 | w76[65] = .0175459372914742 |
---|
236 | w76[11] = .0190554584671906 |
---|
237 | w76[64] = .0190554584671906 |
---|
238 | w76[12] = .020532847967908 |
---|
239 | w76[63] = .020532847967908 |
---|
240 | w76[13] = .0219756145344162 |
---|
241 | w76[62] = .0219756145344162 |
---|
242 | w76[14] = .0233813253070112 |
---|
243 | w76[61] = .0233813253070112 |
---|
244 | w76[15] = .0247476099206597 |
---|
245 | w76[60] = .0247476099206597 |
---|
246 | w76[16] = .026072164497986 |
---|
247 | w76[59] = .026072164497986 |
---|
248 | w76[17] = .0273527555318275 |
---|
249 | w76[58] = .0273527555318275 |
---|
250 | w76[18] = .028587223650054 |
---|
251 | w76[57] = .028587223650054 |
---|
252 | w76[19] = .029773487255905 |
---|
253 | w76[56] = .029773487255905 |
---|
254 | w76[20] = .0309095460374916 |
---|
255 | w76[55] = .0309095460374916 |
---|
256 | w76[21] = .0319934843404216 |
---|
257 | w76[54] = .0319934843404216 |
---|
258 | w76[22] = .0330234743977917 |
---|
259 | w76[53] = .0330234743977917 |
---|
260 | w76[23] = .0339977794120564 |
---|
261 | w76[52] = .0339977794120564 |
---|
262 | w76[24] = .0349147564835508 |
---|
263 | w76[51] = .0349147564835508 |
---|
264 | w76[25] = .0357728593807139 |
---|
265 | w76[50] = .0357728593807139 |
---|
266 | w76[26] = .0365706411473296 |
---|
267 | w76[49] = .0365706411473296 |
---|
268 | w76[27] = .0373067565423816 |
---|
269 | w76[48] = .0373067565423816 |
---|
270 | w76[28] = .0379799643084053 |
---|
271 | w76[47] = .0379799643084053 |
---|
272 | w76[29] = .0385891292645067 |
---|
273 | w76[46] = .0385891292645067 |
---|
274 | w76[30] = .0391332242205184 |
---|
275 | w76[45] = .0391332242205184 |
---|
276 | w76[31] = .0396113317090621 |
---|
277 | w76[44] = .0396113317090621 |
---|
278 | w76[32] = .0400226455325968 |
---|
279 | w76[43] = .0400226455325968 |
---|
280 | w76[33] = .040366472122844 |
---|
281 | w76[42] = .040366472122844 |
---|
282 | w76[34] = .0406422317102947 |
---|
283 | w76[41] = .0406422317102947 |
---|
284 | w76[35] = .0408494593018285 |
---|
285 | w76[40] = .0408494593018285 |
---|
286 | w76[36] = .040987805464794 |
---|
287 | w76[39] = .040987805464794 |
---|
288 | w76[37] = .0410570369162294 |
---|
289 | w76[38] = .0410570369162294 |
---|
290 | |
---|
291 | End //Make76GaussPoints() |
---|
292 | |
---|
293 | //////////// |
---|
294 | Function IntegrateFn5(fcn,loLim,upLim,w,x) |
---|
295 | FUNCREF GenericQuadrature_proto fcn |
---|
296 | Variable loLim,upLim //limits of integration |
---|
297 | Wave w //coefficients of function fcn(w,x) |
---|
298 | Variable x //x-value (q) for the calculation |
---|
299 | |
---|
300 | // local variables |
---|
301 | Variable nord,ii,va,vb,summ,yyy,zi |
---|
302 | Variable answer,dum |
---|
303 | String weightStr,zStr |
---|
304 | |
---|
305 | weightStr = "gauss5wt" |
---|
306 | zStr = "gauss5z" |
---|
307 | nord = 5 |
---|
308 | |
---|
309 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
310 | Make/D/N=5 $weightStr,$zStr |
---|
311 | Wave w5 = $weightStr |
---|
312 | Wave z5 = $zStr // wave references to pass |
---|
313 | Make5GaussPoints(w5,z5) |
---|
314 | else |
---|
315 | if(exists(weightStr) > 1) |
---|
316 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
317 | endif |
---|
318 | Wave w5 = $weightStr |
---|
319 | Wave z5 = $zStr // create the wave references |
---|
320 | endif |
---|
321 | |
---|
322 | //limits of integration are input to function |
---|
323 | va = loLim |
---|
324 | vb = upLim |
---|
325 | // Using 5 Gauss points |
---|
326 | // remember to index from 0,size-1 |
---|
327 | |
---|
328 | summ = 0.0 // initialize integral |
---|
329 | ii=0 // loop counter |
---|
330 | do |
---|
331 | // calculate Gauss points on integration interval (q-value for evaluation) |
---|
332 | zi = ( z5[ii]*(vb-va) + vb + va )/2.0 |
---|
333 | //calculate partial sum for the passed-in model function |
---|
334 | yyy = w5[ii] * fcn(w,x,zi) |
---|
335 | summ += yyy //add to the running total of the quadrature |
---|
336 | ii+=1 |
---|
337 | while (ii<nord) // end of loop over quadrature points |
---|
338 | |
---|
339 | // calculate value of integral to return |
---|
340 | answer = (vb-va)/2.0*summ |
---|
341 | |
---|
342 | Return (answer) |
---|
343 | End |
---|
344 | /////////////////////////////////////////////////////////////// |
---|
345 | |
---|
346 | //////////// |
---|
347 | Function IntegrateFn10(fcn,loLim,upLim,w,x) |
---|
348 | FUNCREF GenericQuadrature_proto fcn |
---|
349 | Variable loLim,upLim //limits of integration |
---|
350 | Wave w //coefficients of function fcn(w,x) |
---|
351 | Variable x //x-value (q) for the calculation |
---|
352 | |
---|
353 | // local variables |
---|
354 | Variable nord,ii,va,vb,summ,yyy,zi |
---|
355 | Variable answer,dum |
---|
356 | String weightStr,zStr |
---|
357 | |
---|
358 | weightStr = "gauss10wt" |
---|
359 | zStr = "gauss10z" |
---|
360 | nord = 10 |
---|
361 | |
---|
362 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
363 | Make/D/N=10 $weightStr,$zStr |
---|
364 | Wave w10 = $weightStr |
---|
365 | Wave z10 = $zStr // wave references to pass |
---|
366 | Make10GaussPoints(w10,z10) |
---|
367 | else |
---|
368 | if(exists(weightStr) > 1) |
---|
369 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
370 | endif |
---|
371 | Wave w10 = $weightStr |
---|
372 | Wave z10 = $zStr // create the wave references |
---|
373 | endif |
---|
374 | |
---|
375 | //limits of integration are input to function |
---|
376 | va = loLim |
---|
377 | vb = upLim |
---|
378 | // Using 10 Gauss points |
---|
379 | // remember to index from 0,size-1 |
---|
380 | |
---|
381 | summ = 0.0 // initialize integral |
---|
382 | ii=0 // loop counter |
---|
383 | do |
---|
384 | // calculate Gauss points on integration interval (q-value for evaluation) |
---|
385 | zi = ( z10[ii]*(vb-va) + vb + va )/2.0 |
---|
386 | //calculate partial sum for the passed-in model function |
---|
387 | yyy = w10[ii] * fcn(w,x,zi) |
---|
388 | summ += yyy //add to the running total of the quadrature |
---|
389 | ii+=1 |
---|
390 | while (ii<nord) // end of loop over quadrature points |
---|
391 | |
---|
392 | // calculate value of integral to return |
---|
393 | answer = (vb-va)/2.0*summ |
---|
394 | |
---|
395 | Return (answer) |
---|
396 | End |
---|
397 | /////////////////////////////////////////////////////////////// |
---|
398 | |
---|
399 | //////////// |
---|
400 | Function IntegrateFn20(fcn,loLim,upLim,w,x) |
---|
401 | FUNCREF GenericQuadrature_proto fcn |
---|
402 | Variable loLim,upLim //limits of integration |
---|
403 | Wave w //coefficients of function fcn(w,x) |
---|
404 | Variable x //x-value (q) for the calculation |
---|
405 | |
---|
406 | // local variables |
---|
407 | Variable nord,ii,va,vb,summ,yyy,zi |
---|
408 | Variable answer,dum |
---|
409 | String weightStr,zStr |
---|
410 | |
---|
411 | weightStr = "gauss20wt" |
---|
412 | zStr = "gauss20z" |
---|
413 | nord = 20 |
---|
414 | |
---|
415 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
416 | Make/D/N=20 $weightStr,$zStr |
---|
417 | Wave w20 = $weightStr |
---|
418 | Wave z20 = $zStr // wave references to pass |
---|
419 | Make20GaussPoints(w20,z20) |
---|
420 | else |
---|
421 | if(exists(weightStr) > 1) |
---|
422 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
423 | endif |
---|
424 | Wave w20 = $weightStr |
---|
425 | Wave z20 = $zStr // create the wave references |
---|
426 | endif |
---|
427 | |
---|
428 | //limits of integration are input to function |
---|
429 | va = loLim |
---|
430 | vb = upLim |
---|
431 | // Using 20 Gauss points |
---|
432 | // remember to index from 0,size-1 |
---|
433 | |
---|
434 | summ = 0.0 // initialize integral |
---|
435 | ii=0 // loop counter |
---|
436 | do |
---|
437 | // calculate Gauss points on integration interval (q-value for evaluation) |
---|
438 | zi = ( z20[ii]*(vb-va) + vb + va )/2.0 |
---|
439 | //calculate partial sum for the passed-in model function |
---|
440 | yyy = w20[ii] * fcn(w,x,zi) |
---|
441 | summ += yyy //add to the running total of the quadrature |
---|
442 | ii+=1 |
---|
443 | while (ii<nord) // end of loop over quadrature points |
---|
444 | |
---|
445 | // calculate value of integral to return |
---|
446 | answer = (vb-va)/2.0*summ |
---|
447 | |
---|
448 | Return (answer) |
---|
449 | End |
---|
450 | /////////////////////////////////////////////////////////////// |
---|
451 | |
---|
452 | Function IntegrateFn76(fcn,loLim,upLim,w,x) |
---|
453 | FUNCREF GenericQuadrature_proto fcn |
---|
454 | Variable loLim,upLim //limits of integration |
---|
455 | Wave w //coefficients of function fcn(w,x) |
---|
456 | Variable x //x-value (q) for the calculation |
---|
457 | |
---|
458 | //**** The coefficient wave is passed into this function and straight through to the unsmeared model function |
---|
459 | |
---|
460 | // local variables |
---|
461 | Variable nord,ii,va,vb,summ,yyy,zi |
---|
462 | Variable answer,dum |
---|
463 | String weightStr,zStr |
---|
464 | |
---|
465 | weightStr = "gauss76wt" |
---|
466 | zStr = "gauss76z" |
---|
467 | nord = 76 |
---|
468 | |
---|
469 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
470 | Make/D/N=76 $weightStr,$zStr |
---|
471 | Wave w76 = $weightStr |
---|
472 | Wave z76 = $zStr // wave references to pass |
---|
473 | Make76GaussPoints(w76,z76) |
---|
474 | else |
---|
475 | if(exists(weightStr) > 1) |
---|
476 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
477 | endif |
---|
478 | Wave w76 = $weightStr |
---|
479 | Wave z76 = $zStr // create the wave references |
---|
480 | endif |
---|
481 | |
---|
482 | //limits of integration are input to function |
---|
483 | va = loLim |
---|
484 | vb = upLim |
---|
485 | // Using 76 Gauss points |
---|
486 | // remember to index from 0,size-1 |
---|
487 | |
---|
488 | summ = 0.0 // initialize integral |
---|
489 | ii=0 // loop counter |
---|
490 | do |
---|
491 | // calculate Gauss points on integration interval (q-value for evaluation) |
---|
492 | zi = ( z76[ii]*(vb-va) + vb + va )/2.0 |
---|
493 | //calculate partial sum for the passed-in model function |
---|
494 | yyy = w76[ii] * fcn(w,x,zi) |
---|
495 | summ += yyy //add to the running total of the quadrature |
---|
496 | ii+=1 |
---|
497 | while (ii<nord) // end of loop over quadrature points |
---|
498 | |
---|
499 | // calculate value of integral to return |
---|
500 | answer = (vb-va)/2.0*summ |
---|
501 | |
---|
502 | Return (answer) |
---|
503 | End |
---|
504 | /////////////////////////////////////////////////////////////// |
---|
505 | |
---|
506 | //////Resolution Smearing Utilities |
---|
507 | |
---|
508 | // To check for the existence of all waves needed for smearing |
---|
509 | // returns 1 if any waves are missing, 0 if all is OK |
---|
510 | Function ResolutionWavesMissing() |
---|
511 | |
---|
512 | SVAR/Z sq = gSig_Q |
---|
513 | SVAR/Z qb = gQ_bar |
---|
514 | SVAR/Z sh = gShadow |
---|
515 | SVAR/Z gQ = gQVals |
---|
516 | |
---|
517 | Variable err=0 |
---|
518 | if(!SVAR_Exists(sq) || !SVAR_Exists(qb) || !SVAR_Exists(sh) || !SVAR_Exists(gQ)) |
---|
519 | DoAlert 0,"Some 6-column QSIG waves are missing. Re-load experimental data with LoadOneDData macro" |
---|
520 | return(1) |
---|
521 | endif |
---|
522 | |
---|
523 | if(WaveExists($sq) == 0) //wave ref does not exist |
---|
524 | err=1 |
---|
525 | endif |
---|
526 | if(WaveExists($qb) == 0) //wave ref does not exist |
---|
527 | err=1 |
---|
528 | endif |
---|
529 | if(WaveExists($sh) == 0) //wave ref does not exist |
---|
530 | err=1 |
---|
531 | endif |
---|
532 | if(WaveExists($gQ) == 0) //wave ref does not exist |
---|
533 | err=1 |
---|
534 | endif |
---|
535 | |
---|
536 | if(err) |
---|
537 | DoAlert 0,"Some 6-column QSIG waves are missing. Re-load experimental data with LoadOneDData macro" |
---|
538 | endif |
---|
539 | |
---|
540 | return(err) |
---|
541 | end |
---|
542 | |
---|
543 | |
---|
544 | //Utility function to smear model function with resolution |
---|
545 | // |
---|
546 | // call as in this example... |
---|
547 | // |
---|
548 | ////Function SmearedSphere_HS(w,x) :FitFunc |
---|
549 | //// Wave w |
---|
550 | //// Variable x |
---|
551 | //// |
---|
552 | ////// Variable timer=StartMSTimer |
---|
553 | //// Variable ans |
---|
554 | //// SVAR sq = gSig_Q |
---|
555 | //// SVAR qb = gQ_bar |
---|
556 | //// SVAR sh = gShadow |
---|
557 | //// SVAR gQ = gQVals |
---|
558 | //// |
---|
559 | //// ans = Smear_Model_20(Sphere_HS,$sq,$qb,$sh,$gQ,w,x) |
---|
560 | //// |
---|
561 | ////// Print "HS elapsed time(s) = ",StopMSTimer(timer)*1e-6 |
---|
562 | //// return(ans) |
---|
563 | ////End |
---|
564 | // |
---|
565 | // |
---|
566 | |
---|
567 | Function Smear_Model_76(fcn,sigq,qbar,shad,qvals,w,x) |
---|
568 | FUNCREF SANSModel_proto fcn |
---|
569 | Wave sigq //std dev of resolution fn |
---|
570 | Wave qbar //mean q-value |
---|
571 | Wave shad //beamstop shadow factor |
---|
572 | Wave qvals //q-values where R(q) is known |
---|
573 | Wave w //coefficients of function fcn(w,x) |
---|
574 | Variable x //x-value (q) for the calculation |
---|
575 | |
---|
576 | //**** The coefficient wave is passed into this function and straight through to the unsmeared model function |
---|
577 | |
---|
578 | // local variables |
---|
579 | Variable nord,ii,va,vb,contr,nden,summ,yyy,zi,q |
---|
580 | Variable answer,Resoln,i_shad,i_qbar,i_sigq |
---|
581 | String weightStr,zStr,weightStr1,zStr1 |
---|
582 | |
---|
583 | // weightStr = "gauss20wt" |
---|
584 | // zStr = "gauss20z" |
---|
585 | weightStr = "gauss76wt" |
---|
586 | zStr = "gauss76z" |
---|
587 | |
---|
588 | // if wt,z waves don't exist, create them (only check for weight, should really check for both) |
---|
589 | // 20 Gauss points (enough for smearing with Gaussian resolution function) |
---|
590 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
591 | Make/D/N=76 $weightStr,$zStr |
---|
592 | Wave w76 = $weightStr |
---|
593 | Wave z76 = $zStr // wave references to pass |
---|
594 | Make76GaussPoints(w76,z76) |
---|
595 | else |
---|
596 | if(exists(weightStr) > 1) |
---|
597 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
598 | endif |
---|
599 | Wave w76 = $weightStr |
---|
600 | Wave z76 = $zStr // create the wave references |
---|
601 | endif |
---|
602 | |
---|
603 | // current x point is the q-value for evaluation |
---|
604 | // |
---|
605 | // * for the input x, the resolution function waves are interpolated to get the correct values for |
---|
606 | // sigq, qbar and shad - since the model x-spacing may not be the same as |
---|
607 | // the experimental QSIG data. This is always the case when curve fitting, since fit_wave is |
---|
608 | // Igor-defined as 200 points and has its own (linear) q-(x)-scaling which will be quite different |
---|
609 | // from experimental data. |
---|
610 | // **note** if the (x) passed in is the experimental q-values, these values are |
---|
611 | // returned from the interpolation (as expected) |
---|
612 | |
---|
613 | i_shad = interp(x,qvals,shad) |
---|
614 | i_qbar = interp(x,qvals,qbar) |
---|
615 | i_sigq = interp(x,qvals,sigq) |
---|
616 | |
---|
617 | // set up the integration |
---|
618 | // number of Gauss Quadrature points |
---|
619 | |
---|
620 | if (i_sigq >= 0) |
---|
621 | |
---|
622 | // end points of integration |
---|
623 | // limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero |
---|
624 | // +/- 3 sigq catches 99.73% of distrubution |
---|
625 | // change limits (and spacing of zi) at each evaluation based on R() |
---|
626 | //integration from va to vb |
---|
627 | |
---|
628 | nord = 76 |
---|
629 | va = -3*i_sigq + i_qbar |
---|
630 | if (va<0) |
---|
631 | va=0 //to avoid numerical error when va<0 (-ve q-value) |
---|
632 | // Print "truncated Gaussian at nominal q = ",x |
---|
633 | endif |
---|
634 | vb = 3*i_sigq + i_qbar |
---|
635 | |
---|
636 | // Using 20 Gauss points |
---|
637 | // remember to index from 0,size-1 |
---|
638 | |
---|
639 | summ = 0.0 // initialize integral |
---|
640 | ii=0 // loop counter |
---|
641 | do |
---|
642 | // calculate Gauss points on integration interval (q-value for evaluation) |
---|
643 | zi = ( z76[ii]*(vb-va) + vb + va )/2.0 |
---|
644 | // calculate resolution function at input q-value (use the interpolated values and zi) |
---|
645 | Resoln = i_shad/sqrt(2*pi*i_sigq*i_sigq) |
---|
646 | Resoln *= exp((-1*(zi - i_qbar)^2)/(2*i_sigq*i_sigq)) |
---|
647 | //calculate partial sum for the passed-in model function |
---|
648 | yyy = w76[ii] * Resoln * fcn(w,zi) |
---|
649 | summ += yyy //add to the running total of the quadrature |
---|
650 | ii+=1 |
---|
651 | while (ii<nord) // end of loop over quadrature points |
---|
652 | |
---|
653 | // calculate value of integral to return |
---|
654 | answer = (vb-va)/2.0*summ |
---|
655 | // all scaling, background addition... etc. is done in the model calculation |
---|
656 | |
---|
657 | else |
---|
658 | //smear with the USANS routine |
---|
659 | // Make global string and local variables |
---|
660 | // now data folder aware, necessary for GlobalFit = FULL path to wave |
---|
661 | String/G gTrap_coefStr = GetWavesDataFolder(w, 2 ) |
---|
662 | Variable maxiter=20, tol=1e-4 |
---|
663 | |
---|
664 | // set up limits for the integration |
---|
665 | va=0 |
---|
666 | vb=abs(i_sigq) |
---|
667 | |
---|
668 | Variable/G gEvalQval = x |
---|
669 | |
---|
670 | // call qtrap to do actual work |
---|
671 | answer = qtrap_USANS(fcn,va,vb,tol,maxiter) |
---|
672 | answer /= vb |
---|
673 | |
---|
674 | endif |
---|
675 | |
---|
676 | Return (answer) |
---|
677 | |
---|
678 | End |
---|
679 | /////////////////////////////////////////////////////////////// |
---|
680 | |
---|
681 | Function Smear_Model_20(fcn,sigq,qbar,shad,qvals,w,x) |
---|
682 | FUNCREF SANSModel_proto fcn |
---|
683 | Wave sigq //std dev of resolution fn |
---|
684 | Wave qbar //mean q-value |
---|
685 | Wave shad //beamstop shadow factor |
---|
686 | Wave qvals //q-values where R(q) is known |
---|
687 | Wave w //coefficients of function fcn(w,x) |
---|
688 | Variable x //x-value (q) for the calculation |
---|
689 | |
---|
690 | //**** The coefficient wave is passed into this function and straight through to the unsmeared model function |
---|
691 | |
---|
692 | // local variables |
---|
693 | Variable nord,ii,va,vb,contr,nden,summ,yyy,zi,q |
---|
694 | Variable answer,Resoln,i_shad,i_qbar,i_sigq |
---|
695 | String weightStr,zStr,weightStr1,zStr1 |
---|
696 | |
---|
697 | weightStr = "gauss20wt" |
---|
698 | zStr = "gauss20z" |
---|
699 | weightStr1 = "gauss76wt" |
---|
700 | zStr1 = "gauss76z" |
---|
701 | |
---|
702 | // if wt,z waves don't exist, create them (only check for weight, should really check for both) |
---|
703 | // 20 Gauss points (enough for smearing with Gaussian resolution function) |
---|
704 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
705 | Make/D/N=20 $weightStr,$zStr |
---|
706 | Wave w20 = $weightStr |
---|
707 | Wave z20 = $zStr // wave references to pass |
---|
708 | Make20GaussPoints(w20,z20) |
---|
709 | else |
---|
710 | if(exists(weightStr) > 1) |
---|
711 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
712 | endif |
---|
713 | Wave w20 = $weightStr |
---|
714 | Wave z20 = $zStr // create the wave references |
---|
715 | endif |
---|
716 | |
---|
717 | // current x point is the q-value for evaluation |
---|
718 | // |
---|
719 | // * for the input x, the resolution function waves are interpolated to get the correct values for |
---|
720 | // sigq, qbar and shad - since the model x-spacing may not be the same as |
---|
721 | // the experimental QSIG data. This is always the case when curve fitting, since fit_wave is |
---|
722 | // Igor-defined as 200 points and has its own (linear) q-(x)-scaling which will be quite different |
---|
723 | // from experimental data. |
---|
724 | // **note** if the (x) passed in is the experimental q-values, these values are |
---|
725 | // returned from the interpolation (as expected) |
---|
726 | |
---|
727 | i_shad = interp(x,qvals,shad) |
---|
728 | i_qbar = interp(x,qvals,qbar) |
---|
729 | i_sigq = interp(x,qvals,sigq) |
---|
730 | |
---|
731 | // set up the integration |
---|
732 | // number of Gauss Quadrature points |
---|
733 | |
---|
734 | if (i_sigq >= 0) |
---|
735 | |
---|
736 | // end points of integration |
---|
737 | // limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero |
---|
738 | // +/- 3 sigq catches 99.73% of distrubution |
---|
739 | // change limits (and spacing of zi) at each evaluation based on R() |
---|
740 | //integration from va to vb |
---|
741 | |
---|
742 | nord = 20 |
---|
743 | va = -3*i_sigq + i_qbar |
---|
744 | if (va<0) |
---|
745 | va=0 //to avoid numerical error when va<0 (-ve q-value) |
---|
746 | // Print "truncated Gaussian at nominal q = ",x |
---|
747 | endif |
---|
748 | vb = 3*i_sigq + i_qbar |
---|
749 | |
---|
750 | // Using 20 Gauss points |
---|
751 | // remember to index from 0,size-1 |
---|
752 | |
---|
753 | summ = 0.0 // initialize integral |
---|
754 | ii=0 // loop counter |
---|
755 | do |
---|
756 | // calculate Gauss points on integration interval (q-value for evaluation) |
---|
757 | zi = ( z20[ii]*(vb-va) + vb + va )/2.0 |
---|
758 | // calculate resolution function at input q-value (use the interpolated values and zi) |
---|
759 | Resoln = i_shad/sqrt(2*pi*i_sigq*i_sigq) |
---|
760 | Resoln *= exp((-1*(zi - i_qbar)^2)/(2*i_sigq*i_sigq)) |
---|
761 | //calculate partial sum for the passed-in model function |
---|
762 | yyy = w20[ii] * Resoln * fcn(w,zi) |
---|
763 | summ += yyy //add to the running total of the quadrature |
---|
764 | ii+=1 |
---|
765 | while (ii<nord) // end of loop over quadrature points |
---|
766 | |
---|
767 | // calculate value of integral to return |
---|
768 | answer = (vb-va)/2.0*summ |
---|
769 | // all scaling, background addition... etc. is done in the model calculation |
---|
770 | |
---|
771 | else |
---|
772 | //smear with the USANS routine |
---|
773 | // Make global string and local variables |
---|
774 | // now data folder aware, necessary for GlobalFit = FULL path to wave |
---|
775 | String/G gTrap_coefStr = GetWavesDataFolder(w, 2 ) |
---|
776 | Variable maxiter=20, tol=1e-4 |
---|
777 | |
---|
778 | // set up limits for the integration |
---|
779 | va=0 |
---|
780 | vb=abs(i_sigq) |
---|
781 | |
---|
782 | Variable/G gEvalQval = x |
---|
783 | |
---|
784 | // call qtrap to do actual work |
---|
785 | answer = qtrap_USANS(fcn,va,vb,tol,maxiter) |
---|
786 | answer /= vb |
---|
787 | |
---|
788 | endif |
---|
789 | |
---|
790 | Return (answer) |
---|
791 | |
---|
792 | End |
---|
793 | /////////////////////////////////////////////////////////////// |
---|
794 | |
---|
795 | |
---|
796 | //resolution smearing, using only 5 Gauss points |
---|
797 | // |
---|
798 | Function Smear_Model_5(fcn,sigq,qbar,shad,qvals,w,x) |
---|
799 | FUNCREF SANSModel_proto fcn |
---|
800 | Wave sigq //std dev of resolution fn |
---|
801 | Wave qbar //mean q-value |
---|
802 | Wave shad //beamstop shadow factor |
---|
803 | Wave qvals //q-values where R(q) is known |
---|
804 | Wave w //coefficients of function fcn(w,x) |
---|
805 | Variable x //x-value (q) for the calculation |
---|
806 | |
---|
807 | //**** The coefficient wave is passed into this function and straight through to the unsmeared model function |
---|
808 | |
---|
809 | // local variables |
---|
810 | Variable nord,ii,va,vb,contr,nden,summ,yyy,zi,q |
---|
811 | Variable answer,Resoln,i_shad,i_qbar,i_sigq |
---|
812 | String weightStr,zStr |
---|
813 | |
---|
814 | weightStr = "gauss5wt" |
---|
815 | zStr = "gauss5z" |
---|
816 | |
---|
817 | // if wt,z waves don't exist, create them (only check for weight, should really check for both) |
---|
818 | // 5 Gauss points |
---|
819 | if (WaveExists($weightStr) == 0) // wave reference is not valid, |
---|
820 | Make/D/N=5 $weightStr,$zStr |
---|
821 | Wave w5 = $weightStr |
---|
822 | Wave z5 = $zStr // wave references to pass |
---|
823 | Make5GaussPoints(w5,z5) |
---|
824 | else |
---|
825 | if(exists(weightStr) > 1) |
---|
826 | Abort "wave name is already in use" //executed only if name is in use elsewhere |
---|
827 | endif |
---|
828 | Wave w5 = $weightStr |
---|
829 | Wave z5 = $zStr // create the wave references |
---|
830 | endif |
---|
831 | |
---|
832 | // current x point is the q-value for evaluation |
---|
833 | // |
---|
834 | // * for the input x, the resolution function waves are interpolated to get the correct values for |
---|
835 | // sigq, qbar and shad - since the model x-spacing may not be the same as |
---|
836 | // the experimental QSIG data. This is always the case when curve fitting, since fit_wave is |
---|
837 | // Igor-defined as 200 points and has its own (linear) q-(x)-scaling which will be quite different |
---|
838 | // from experimental data. |
---|
839 | // **note** if the (x) passed in is the experimental q-values, these values are |
---|
840 | // returned from the interpolation (as expected) |
---|
841 | |
---|
842 | i_shad = interp(x,qvals,shad) |
---|
843 | i_qbar = interp(x,qvals,qbar) |
---|
844 | i_sigq = interp(x,qvals,sigq) |
---|
845 | |
---|
846 | if (i_sigq>=0) |
---|
847 | // set up the integration |
---|
848 | // number of Gauss Quadrature points |
---|
849 | nord = 5 |
---|
850 | |
---|
851 | // end points of integration |
---|
852 | // limits are technically 0-inf, but wisely choose interesting region of q where R() is nonzero |
---|
853 | // +/- 3 sigq catches 99.73% of distrubution |
---|
854 | // change limits (and spacing of zi) at each evaluation based on R() |
---|
855 | //integration from va to vb |
---|
856 | va = -3*i_sigq + i_qbar |
---|
857 | if (va<0) |
---|
858 | va=0 //to avoid numerical error when va<0 (-ve q-value) |
---|
859 | endif |
---|
860 | vb = 3*i_sigq + i_qbar |
---|
861 | // Using 5 Gauss points |
---|
862 | // remember to index from 0,size-1 |
---|
863 | |
---|
864 | summ = 0.0 // initialize integral |
---|
865 | ii=0 // loop counter |
---|
866 | do |
---|
867 | // calculate Gauss points on integration interval (q-value for evaluation) |
---|
868 | zi = ( z5[ii]*(vb-va) + vb + va )/2.0 |
---|
869 | // calculate resolution function at input q-value (use the interpolated values and zi) |
---|
870 | Resoln = i_shad/sqrt(2*pi*i_sigq*i_sigq) |
---|
871 | Resoln *= exp((-1*(zi - i_qbar)^2)/(2*i_sigq*i_sigq)) |
---|
872 | //calculate partial sum for the passed-in model function |
---|
873 | yyy = w5[ii] * Resoln * fcn(w,zi) |
---|
874 | summ += yyy //add to the running total of the quadrature |
---|
875 | ii+=1 |
---|
876 | while (ii<nord) // end of loop over quadrature points |
---|
877 | |
---|
878 | // calculate value of integral to return |
---|
879 | answer = (vb-va)/2.0*summ |
---|
880 | // all scaling, background addition... etc. is done in the model calculation |
---|
881 | |
---|
882 | else |
---|
883 | // smear with the USANS routine |
---|
884 | // Make global string and local variables |
---|
885 | // now data folder aware, necessary for GlobalFit = FULL path to wave |
---|
886 | String/G gTrap_coefStr = GetWavesDataFolder(w, 2 ) |
---|
887 | Variable maxiter=10, tol=1e-2 |
---|
888 | |
---|
889 | // set up limits for the integration |
---|
890 | va=0 |
---|
891 | vb=abs(i_sigq) |
---|
892 | |
---|
893 | Variable/G gEvalQval = x |
---|
894 | |
---|
895 | // call qtrap to do actual work |
---|
896 | answer = qtrap_USANS(fcn,va,vb,tol,maxiter) |
---|
897 | answer /= vb |
---|
898 | |
---|
899 | endif |
---|
900 | |
---|
901 | Return (answer) |
---|
902 | |
---|
903 | End |
---|
904 | /////////////////////////////////////////////////////////////// |
---|
905 | |
---|
906 | Function GenericQuadrature_proto(w,x,dum) |
---|
907 | Wave w |
---|
908 | Variable x,dum |
---|
909 | |
---|
910 | Print "in GenericQuadrature_proto function" |
---|
911 | return(1) |
---|
912 | end |
---|
913 | |
---|
914 | // prototype function for smearing routines, Smear_Model_N |
---|
915 | // and trapzd_USANS() and qtrap_USANS() |
---|
916 | // it intentionally does nothing |
---|
917 | Function SANSModel_proto(w,x) |
---|
918 | Wave w |
---|
919 | Variable x |
---|
920 | |
---|
921 | Print "in SANSModel_proto function" |
---|
922 | return(1) |
---|
923 | end |
---|
924 | |
---|
925 | //Numerical Recipes routine to calculate the nn(th) stage |
---|
926 | //refinement of a trapezoid integration |
---|
927 | // |
---|
928 | //must be called sequentially from nn=1...n from qtrap() |
---|
929 | // to cumulatively refine the integration value |
---|
930 | // |
---|
931 | // in the conversion: |
---|
932 | // -- s was replaced with sVal and declared global (rather than static) |
---|
933 | // so that the nn-1 value would be available during the nn(th) call |
---|
934 | // |
---|
935 | // -- the specific coefficient wave for func() is passed in as a |
---|
936 | // global string (then converted to a wave reference), since |
---|
937 | // func() will eventually call sphereForm() |
---|
938 | // |
---|
939 | Function trapzd_USANS(fcn,aa,bb,nn) |
---|
940 | FUNCREF SANSModel_proto fcn |
---|
941 | Variable aa,bb,nn |
---|
942 | |
---|
943 | Variable xx,tnm,summ,del |
---|
944 | Variable it,jj,arg1,arg2 |
---|
945 | NVAR sVal=sVal //calling function must initialize this global |
---|
946 | NVAR qval = gEvalQval |
---|
947 | SVAR cwStr = gTrap_CoefStr //pass in the coefficient wave (string) |
---|
948 | Wave cw=$cwStr |
---|
949 | Variable temp=0 |
---|
950 | if(nn==1) |
---|
951 | arg1 = qval^2 + aa^2 |
---|
952 | arg2 = qval^2 + bb^2 |
---|
953 | temp = 0.5*(bb-aa)*(fcn(cw,sqrt(arg1)) + fcn(cw,sqrt(arg2))) |
---|
954 | sval = temp |
---|
955 | return(sVal) |
---|
956 | else |
---|
957 | it=1 |
---|
958 | it= 2^(nn-2) //done in NR with a bit shift <<= |
---|
959 | tnm = it |
---|
960 | del = (bb - aa)/tnm //this is the spacing of points to add |
---|
961 | xx = aa+0.5*del |
---|
962 | summ=0 |
---|
963 | for(jj=1;jj<=it;jj+=1) |
---|
964 | arg1 = qval^2 + xx^2 |
---|
965 | summ += fcn(cw,sqrt(arg1)) |
---|
966 | xx += del |
---|
967 | endfor |
---|
968 | sval = 0.5*(sval+(bb-aa)*summ/tnm) //replaces sval with its refined value |
---|
969 | return (sval) |
---|
970 | endif |
---|
971 | |
---|
972 | End |
---|
973 | |
---|
974 | |
---|
975 | // Numerical Recipes routine to calculate the integral of a |
---|
976 | // specified function, trapezoid rule is used to a user-specified |
---|
977 | // level of refinement using sequential calls to trapzd() |
---|
978 | // |
---|
979 | // in NR, eps and maxIt were global, pass them in here... |
---|
980 | // eps typically 1e-5 |
---|
981 | // maxit typically 20 |
---|
982 | Function qtrap_USANS(fcn,aa,bb,eps,maxIt) |
---|
983 | FUNCREF SANSModel_proto fcn |
---|
984 | Variable aa,bb,eps,maxit |
---|
985 | |
---|
986 | Variable/G sVal=0 //create and initialize what trapzd will return |
---|
987 | Variable jj,ss,olds |
---|
988 | |
---|
989 | olds = -1e30 //any number that is not the avg. of endpoints of the funciton |
---|
990 | for(jj=1;jj<=maxit;jj+=1) //call trapzd() repeatedly until convergence |
---|
991 | ss = trapzd_USANS(fcn,aa,bb,jj) |
---|
992 | if( abs(ss-olds) < eps*abs(olds) ) // good enough, stop now |
---|
993 | return ss |
---|
994 | endif |
---|
995 | if( (ss == 0.0) && (olds == 0.0) && (jj>6) ) //no progress? |
---|
996 | return ss |
---|
997 | endif |
---|
998 | olds = ss |
---|
999 | endfor |
---|
1000 | |
---|
1001 | Print "Maxit exceeded in qtrap. If you're here, there was an error in qtrap" |
---|
1002 | return(ss) //should never get here if function is well-behaved |
---|
1003 | |
---|
1004 | End |
---|
1005 | |
---|
1006 | Proc RRW() |
---|
1007 | ResetResolutionWaves() |
---|
1008 | End |
---|
1009 | |
---|
1010 | //utility procedures that are currently untied to any actions, although useful... |
---|
1011 | Proc ResetResolutionWaves(str) |
---|
1012 | String Str |
---|
1013 | Prompt str,"Pick the intensity wave with the resolution you want",popup,WaveList("*_i",";","") |
---|
1014 | |
---|
1015 | String/G root:gQvals = str[0,strlen(str)-3]+"_q" |
---|
1016 | String/G root:gSig_Q = str[0,strlen(str)-3]+"sq" |
---|
1017 | String/G root:gQ_bar = str[0,strlen(str)-3]+"qb" |
---|
1018 | String/G root:gShadow = str[0,strlen(str)-3]+"fs" |
---|
1019 | |
---|
1020 | |
---|
1021 | |
---|
1022 | //touch everything to make sure that the dependencies are |
---|
1023 | //properly updated - especially the $gQvals reference in the |
---|
1024 | // dependency assignment |
---|
1025 | fKillDependencies("Smear*") |
---|
1026 | |
---|
1027 | //replace the q-values and intensity (so they're the right length) |
---|
1028 | fResetSmearedModels("Smear*",root:gQvals) |
---|
1029 | |
---|
1030 | fRestoreDependencies("Smear*") |
---|
1031 | End |
---|
1032 | |
---|
1033 | // pass "*" as the matchString to do ALL dependent waves |
---|
1034 | // or "abc*" to get just the matching waves |
---|
1035 | // |
---|
1036 | Function fKillDependencies(matchStr) |
---|
1037 | String matchStr |
---|
1038 | |
---|
1039 | String str=WaveList(matchStr, ";", "" ),item,formula |
---|
1040 | Variable ii |
---|
1041 | |
---|
1042 | for(ii=0;ii<ItemsInList(str ,";");ii+=1) |
---|
1043 | item = StringFromList(ii, str ,";") |
---|
1044 | formula = GetFormula($item) |
---|
1045 | if(cmpstr("",formula)!=0) |
---|
1046 | Printf "wave %s had the formula %s removed\r",item,formula |
---|
1047 | Note $item, "FORMULA:"+formula |
---|
1048 | SetFormula $item, "" //clears the formula |
---|
1049 | endif |
---|
1050 | endfor |
---|
1051 | return(0) |
---|
1052 | end |
---|
1053 | |
---|
1054 | // pass "*" as the matchString to do ALL dependent waves |
---|
1055 | // or "abc*" to get just the matching waves |
---|
1056 | // |
---|
1057 | Function fRestoreDependencies(matchStr) |
---|
1058 | String matchStr |
---|
1059 | |
---|
1060 | String str=WaveList(matchStr, ";", "" ),item,formula |
---|
1061 | Variable ii |
---|
1062 | |
---|
1063 | for(ii=0;ii<ItemsInList(str ,";");ii+=1) |
---|
1064 | item = StringFromList(ii, str ,";") |
---|
1065 | formula = StringByKey("FORMULA", note($item),":",";") |
---|
1066 | if(cmpstr("",formula)!=0) |
---|
1067 | Printf "wave %s had the formula %s restored\r",item,formula |
---|
1068 | Note/K $item |
---|
1069 | SetFormula $item, formula //restores the formula |
---|
1070 | endif |
---|
1071 | endfor |
---|
1072 | return(0) |
---|
1073 | end |
---|
1074 | |
---|
1075 | Function fResetSmearedModels(matchStr,qStr) |
---|
1076 | String matchStr,qStr |
---|
1077 | |
---|
1078 | Duplicate/O $qStr root:smeared_qvals |
---|
1079 | |
---|
1080 | String str=WaveList(matchStr, ";", "" ),item,formula |
---|
1081 | Variable ii |
---|
1082 | |
---|
1083 | for(ii=0;ii<ItemsInList(str ,";");ii+=1) |
---|
1084 | item = StringFromList(ii, str ,";") |
---|
1085 | formula = StringByKey("FORMULA", note($item),":",";") |
---|
1086 | if(cmpstr("",formula)!=0) |
---|
1087 | Printf "wave %s has been duplicated to gQvals\r",item |
---|
1088 | Duplicate/O $qStr $item |
---|
1089 | Note $item, "FORMULA:"+formula //be sure to keep the formula note |
---|
1090 | Print " and the formula is",formula |
---|
1091 | endif |
---|
1092 | endfor |
---|
1093 | return(0) |
---|
1094 | end |
---|