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