Context Navigation

GaussUtils_v40.ipf @ 794

Last change on this file since 794 was 794, checked in by srkline, 12 years ago

Lots of changes:
-2D resolution smearing
-2D error propagation

1) 2D resolution smearing has been corrected to use sigma (perp) correctly
rather than phi. This shows up in the quadrature loop in all of the 2D models
and in the Smear_2D "generic" function.

2) 1D resolutionn smearing is now corrected to account for the integration
range of +/- 3 sigma (99.73% of distribution). Now the result is divided by
0.9973 to rescale it to the correct value.

3) Smeared models are now AAO to improve speed and to allow easier use with
functions that are inherently AAO. No changes are needed, since the call is
behind the scenes, replacing Smear_Model_N() with Smear_Model_N_AAO().

4) in PlotUtils_2D added functions to re-bin the QxQy? data into a 1D format
BinQxQy_to_1D(). This also re-bins the errors in two ways, adding the per-pixel
errors in quadrature, or the deviation from the mean of the intensity. Manually
editing the intensity declaration allows 2D->1D binning of smeared models.

5) Per-pixel propagation of errors has been added through the entire data
reduction sequence. Data errors are generated on loading using Poisson
statistics (specifically tailored for accuracy at low counts) and then is
propagated through each manipulation of the data using standard error
propagation. The error matrix is linear_data_error. As a by-product, all of
the math operations on data are explicitly done on linear_data, to avoid
any potential mistakes of log/linear scaling. Results of this propagation
largely match J. Barker's /ERROR routines from the VAX, with some differences
at low data count values (as expected) and at higher count values near the
beam stop (somewhat unexpected). This per-pixel error is ouput in the QxQy_ASCII
data files. NO CHANGE has been made to the 1D data, which uses the deviation from
the mean as the error - since this is correct.

6) Added tables for the uncertainty in attenuator transmission (from JGB)

7) added two more REAL values to the VAX header to store information
necessary for error propagation. These are couting error that are part of
the transmission error and of the absolute scaling error. These add Read/Write?
functions in NCNR_DataReadWrite

The transmission error (one standard deviation) is written at byte 396 (4 bytes)

RealsRead?[41]

The Box Sum error (one standard deviation) is written at byte 400 (4 bytes)

RealsRead?[42]

File size: 50.4 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format