Context Navigation

V_DetectorBinning.ipf @ 947

Last change on this file since 947 was 947, checked in by srkline, 8 years ago

Changes to VCALC files, moving detector information into sub-folders in anticipation of the possible folder structure that I will use for the actual data reduction. But purely speculative at this point.

Fixing some typos in other procedures.

File size: 56.6 KB

Line
1	#pragma rtGlobals=1 // Use modern global access method.
2
3
4	///// Procedures for:
5	// Generating the detector panels
6	// Filling the panels with Qtot, QxQyQz
7	// Filling the "data" with a model function
8	// Averaging the panels (independently) into I(Q)
9	//
10	//
11	//There are some things in the current circular averaging that don't make any sense
12	//and don't seem to really do anything at all, so trim them out.
13	//1) subdividing pixels near the beam stop into 9 sub-pixels
14	//2) non-linear correction (maybe keep this as a dummy). I may need to put some sort
15	//of correction here for the offset, or it may already be done as I receive it.
16	//
17	//
18	//
19	//Do I separate out the circular, sector, rectangular, annular averaging into
20	//separate routines?
21	//
22	//
23
24
25	Proc PlotFrontPanels()
26	fPlotFrontPanels()
27	End
28
29	// to plot I(q) for the 4 front panels
30	//
31	// *** Call this function when front panels are adjusted, or wavelength, etc. changed
32	//
33	Function fPlotFrontPanels()
34
35	// space is allocated for all of the detectors and Q's on initialization
36	// calculate Qtot, qxqyqz arrays from geometry
37	V_CalculateQFrontPanels()
38
39	// fill the panels with fake sphere scattering data
40	// TODO: am I in the right data folder??
41	SetDataFolder root:Packages:NIST:VSANS:VCALC:Front
42
43	WAVE det_FL = det_FL
44	WAVE det_FR = det_FR
45	WAVE det_FT = det_FT
46	WAVE det_FB = det_FB
47
48	WAVE qTot_FL = qTot_FL
49	WAVE qTot_FR = qTot_FR
50	WAVE qTot_FT = qTot_FT
51	WAVE qTot_FB = qTot_FB
52
53	FillPanel_wModelData(det_FL,qTot_FL,"FL")
54	FillPanel_wModelData(det_FR,qTot_FR,"FR")
55	FillPanel_wModelData(det_FT,qTot_FT,"FT")
56	FillPanel_wModelData(det_FB,qTot_FB,"FB")
57	// det_FL = V_SphereForm(1,60,1e-6,1,qTot_FL[p][q])
58	// det_FR = V_SphereForm(1,60,1e-6,1,qTot_FR[p][q])
59	// det_FT = V_SphereForm(1,60,1e-6,1,qTot_FT[p][q])
60	// det_FB = V_SphereForm(1,60,1e-6,1,qTot_FB[p][q])
61
62	SetDataFolder root:
63
64	// set any "shadowed" area of the T/B detectors to NaN to get a realitic
65	// view of how much of the detectors are actually collecting data
66	// -- I can get the separation L/R from the panel - only this "open" width is visible.
67	//TODO - make this a proper shadow - TB extent of the LR panels matters too, not just the LR separation
68	V_SetShadow_TopBottom("Front","FT") // TODO: -- be sure the data folder is properly set (within the function...)
69	V_SetShadow_TopBottom("Front","FB")
70
71	// do the q-binning for each of the panels to get I(Q)
72	Execute "BinAllFrontPanels()"
73
74	// plot the results
75	Execute "Front_IQ_Graph()"
76	Execute "FrontPanels_AsQ()"
77	End
78
79	// TODO: hard wired for a sphere - change this to allow minimal selections and altering of coefficients
80	// TODO: add the "fake" 2D simulation to fill the panels which are then later averaged as I(Q)
81	Function FillPanel_wModelData(det,qTot,type)
82	Wave det,qTot
83	String type
84
85	SetDataFolder root:Packages:NIST:VSANS:VCALC:Front
86
87	// q-values and detector arrays already allocated and calculated
88	Duplicate/O det tmpInten,tmpSig,prob_i
89
90	Variable imon,trans,thick,sdd,pixSizeX,pixSizeY,sdd_offset
91
92	//imon = V_BeamIntensity()*CountTime
93	imon = VCALC_getImon() //TODO: currently from the panel, not calculated
94	trans = 0.8
95	thick = 0.1
96
97	// need SDD
98	// need pixel dimensions
99	// nominal sdd in meters, offset in mm, want result in cm !
100	sdd = VCALC_getSDD(type)*100 + VSANS_getTopBottomSDDOffset(type) / 10 // result is sdd in [cm]
101
102	pixSizeX = VCALC_getPixSizeX(type) // cm
103	pixSizeY = VCALC_getPixSizeY(type)
104
105
106	//?? pick the function from a popup on the panel? (bypass the analysis panel, or maybe it's better to
107	// keep the panel to keep people used to using it.)
108	String funcStr = VCALC_getModelFunctionStr()
109	strswitch(funcStr)
110	case "Big Debye":
111	tmpInten = V_Debye(10,3000,0.0001,qTot[p][q])
112	break
113	case "Big Sphere":
114	tmpInten = V_SphereForm(1,900,1e-6,0.01,qTot[p][q])
115	break
116	case "Debye":
117	tmpInten = V_Debye(10,300,0.0001,qTot[p][q])
118	break
119	case "Sphere":
120	tmpInten = V_SphereForm(1,60,1e-6,0.001,qTot[p][q])
121	break
122	default:
123	tmpInten = V_Debye(10,300,0.1,qTot[p][q])
124	endswitch
125	// tmpInten = V_SphereForm(1,100,1e-6,1,qTot[p][q])
126
127	// tmpInten = V_Debye(scale,rg,bkg,x)
128	// tmpInten = V_Debye(10,300,0.1,qTot[p][q])
129
130
131	///////////////
132	// // calculate the scattering cross section simply to be able to estimate the transmission
133	// Variable sig_sas=0
134	//
135	// // remember that the random deviate is the coherent portion ONLY - the incoherent background is
136	// // subtracted before the calculation.
137	// CalculateRandomDeviate(funcUnsmeared,$coefStr,wavelength,"root:Packages:NIST:SAS:ran_dev",sig_sas)
138	//
139	// if(sig_sas > 100)
140	// DoAlert 0,"SAS cross section > 100. Estimates of multiple scattering are unreliable. Choosing a model with a well-defined Rg may help"
141	// endif
142	//
143	// // calculate the multiple scattering fraction for display (10/2009)
144	// Variable ii,nMax=10,tau
145	// mScat=0
146	// tau = thick*sig_sas
147	// // this sums the normalized scattering P', so the result is the fraction of multiply coherently scattered
148	// // neutrons out of those that were scattered
149	// for(ii=2;ii<nMax;ii+=1)
150	// mScat += tau^(ii)/factorial(ii)
151	//// print tau^(ii)/factorial(ii)
152	// endfor
153	// estTrans = exp(-1thicksig_sas) //thickness and sigma both in units of cm
154	// mscat *= (estTrans)/(1-estTrans)
155	//
156	//// if(mScat > 0.1) // Display warning
157	//
158	// Print "Sig_sas = ",sig_sas
159	////////////////////
160
161	prob_i = transthickpixSizeXpixSizeY/(sdd)^2tmpInten //probability of a neutron in q-bin(i)
162
163	tmpInten = (imon)*prob_i //tmpInten is not the model calculation anymore!!
164
165
166	/// **** can I safely assume a Gaussian error in the count rate??
167	tmpSig = sqrt(tmpInten) // corrected based on John's memo, from 8/9/99
168
169	tmpInten += gnoise(tmpSig)
170	tmpInten = (tmpInten[p][q] < 0) ? 0 : tmpInten[p][q] // MAR 2013 -- is this the right thing to do
171	tmpInten = trunc(tmpInten)
172
173
174	det = tmpInten
175
176	// if I want "absolute" scale -- then I lose the integer nature of the detector (but keep the random)
177	det /= transthickpixSizeXpixSizeY/(sdd)^2imon
178
179	KillWaves/Z tmpInten,tmpSig,prob_i
180	SetDataFolder root:
181
182	return(0)
183	End
184
185
186
187	// works for Left, works for Right... works for T/B too.
188	//
189	// - TODO: be sure that the Q's are calculated correctly even when the beam is off of the
190	// detector, and on different sides of the detector (or T/B) - since it will be in a different
191	// relative postion to 0,0 on the detector. If the postions are symmetric, then the Q's should be identical.
192	// --- test this...
193	// TODO -- be sure I'm in the right data folder. nothing is set correctly right now
194	//
195	// TODO: make all detector parameters global, not hard-wired
196	//
197	//
198	// --- Panels are all allocated in the initialization. Here, only the q-values are calculated
199	// when anything changes
200	//
201	Function V_CalculateQFrontPanels()
202
203	Variable xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY
204	Variable F_LR_sep,F_TB_sep,F_offset,F_sdd_offset
205
206	// get the values from the panel + constants
207	F_LR_sep = VCALC_getPanelSeparation("FLR")
208	F_TB_sep = VCALC_getPanelSeparation("FTB")
209	F_offset = VCALC_getLateralOffset("FL")
210
211	SDD = VCALC_getSDD("FL") //nominal SDD - need offset for TB
212	lam = VCALC_getWavelength()
213
214	//separations are in mm -- need to watch the units, convert to cm
215	F_LR_sep /= 10
216	F_TB_sep /= 10
217	// TODO - I'm treating the separation as the TOTAL width - so the difference
218	// from the "center" to the edge is 1/2 of the separation
219
220	// TODO (make the N along the tube length a variable, since this can be reset @ acquisition)
221
222	F_sdd_offset = VSANS_getTopBottomSDDOffset("FT") //T/B are 300 mm farther back //TODO: make all detector parameters global, not hard-wired
223
224	SetDataFolder root:Packages:NIST:VSANS:VCALC:Front
225	Wave det_FL,det_FR // these are (48,256)
226	Wave det_FT,det_FB // these are (128,48)
227
228	//FRONT/LEFT
229	WAVE qTot_FL,qx_FL,qy_FL,qz_FL
230	qTot_FL = 0
231	qx_FL = 0
232	qy_FL = 0
233	qz_FL = 0
234
235	// TODO - these are to be set from globals, not hard-wired. N and pixelSixze will be known (or pre-measured)
236	// pixel sizes are in cm
237	pixSizeX = VCALC_getPixSizeX("FL")
238	pixSizeY = VCALC_getPixSizeY("FL")
239	// pixSizeX = 0.8 // 0.8 cm/pixel along width
240	// pixSizeY = 0.4 // approx 0.4 cm/pixel along length
241
242	xCtr = 48+(F_LR_sep/2/pixSizeX) // TODO -- check -- starting from 47 rather than 48 (but I'm in pixel units for centers)??
243	yCtr = 127
244	V_Detector_2Q(det_FL,qTot_FL,qx_FL,qy_FL,qz_FL,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
245	// Print "xy for FL = ",xCtr,yCtr
246
247	//set the wave scaling for the detector image so that it can be plotted in q-space
248	// TODO: this is only approximate - since the left "edge" is not the same from top to bottom, so I crudely
249	// take the middle value. At very small angles, OK, at 1m, this is a crummy approximation.
250	// since qTot is magnitude only, I need to put in the (-ve)
251	SetScale/I x WaveMin(qx_FL),WaveMax(qx_FL),"", det_FL //this sets the left and right ends of the data scaling
252	SetScale/I y WaveMin(qy_FL),WaveMax(qy_FL),"", det_FL
253	//////////////////
254
255	//FRONT/RIGHT
256	SetDataFolder root:Packages:NIST:VSANS:VCALC:Front
257	WAVE qTot_FR,qx_FR,qy_FR,qz_FR
258	qTot_FR = 0
259	qx_FR = 0
260	qy_FR = 0
261	qz_FR = 0
262
263	// TODO - these are to be set from globals, not hard-wired
264	// pixel sizes are in cm
265	pixSizeX = VCALC_getPixSizeX("FR")
266	pixSizeY = VCALC_getPixSizeY("FR")
267
268	xCtr = -(F_LR_sep/2/pixSizeX)-1
269	yCtr = 127
270	V_Detector_2Q(det_FR,qTot_FR,qx_FR,qy_FR,qz_FR,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
271	// Print "xy for FR = ",xCtr,yCtr
272	SetScale/I x WaveMin(qx_FR),WaveMax(qx_FR),"", det_FR //this sets the left and right ends of the data scaling
273	SetScale/I y WaveMin(qy_FR),WaveMax(qy_FR),"", det_FR
274	/////////////////
275
276	//FRONT/TOP
277	SetDataFolder root:Packages:NIST:VSANS:VCALC:Front
278	WAVE qTot_FT,qx_FT,qy_FT,qz_FT
279	qTot_FT = 0
280	qx_FT = 0
281	qy_FT = 0
282	qz_FT = 0
283
284	// TODO - these are to be set from globals, not hard-wired
285	// pixel sizes are in cm
286	pixSizeX = VCALC_getPixSizeX("FT")
287	pixSizeY = VCALC_getPixSizeY("FT")
288
289	xCtr = 64
290	yCtr = -(F_TB_sep/2/pixSizeY)-1
291	// global sdd_offset is in (mm), convert to meters here for the Q-calculation
292	V_Detector_2Q(det_FT,qTot_FT,qx_FT,qy_FT,qz_FT,xCtr,yCtr,sdd+F_sdd_offset/1000,lam,pixSizeX,pixSizeY)
293	// Print "xy for FT = ",xCtr,yCtr
294	SetScale/I x WaveMin(qx_FT),WaveMax(qx_FT),"", det_FT //this sets the left and right ends of the data scaling
295	SetScale/I y WaveMin(qy_FT),WaveMax(qy_FT),"", det_FT
296	//////////////////
297
298	//FRONT/BOTTOM
299	SetDataFolder root:Packages:NIST:VSANS:VCALC:Front
300	WAVE qTot_FB,qx_FB,qy_FB,qz_FB
301	qTot_FB = 0
302	qx_FB = 0
303	qy_FB = 0
304	qz_FB = 0
305
306	// TODO - these are to be set from globals, not hard-wired
307	// pixel sizes are in cm
308	pixSizeX = VCALC_getPixSizeX("FB")
309	pixSizeY = VCALC_getPixSizeY("FB")
310
311	xCtr = 64
312	yCtr = 48+(F_TB_sep/2/pixSizeY) // TODO -- check -- starting from 47 rather than 48 (but I'm in pixel units for centers)??
313	// global sdd_offset is in (mm), convert to meters here for the Q-calculation
314	V_Detector_2Q(det_FB,qTot_FB,qx_FB,qy_FB,qz_FB,xCtr,yCtr,sdd+F_sdd_offset/1000,lam,pixSizeX,pixSizeY)
315	// Print "xy for FB = ",xCtr,yCtr
316	SetScale/I x WaveMin(qx_FB),WaveMax(qx_FB),"", det_FB //this sets the left and right ends of the data scaling
317	SetScale/I y WaveMin(qy_FB),WaveMax(qy_FB),"", det_FB
318	/////////////////
319
320	SetDataFolder root:
321
322	return(0)
323	End
324
325
326	// TODO" - this doesn't quite mask things out as they should be (too much is masked L/R of center)
327	// and the outer edges of the detector are masked even if the TB panels extend past the TB of the LR panels.
328	// ? skip the masking? but then I bin the detector data directly to get I(q), skipping the masked NaN values...
329	//
330	Function V_SetShadow_TopBottom(folderStr,type)
331	String folderStr,type
332
333	Variable LR_sep,nPix,xCtr,ii,jj,numCol,pixSizeX,pixSizeY
334
335	/// !! type passed in will be FT, FB, MT, MB, so I can't ask for the panel separation -- or I'll get the TB separation...
336	if(cmpstr(type[0],"F")==0)
337	//front
338	ControlInfo/W=VCALC VCALCCtrl_2a
339	LR_sep = V_Value
340	else
341	//middle
342	ControlInfo/W=VCALC VCALCCtrl_3a
343	LR_sep = V_Value
344	endif
345	//separations on panel are in mm -- need to watch the units, convert to cm
346	LR_sep /= 10
347
348	//detector data
349	Wave det = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+":det_"+type)
350
351	// TODO - these are to be set from globals, not hard-wired
352	// pixel sizes are in cm for T/B detector
353	// TODO: the "FT" check is hard wired for FRONT -- get rid of this...
354
355	pixSizeX = VCALC_getPixSizeX(type)
356	pixSizeY = VCALC_getPixSizeY(type)
357
358	//TODO -- get this from a global
359	xCtr = 64
360	nPix = (LR_sep/2/pixSizeX) // # of pixels Left/right of center that are not obscured by L/R panels
361
362	numCol = DimSize(det,0) // x dim (columns)
363	for(ii=0;ii<(xCtr-nPix-2);ii+=1)
364	det[ii][] = NaN
365	endfor
366	for(ii=(xCtr+nPix+2);ii<numCol;ii+=1)
367	det[ii][] = NaN
368	endfor
369
370	return(0)
371	end
372
373
374	Window FrontPanels_AsQ() : Graph
375	// DoWindow/F FrontPanels_AsQ
376	// if(V_flag == 0)
377	// PauseUpdate; Silent 1 // building window...
378	// Display /W=(1477,44,1978,517)
379
380	SetDataFolder root:Packages:NIST:VSANS:VCALC:Front
381
382	CheckDisplayed /W=VCALC#Panels_Q det_FB
383	if(V_flag == 0)
384	AppendImage/W=VCALC#Panels_Q det_FB
385	ModifyImage/W=VCALC#Panels_Q det_FB ctab= {,,ColdWarm,0}
386	AppendImage/W=VCALC#Panels_Q det_FT
387	ModifyImage/W=VCALC#Panels_Q det_FT ctab= {,,ColdWarm,0}
388	AppendImage/W=VCALC#Panels_Q det_FL
389	ModifyImage/W=VCALC#Panels_Q det_FL ctab= {,,ColdWarm,0}
390	AppendImage/W=VCALC#Panels_Q det_FR
391	ModifyImage/W=VCALC#Panels_Q det_FR ctab= {,,ColdWarm,0}
392	endif
393
394	Variable dval
395	ControlInfo/W=VCALC setVar_b
396	dval = V_Value
397
398	SetAxis/W=VCALC#Panels_Q left -dval,dval
399	SetAxis/W=VCALC#Panels_Q bottom -dval,dval
400
401	ControlInfo/W=VCALC check_0a
402	// V_Value == 1 if checked
403	ModifyImage/W=VCALC#Panels_Q det_FB log=V_Value
404	ModifyImage/W=VCALC#Panels_Q det_FT log=V_Value
405	ModifyImage/W=VCALC#Panels_Q det_FL log=V_Value
406	ModifyImage/W=VCALC#Panels_Q det_FR log=V_Value
407
408
409	SetDataFolder root:
410
411	// ModifyGraph width={Aspect,1},height={Aspect,1},gbRGB=(56797,56797,56797)
412	// ModifyGraph grid=2
413	// ModifyGraph mirror=2
414	// SetAxis left -0.2,0.2
415	// SetAxis bottom -0.2,0.2
416	// endif
417	EndMacro
418
419
420	// For a given detector panel, calculate the q-values
421	// -work with everything as arrays
422	// Input needed:
423	// detector data
424	// detector type (LRTB?)
425	// beam center (may be off the detector)
426	// SDD
427	// lambda
428	//
429	// pixel dimensions for detector type (global constants)
430	// - data dimensions read directly from array
431	//
432	// --What is calculated:
433	// array of Q
434	// array of qx,qy,qz
435	// array of error already exists
436	//
437	//
438	// -- sdd in meters
439	// -- lambda in Angstroms
440	Function V_Detector_2Q(data,qTot,qx,qy,qz,xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY)
441	Wave data,qTot,qx,qy,qz
442	Variable xCtr,yCtr,sdd,lam,pixSizeX,pixSizeY
443
444	// loop over the array and calculate the values - this is done as a wave assignment
445	// TODO -- be sure that it's p,q -- or maybe p+1,q+1 as used in WriteQIS.ipf
446	qTot = V_CalcQval(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
447	qx = V_CalcQX(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
448	qy = V_CalcQY(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
449	qz = V_CalcQZ(p,q,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
450
451	return(0)
452	End
453
454
455	//////////////////////
456	// NOTE: The Q calculations are different than what is in GaussUtils in that they take into
457	// accout the different x/y pixel sizes and the beam center not being on the detector -
458	// off a different edge for each LRTB type
459	/////////////////////
460
461	//function to calculate the overall q-value, given all of the necesary trig inputs
462	//NOTE: detector locations passed in are pixels = 0.5cm real space on the detector
463	//and are in detector coordinates (1,128) rather than axis values
464	//the pixel locations need not be integers, reals are ok inputs
465	//sdd is in meters
466	//wavelength is in Angstroms
467	//
468	//returned magnitude of Q is in 1/Angstroms
469	//
470	Function V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
471	Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
472
473	Variable dx,dy,qval,two_theta,dist
474
475	sdd *=100 //convert to cm
476	dx = (xaxval - xctr)*pixSizeX //delta x in cm
477	dy = (yaxval - yctr)*pixSizeY //delta y in cm
478	dist = sqrt(dx^2 + dy^2)
479
480	two_theta = atan(dist/sdd)
481
482	qval = 4Pi/lamsin(two_theta/2)
483
484	return qval
485	End
486
487	//calculates just the q-value in the x-direction on the detector
488	//input/output is the same as CalcQval()
489	//ALL inputs are in detector coordinates
490	//
491	//NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector
492	//sdd is in meters
493	//wavelength is in Angstroms
494	//
495	// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
496	// now properly accounts for qz
497	//
498	Function V_CalcQX(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
499	Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
500
501	Variable qx,qval,phi,dx,dy,dist,two_theta
502
503	qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
504
505	sdd *=100 //convert to cm
506	dx = (xaxval - xctr)*pixSizeX //delta x in cm
507	dy = (yaxval - yctr)*pixSizeY //delta y in cm
508	phi = V_FindPhi(dx,dy)
509
510	//get scattering angle to project onto flat detector => Qr = qval*cos(theta)
511	dist = sqrt(dx^2 + dy^2)
512	two_theta = atan(dist/sdd)
513
514	qx = qvalcos(two_theta/2)cos(phi)
515
516	return qx
517	End
518
519	//calculates just the q-value in the y-direction on the detector
520	//input/output is the same as CalcQval()
521	//ALL inputs are in detector coordinates
522	//NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector
523	//sdd is in meters
524	//wavelength is in Angstroms
525	//
526	// repaired incorrect qx and qy calculation 3 dec 08 SRK (Lionel and C. Dewhurst)
527	// now properly accounts for qz
528	//
529	Function V_CalcQY(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
530	Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
531
532	Variable dy,qval,dx,phi,qy,dist,two_theta
533
534	qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
535
536	sdd *=100 //convert to cm
537	dx = (xaxval - xctr)*pixSizeX //delta x in cm
538	dy = (yaxval - yctr)*pixSizeY //delta y in cm
539	phi = V_FindPhi(dx,dy)
540
541	//get scattering angle to project onto flat detector => Qr = qval*cos(theta)
542	dist = sqrt(dx^2 + dy^2)
543	two_theta = atan(dist/sdd)
544
545	qy = qvalcos(two_theta/2)sin(phi)
546
547	return qy
548	End
549
550	//calculates just the z-component of the q-vector, not measured on the detector
551	//input/output is the same as CalcQval()
552	//ALL inputs are in detector coordinates
553	//NOTE: detector locations passed in are pixel = 0.5cm real space on the Ordela detector
554	//sdd is in meters
555	//wavelength is in Angstroms
556	//
557	// not actually used, but here for completeness if anyone asks
558	//
559	Function V_CalcQZ(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
560	Variable xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY
561
562	Variable dy,qval,dx,phi,qz,dist,two_theta
563
564	qval = V_CalcQval(xaxval,yaxval,xctr,yctr,sdd,lam,pixSizeX,pixSizeY)
565
566	sdd *=100 //convert to cm
567
568	//get scattering angle to project onto flat detector => Qr = qval*cos(theta)
569	dx = (xaxval - xctr)*pixSizeX //delta x in cm
570	dy = (yaxval - yctr)*pixSizeY //delta y in cm
571	dist = sqrt(dx^2 + dy^2)
572	two_theta = atan(dist/sdd)
573
574	qz = qval*sin(two_theta/2)
575
576	return qz
577	End
578
579	//phi is defined from +x axis, proceeding CCW around [0,2Pi]
580	Threadsafe Function V_FindPhi(vx,vy)
581	variable vx,vy
582
583	variable phi
584
585	phi = atan(vy/vx) //returns a value from -pi/2 to pi/2
586
587	// special cases
588	if(vx==0 && vy > 0)
589	return(pi/2)
590	endif
591	if(vx==0 && vy < 0)
592	return(3*pi/2)
593	endif
594	if(vx >= 0 && vy == 0)
595	return(0)
596	endif
597	if(vx < 0 && vy == 0)
598	return(pi)
599	endif
600
601
602	if(vx > 0 && vy > 0)
603	return(phi)
604	endif
605	if(vx < 0 && vy > 0)
606	return(phi + pi)
607	endif
608	if(vx < 0 && vy < 0)
609	return(phi + pi)
610	endif
611	if( vx > 0 && vy < 0)
612	return(phi + 2*pi)
613	endif
614
615	return(phi)
616	end
617
618	Function V_SphereForm(scale,radius,delrho,bkg,x)
619	Variable scale,radius,delrho,bkg
620	Variable x
621
622	// variables are:
623	//[0] scale
624	//[1] radius (A)
625	//[2] delrho (A-2)
626	//[3] background (cm-1)
627
628	// Variable scale,radius,delrho,bkg
629	// scale = w[0]
630	// radius = w[1]
631	// delrho = w[2]
632	// bkg = w[3]
633
634
635	// calculates scale * f^2/Vol where f=Vol3delrho*((sin(qr)-qrcos(qr))/qr^3
636	// and is rescaled to give [=] cm^-1
637
638	Variable bes,f,vol,f2
639	//
640	//handle q==0 separately
641	If(x==0)
642	f = 4/3piradius^3delrhodelrhoscale1e8 + bkg
643	return(f)
644	Endif
645
646	// bes = 3(sin(xradius)-xradiuscos(x*radius))/x^3/radius^3
647
648	bes = 3sqrt(pi/(2xradius))BesselJ(1.5,xradius)/(xradius)
649
650	vol = 4pi/3radius^3
651	f = volbesdelrho // [=] A
652	// normalize to single particle volume, convert to 1/cm
653	f2 = f * f / vol * 1.0e8 // [=] 1/cm
654
655	return (scale*f2+bkg) // Scale, then add in the background
656
657	End
658
659	Function V_Debye(scale,rg,bkg,x)
660	Variable scale,rg,bkg
661	Variable x
662
663	// variables are:
664	//[0] scale factor
665	//[1] radius of gyration [A]
666	//[2] background [cm-1]
667
668	// calculates (scale*debye)+bkg
669	Variable Pq,qr2
670
671	qr2=(x*rg)^2
672	Pq = 2*(exp(-(qr2))-1+qr2)/qr2^2
673
674	//scale
675	Pq *= scale
676	// then add in the background
677	return (Pq+bkg)
678	End
679
680
681
682
683
684
685
686
687
688	//
689	// these routines bin the 2D q data to 1D I(q). Currently the Qtot is magnitude only, no sign (since
690	// it's being binned to I(Q), having a sign makes no sense. If you want the sign, work from qxqyqz
691	//
692	// first - the DeltaQ step is set as the smaller detector resolution (along tube)
693	// which is different for LR / TB geometry. This is not set in stone.
694	//
695	// second - each detector is binned separately
696	//
697	// -- like the routines in CircSectAve, start with 500 points, and trim after binning is done.
698	// you'l end up with < 200 points.
699	//
700	// the results are in iBin_qxqy, qBin_qxqy, and eBin_qxqy, in the folder passed
701	//
702	Proc BinAllFrontPanels()
703
704	SetDeltaQ("Front","FL")
705	SetDeltaQ("Front","FT")
706
707	Variable binType
708	ControlInfo/W=VCALC popup_b
709	binType = V_Value // V_value counts menu items from 1, so 1=1, 2=2, 3=4
710
711	if(binType == 1)
712	V_BinQxQy_to_1D("","FL")
713	V_BinQxQy_to_1D("","FR")
714	V_BinQxQy_to_1D("","FT")
715	V_BinQxQy_to_1D("","FB")
716	endif
717
718	if(binType == 2)
719	V_BinQxQy_to_1D("","FLR")
720	V_BinQxQy_to_1D("","FTB")
721	endif
722
723	if(binType == 3)
724	V_BinQxQy_to_1D("","FLRTB")
725	endif
726
727	// TODO -- this is only a temporary fix for slit mode
728	if(binType == 4)
729	/// this is for a tall, narrow slit mode
730	V_fBinDetector_byRows("Front","FL")
731	V_fBinDetector_byRows("Front","FR")
732	V_fBinDetector_byRows("Front","FT")
733	V_fBinDetector_byRows("Front","FB")
734	endif
735
736	End
737
738
739	Function SetDeltaQ(folderStr,type)
740	String folderStr,type
741
742	WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + folderStr + ":det_"+type) // 2D detector data
743
744	Variable xDim,yDim,delQ
745
746	xDim=DimSize(inten,0)
747	yDim=DimSize(inten,1)
748
749	if(xDim<yDim)
750	WAVE qx = $("root:Packages:NIST:VSANS:VCALC:" + folderStr + ":qx_"+type)
751	delQ = abs(qx[0][0] - qx[1][0])/2
752	else
753	WAVE qy = $("root:Packages:NIST:VSANS:VCALC:" + folderStr + ":qy_"+type)
754	delQ = abs(qy[0][1] - qy[0][0])/2
755	endif
756
757	// set the global
758	Variable/G $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_"+type) = delQ
759	// Print "SET delQ = ",delQ," for ",type
760
761	return(0)
762	end
763
764
765	//TODO -- need a switch here to dispatch to the averaging type
766	Proc V_BinQxQy_to_1D(folderStr,type)
767	String folderStr
768	String type
769	// Prompt folderStr,"Pick the data folder containing 2D data",popup,getAList(4)
770	// Prompt type,"detector identifier"
771
772
773	V_fDoBinning_QxQy2D(folderStr, type)
774
775
776	/// this is for a tall, narrow slit mode
777	// V_fBinDetector_byRows(folderStr,type)
778
779	End
780
781
782	Proc V_Graph_1D_detType(folderStr,type)
783	String folderStr,type
784
785	SetDataFolder root:Packages:NIST:VSANS:VCALC
786
787	Display $("iBin_qxqy"+"_"+type) vs $("qBin_qxqy"+"_"+type)
788	ModifyGraph mirror=2,grid=1,log=1
789	ModifyGraph mode=4,marker=19,msize=2
790	// ErrorBars/T=0 iBin_qxqy Y,wave=(eBin2D_qxqy,eBin2D_qxqy) // for simulations, I don't have 2D uncertainty
791	ErrorBars/T=0 $("iBin_qxqy"+"_"+type) Y,wave=($("eBin_qxqy"+"_"+type),$("eBin_qxqy"+"_"+type))
792	legend
793
794	SetDataFolder root:
795
796	End
797
798
799	// see the equivalent function in PlotUtils2D_v40.ipf
800	//
801	//Function fDoBinning_QxQy2D(inten,qx,qy,qz)
802	//
803	// this has been modified to accept different detector panels and to take arrays
804	// -- type = FL or FR or...other panel identifiers
805	//
806	// TODO "iErr" is all messed up since it doesn't really apply here for data that is not 2D simulation
807	//
808	//
809	Function V_fDoBinning_QxQy2D(folderStr,type)
810	String folderStr,type
811
812	Variable nSets = 0
813	Variable xDim,yDim
814	Variable ii,jj
815	Variable qVal,nq,var,avesq,aveisq
816	Variable binIndex,val
817
818
819	SetDataFolder root:Packages:NIST:VSANS:VCALC
820
821	// now switch on the type to determine which waves to declare and create
822	// since there may be more than one panel to step through. There may be two, there may be four
823	//
824
825	strswitch(type) // string switch
826	case "FL": // execute if case matches expression
827	case "FR":
828	NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_FL")
829	WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+type) // 2D detector data
830	WAVE/Z iErr = $("iErr_"+type) // 2D errors -- may not exist, especially for simulation
831	Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+type) // 2D q-values
832	nSets = 1
833	break
834
835	case "FT":
836	case "FB":
837	NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_FT")
838	WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+type) // 2D detector data
839	WAVE/Z iErr = $("iErr_"+type) // 2D errors -- may not exist, especially for simulation
840	Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+type) // 2D q-values
841	nSets = 1
842	break
843
844	case "ML":
845	case "MR":
846	NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_ML")
847	WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+type) // 2D detector data
848	WAVE/Z iErr = $("iErr_"+type) // 2D errors -- may not exist, especially for simulation
849	Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+type) // 2D q-values
850	nSets = 1
851	break
852
853	case "MT":
854	case "MB":
855	NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_MT")
856	WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+type) // 2D detector data
857	WAVE/Z iErr = $("iErr_"+type) // 2D errors -- may not exist, especially for simulation
858	Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+type) // 2D q-values
859	nSets = 1
860	break
861
862	case "B":
863	NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_B")
864	WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Back" + ":det_"+type) // 2D detector data
865	WAVE/Z iErr = $("iErr_"+type) // 2D errors -- may not exist, especially for simulation
866	Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Back" +":qTot_"+type) // 2D q-values
867	nSets = 1
868	break
869
870	case "FLR":
871	NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_FL")
872	WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FL") // 2D detector data
873	WAVE/Z iErr = $("iErr_"+"FL") // 2D errors -- may not exist, especially for simulation
874	Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FL") // 2D q-values
875	WAVE inten2 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FR") // 2D detector data
876	WAVE/Z iErr2 = $("iErr_"+"FR") // 2D errors -- may not exist, especially for simulation
877	Wave qTotal2 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FR") // 2D q-values
878	nSets = 2
879	break
880
881	case "FTB":
882	NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_FT")
883	WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FT") // 2D detector data
884	WAVE/Z iErr = $("iErr_"+"FT") // 2D errors -- may not exist, especially for simulation
885	Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FT") // 2D q-values
886	WAVE inten2 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FB") // 2D detector data
887	WAVE/Z iErr2 = $("iErr_"+"FB") // 2D errors -- may not exist, especially for simulation
888	Wave qTotal2 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FB") // 2D q-values
889	nSets = 2
890	break
891
892	case "FLRTB":
893	NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_FL")
894	WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FL") // 2D detector data
895	WAVE/Z iErr = $("iErr_"+"FL") // 2D errors -- may not exist, especially for simulation
896	Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FL") // 2D q-values
897	WAVE inten2 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FR") // 2D detector data
898	WAVE/Z iErr2 = $("iErr_"+"FR") // 2D errors -- may not exist, especially for simulation
899	Wave qTotal2 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FR") // 2D q-values
900	WAVE inten3 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FT") // 2D detector data
901	WAVE/Z iErr3 = $("iErr_"+"FT") // 2D errors -- may not exist, especially for simulation
902	Wave qTotal3 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FT") // 2D q-values
903	WAVE inten4 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" + ":det_"+"FB") // 2D detector data
904	WAVE/Z iErr4 = $("iErr_"+"FB") // 2D errors -- may not exist, especially for simulation
905	Wave qTotal4 = $("root:Packages:NIST:VSANS:VCALC:" + "Front" +":qTot_"+"FB") // 2D q-values
906	nSets = 4
907	break
908
909
910	case "MLR":
911	NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_ML")
912	WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"ML") // 2D detector data
913	WAVE/Z iErr = $("iErr_"+"ML") // 2D errors -- may not exist, especially for simulation
914	Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"ML") // 2D q-values
915	WAVE inten2 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"MR") // 2D detector data
916	WAVE/Z iErr2 = $("iErr_"+"MR") // 2D errors -- may not exist, especially for simulation
917	Wave qTotal2 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"MR") // 2D q-values
918	nSets = 2
919	break
920
921	case "MTB":
922	NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_MT")
923	WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"MT") // 2D detector data
924	WAVE/Z iErr = $("iErr_"+"MT") // 2D errors -- may not exist, especially for simulation
925	Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"MT") // 2D q-values
926	WAVE inten2 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"MB") // 2D detector data
927	WAVE/Z iErr2 = $("iErr_"+"MB") // 2D errors -- may not exist, especially for simulation
928	Wave qTotal2 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"MB") // 2D q-values
929	nSets = 2
930	break
931
932	case "MLRTB":
933	NVAR delQ = $("root:Packages:NIST:VSANS:VCALC:" + "gDelQ_ML")
934	WAVE inten = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"ML") // 2D detector data
935	WAVE/Z iErr = $("iErr_"+"ML") // 2D errors -- may not exist, especially for simulation
936	Wave qTotal = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"ML") // 2D q-values
937	WAVE inten2 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"MR") // 2D detector data
938	WAVE/Z iErr2 = $("iErr_"+"MR") // 2D errors -- may not exist, especially for simulation
939	Wave qTotal2 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"MR") // 2D q-values
940	WAVE inten3 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"MT") // 2D detector data
941	WAVE/Z iErr3 = $("iErr_"+"MT") // 2D errors -- may not exist, especially for simulation
942	Wave qTotal3 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"MT") // 2D q-values
943	WAVE inten4 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" + ":det_"+"MB") // 2D detector data
944	WAVE/Z iErr4 = $("iErr_"+"MB") // 2D errors -- may not exist, especially for simulation
945	Wave qTotal4 = $("root:Packages:NIST:VSANS:VCALC:" + "Middle" +":qTot_"+"MB") // 2D q-values
946	nSets = 4
947	break
948
949	default:
950	nSets = 0 // optional default expression executed
951	Print "ERROR ---- type is not recognized "
952	endswitch
953
954	// Print "delQ = ",delQ," for ",type
955
956	if(nSets == 0)
957	return(0)
958	endif
959
960
961	//TODO: properly define the errors here - I'll have this if I do the simulation
962	if(WaveExists(iErr)==0)
963	Duplicate/O inten,iErr
964	Wave iErr=iErr
965	// iErr = 1+sqrt(inten+0.75) // can't use this -- it applies to counts, not intensity (already a count rate...)
966	iErr = sqrt(inten+0.75) // TODO -- here I'm just using some fictional value
967	endif
968
969	nq = 600
970
971	// note that the back panel of 320x320 (1mm res) results in 447 data points!
972	// - so I upped nq to 600
973
974	//****TODO**** -- where to put the averaged data -- right now, folderStr is forced to ""
975	// SetDataFolder $("root:"+folderStr) //should already be here, but make sure...
976	Make/O/D/N=(nq) $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"iBin_qxqy"+"_"+type)
977	Make/O/D/N=(nq) $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"qBin_qxqy"+"_"+type)
978	Make/O/D/N=(nq) $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"nBin_qxqy"+"_"+type)
979	Make/O/D/N=(nq) $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"iBin2_qxqy"+"_"+type)
980	Make/O/D/N=(nq) $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"eBin_qxqy"+"_"+type)
981	Make/O/D/N=(nq) $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"eBin2D_qxqy"+"_"+type)
982
983	Wave iBin_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"iBin_qxqy_"+type)
984	Wave qBin_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"qBin_qxqy"+"_"+type)
985	Wave nBin_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"nBin_qxqy"+"_"+type)
986	Wave iBin2_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"iBin2_qxqy"+"_"+type)
987	Wave eBin_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"eBin_qxqy"+"_"+type)
988	Wave eBin2D_qxqy = $("root:Packages:NIST:VSANS:VCALC:"+folderStr+"eBin2D_qxqy"+"_"+type)
989
990
991	// delQ = abs(sqrt(qx[2]^2+qy[2]^2+qz[2]^2) - sqrt(qx[1]^2+qy[1]^2+qz[1]^2)) //use bins of 1 pixel width
992	// TODO: not sure if I want to so dQ in x or y direction...
993	// the short dimension is the 8mm tubes, use this direction as dQ?
994	// but don't use the corner of the detector, since dQ will be very different on T/B or L/R due to the location of [0,0]
995	// WRT the beam center. use qx or qy directly. Still not happy with this way...
996
997
998	qBin_qxqy[] = p*delQ
999	SetScale/P x,0,delQ,"",qBin_qxqy //allows easy binning
1000
1001	iBin_qxqy = 0
1002	iBin2_qxqy = 0
1003	eBin_qxqy = 0
1004	eBin2D_qxqy = 0
1005	nBin_qxqy = 0 //number of intensities added to each bin
1006
1007	// now there are situations of:
1008	// 1 panel
1009	// 2 panels
1010	// 4 panels
1011	//
1012	// this needs to be a double loop now...
1013
1014	// use set 1 (no number) only
1015	if(nSets >= 1)
1016	xDim=DimSize(inten,0)
1017	yDim=DimSize(inten,1)
1018
1019	for(ii=0;ii<xDim;ii+=1)
1020	for(jj=0;jj<yDim;jj+=1)
1021	//qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1022	qVal = qTotal[ii][jj]
1023	binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1024	val = inten[ii][jj]
1025	if (numType(val)==0) //count only the good points, ignore Nan or Inf
1026	iBin_qxqy[binIndex] += val
1027	iBin2_qxqy[binIndex] += val*val
1028	eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
1029	nBin_qxqy[binIndex] += 1
1030	endif
1031	endfor
1032	endfor
1033
1034	endif
1035
1036	// add in set 2 (set 1 already done)
1037	if(nSets >= 2)
1038	xDim=DimSize(inten2,0)
1039	yDim=DimSize(inten2,1)
1040
1041	for(ii=0;ii<xDim;ii+=1)
1042	for(jj=0;jj<yDim;jj+=1)
1043	//qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1044	qVal = qTotal2[ii][jj]
1045	binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1046	val = inten2[ii][jj]
1047	if (numType(val)==0) //count only the good points, ignore Nan or Inf
1048	iBin_qxqy[binIndex] += val
1049	iBin2_qxqy[binIndex] += val*val
1050	eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
1051	nBin_qxqy[binIndex] += 1
1052	endif
1053	endfor
1054	endfor
1055
1056	endif
1057
1058	// add in set 3 and 4 (set 1 and 2already done)
1059	if(nSets == 4)
1060	xDim=DimSize(inten3,0)
1061	yDim=DimSize(inten3,1)
1062
1063	for(ii=0;ii<xDim;ii+=1)
1064	for(jj=0;jj<yDim;jj+=1)
1065	//qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1066	qVal = qTotal3[ii][jj]
1067	binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1068	val = inten3[ii][jj]
1069	if (numType(val)==0) //count only the good points, ignore Nan or Inf
1070	iBin_qxqy[binIndex] += val
1071	iBin2_qxqy[binIndex] += val*val
1072	eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
1073	nBin_qxqy[binIndex] += 1
1074	endif
1075	endfor
1076	endfor
1077
1078
1079	xDim=DimSize(inten4,0)
1080	yDim=DimSize(inten4,1)
1081
1082	for(ii=0;ii<xDim;ii+=1)
1083	for(jj=0;jj<yDim;jj+=1)
1084	//qTot = sqrt(qx[ii]^2 + qy[ii]^2+ qz[ii]^2)
1085	qVal = qTotal4[ii][jj]
1086	binIndex = trunc(x2pnt(qBin_qxqy, qVal))
1087	val = inten4[ii][jj]
1088	if (numType(val)==0) //count only the good points, ignore Nan or Inf
1089	iBin_qxqy[binIndex] += val
1090	iBin2_qxqy[binIndex] += val*val
1091	eBin2D_qxqy[binIndex] += iErr[ii][jj]*iErr[ii][jj]
1092	nBin_qxqy[binIndex] += 1
1093	endif
1094	endfor
1095	endfor
1096
1097	endif
1098
1099
1100	// after looping through all of the data on the panels, calculate errors on I(q),
1101	// just like in CircSectAve.ipf
1102	for(ii=0;ii<nq;ii+=1)
1103	if(nBin_qxqy[ii] == 0)
1104	//no pixels in annuli, data unknown
1105	iBin_qxqy[ii] = 0
1106	eBin_qxqy[ii] = 1
1107	eBin2D_qxqy[ii] = NaN
1108	else
1109	if(nBin_qxqy[ii] <= 1)
1110	//need more than one pixel to determine error
1111	iBin_qxqy[ii] /= nBin_qxqy[ii]
1112	eBin_qxqy[ii] = 1
1113	eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1114	else
1115	//assume that the intensity in each pixel in annuli is normally distributed about mean...
1116	iBin_qxqy[ii] /= nBin_qxqy[ii]
1117	avesq = iBin_qxqy[ii]^2
1118	aveisq = iBin2_qxqy[ii]/nBin_qxqy[ii]
1119	var = aveisq-avesq
1120	if(var<=0)
1121	eBin_qxqy[ii] = 1e-6
1122	else
1123	eBin_qxqy[ii] = sqrt(var/(nBin_qxqy[ii] - 1))
1124	endif
1125	// and calculate as it is propagated pixel-by-pixel
1126	eBin2D_qxqy[ii] /= (nBin_qxqy[ii])^2
1127	endif
1128	endif
1129	endfor
1130
1131	eBin2D_qxqy = sqrt(eBin2D_qxqy) // as equation (3) of John's memo
1132
1133	// find the last non-zero point, working backwards
1134	val=nq
1135	do
1136	val -= 1
1137	while((nBin_qxqy[val] == 0) && val > 0)
1138
1139	// print val, nBin_qxqy[val]
1140	DeletePoints val, nq-val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1141
1142	if(val == 0)
1143	// all the points were deleted
1144	return(0)
1145	endif
1146
1147
1148	// since the beam center is not always on the detector, many of the low Q bins will have zero pixels
1149	// find the first non-zero point, working forwards
1150	val = -1
1151	do
1152	val += 1
1153	while(nBin_qxqy[val] == 0)
1154	DeletePoints 0, val, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1155
1156	// ?? there still may be a point in the q-range that gets zero pixel contribution - so search this out and get rid of it
1157	val = numpnts(nBin_qxqy)-1
1158	do
1159	if(nBin_qxqy[val] == 0)
1160	DeletePoints val, 1, iBin_qxqy,qBin_qxqy,nBin_qxqy,iBin2_qxqy,eBin_qxqy,eBin2D_qxqy
1161	endif
1162	val -= 1
1163	while(val>0)
1164
1165	SetDataFolder root:
1166
1167	return(0)
1168	End
1169
1170	////////////to plot the (4) 2D panels and to plot the I(Q) data on the same plot
1171	//
1172	// ** but now I need to check and see if these waves exist before trying to append them
1173	// since the panels may bave been combined when binned - rather than all separate.
1174	//
1175	// TODO
1176	// -- so maybe I want to clear the traces from the graph?
1177	// -- set a flag on the panel to know how the binning is applied?
1178	//
1179	Window Front_IQ_Graph() : Graph
1180
1181	Variable binType
1182
1183	ControlInfo/W=VCALC popup_b
1184	binType = V_Value // V_value counts menu items from 1, so 1=1, 2=2, 3=4
1185
1186	SetDataFolder root:Packages:NIST:VSANS:VCALC
1187
1188	if(binType==1)
1189	ClearIQIfDisplayed("FLRTB")
1190	ClearIQIfDisplayed("FLR")
1191	ClearIQIfDisplayed("FTB")
1192
1193	SetDataFolder root:Packages:NIST:VSANS:VCALC
1194	CheckDisplayed/W=VCALC#Panels_IQ iBin_qxqy_FL
1195
1196	if(V_flag==0)
1197	AppendtoGraph/W=VCALC#Panels_IQ iBin_qxqy_FL vs qBin_qxqy_FL
1198	AppendToGraph/W=VCALC#Panels_IQ iBin_qxqy_FR vs qBin_qxqy_FR
1199	AppendToGraph/W=VCALC#Panels_IQ iBin_qxqy_FT vs qBin_qxqy_FT
1200	AppendToGraph/W=VCALC#Panels_IQ iBin_qxqy_FB vs qBin_qxqy_FB
1201	ModifyGraph/W=VCALC#Panels_IQ mode=4
1202	ModifyGraph/W=VCALC#Panels_IQ marker=19
1203	ModifyGraph/W=VCALC#Panels_IQ rgb(iBin_qxqy_FL)=(39321,26208,1),rgb(iBin_qxqy_FB)=(2,39321,1),rgb(iBin_qxqy_FR)=(39321,26208,1),rgb(iBin_qxqy_FT)=(2,39321,1)
1204	ModifyGraph/W=VCALC#Panels_IQ msize=2
1205	ModifyGraph/W=VCALC#Panels_IQ muloffset(iBin_qxqy_FL)={0,4},muloffset(iBin_qxqy_FB)={0,2},muloffset(iBin_qxqy_FR)={0,8}
1206	ModifyGraph/W=VCALC#Panels_IQ grid=1
1207	ModifyGraph/W=VCALC#Panels_IQ log=1
1208	ModifyGraph/W=VCALC#Panels_IQ mirror=2
1209	Label/W=VCALC#Panels_IQ left "Intensity (1/cm)"
1210	Label/W=VCALC#Panels_IQ bottom "Q (1/A)"
1211	endif
1212
1213	endif
1214
1215	if(binType==2)
1216	ClearIQIfDisplayed("FLRTB")
1217	ClearIQIfDisplayed("FT")
1218	ClearIQIfDisplayed("FL")
1219	ClearIQIfDisplayed("FR")
1220	ClearIQIfDisplayed("FB")
1221
1222
1223	SetDataFolder root:Packages:NIST:VSANS:VCALC
1224	CheckDisplayed/W=VCALC#Panels_IQ iBin_qxqy_FLR
1225
1226	if(V_flag==0)
1227	AppendtoGraph/W=VCALC#Panels_IQ iBin_qxqy_FLR vs qBin_qxqy_FLR
1228	AppendToGraph/W=VCALC#Panels_IQ iBin_qxqy_FTB vs qBin_qxqy_FTB
1229	ModifyGraph/W=VCALC#Panels_IQ mode=4
1230	ModifyGraph/W=VCALC#Panels_IQ marker=19
1231	ModifyGraph/W=VCALC#Panels_IQ rgb(iBin_qxqy_FLR)=(39321,26208,1),rgb(iBin_qxqy_FTB)=(2,39321,1)
1232	ModifyGraph/W=VCALC#Panels_IQ msize=2
1233	ModifyGraph/W=VCALC#Panels_IQ muloffset(iBin_qxqy_FLR)={0,2}
1234	ModifyGraph/W=VCALC#Panels_IQ grid=1
1235	ModifyGraph/W=VCALC#Panels_IQ log=1
1236	ModifyGraph/W=VCALC#Panels_IQ mirror=2
1237	Label/W=VCALC#Panels_IQ left "Intensity (1/cm)"
1238	Label/W=VCALC#Panels_IQ bottom "Q (1/A)"
1239	endif
1240
1241	endif
1242
1243	if(binType==3)
1244	ClearIQIfDisplayed("FLR")
1245	ClearIQIfDisplayed("FTB")
1246	ClearIQIfDisplayed("FT")
1247	ClearIQIfDisplayed("FL")
1248	ClearIQIfDisplayed("FR")
1249	ClearIQIfDisplayed("FB")
1250
1251	SetDataFolder root:Packages:NIST:VSANS:VCALC
1252	CheckDisplayed/W=VCALC#Panels_IQ iBin_qxqy_FLRTB
1253
1254	if(V_flag==0)
1255	AppendtoGraph/W=VCALC#Panels_IQ iBin_qxqy_FLRTB vs qBin_qxqy_FLRTB
1256	ModifyGraph/W=VCALC#Panels_IQ mode=4
1257	ModifyGraph/W=VCALC#Panels_IQ marker=19
1258	ModifyGraph/W=VCALC#Panels_IQ rgb(iBin_qxqy_FLRTB)=(39321,26208,1)
1259	ModifyGraph/W=VCALC#Panels_IQ msize=2
1260	ModifyGraph/W=VCALC#Panels_IQ grid=1
1261	ModifyGraph/W=VCALC#Panels_IQ log=1
1262	ModifyGraph/W=VCALC#Panels_IQ mirror=2
1263	Label/W=VCALC#Panels_IQ left "Intensity (1/cm)"
1264	Label/W=VCALC#Panels_IQ bottom "Q (1/A)"
1265	endif
1266
1267	endif
1268
1269
1270	if(binType==4) //slit mode
1271	ClearIQIfDisplayed("FLRTB")
1272	ClearIQIfDisplayed("FLR")
1273	ClearIQIfDisplayed("FTB")
1274
1275	SetDataFolder root:Packages:NIST:VSANS:VCALC
1276	CheckDisplayed/W=VCALC#Panels_IQ iBin_qxqy_FL
1277
1278	if(V_flag==0)
1279	AppendtoGraph/W=VCALC#Panels_IQ iBin_qxqy_FL vs qBin_qxqy_FL
1280	AppendToGraph/W=VCALC#Panels_IQ iBin_qxqy_FR vs qBin_qxqy_FR
1281	AppendToGraph/W=VCALC#Panels_IQ iBin_qxqy_FT vs qBin_qxqy_FT
1282	AppendToGraph/W=VCALC#Panels_IQ iBin_qxqy_FB vs qBin_qxqy_FB
1283	ModifyGraph/W=VCALC#Panels_IQ mode=4
1284	ModifyGraph/W=VCALC#Panels_IQ marker=19
1285	ModifyGraph/W=VCALC#Panels_IQ rgb(iBin_qxqy_FL)=(39321,26208,1),rgb(iBin_qxqy_FB)=(2,39321,1),rgb(iBin_qxqy_FR)=(39321,26208,1),rgb(iBin_qxqy_FT)=(2,39321,1)
1286	ModifyGraph/W=VCALC#Panels_IQ msize=2
1287	ModifyGraph/W=VCALC#Panels_IQ muloffset(iBin_qxqy_FL)={0,4},muloffset(iBin_qxqy_FB)={0,2},muloffset(iBin_qxqy_FR)={0,8}
1288	ModifyGraph/W=VCALC#Panels_IQ grid=1
1289	ModifyGraph/W=VCALC#Panels_IQ log=1
1290	ModifyGraph/W=VCALC#Panels_IQ mirror=2
1291	Label/W=VCALC#Panels_IQ left "Intensity (1/cm)"
1292	Label/W=VCALC#Panels_IQ bottom "Q (1/A)"
1293	endif
1294
1295	endif
1296
1297
1298	SetDataFolder root:
1299
1300	EndMacro
1301
1302	Function ClearIQIfDisplayed(type)
1303	String type
1304
1305	SetDataFolder root:Packages:NIST:VSANS:VCALC
1306	CheckDisplayed/W=VCALC#Panels_IQ $("iBin_qxqy_"+type)
1307	if(V_flag==1)
1308	RemoveFromGraph/W=VCALC#Panels_IQ $("iBin_qxqy_"+type)
1309	endif
1310	SetDataFolder root:
1311
1312	return(0)
1313	end
1314
1315	Window Table_of_QBins() : Table
1316	PauseUpdate; Silent 1 // building window...
1317	String fldrSav0= GetDataFolder(1)
1318	SetDataFolder root:Packages:NIST:VSANS:VCALC:
1319	Edit/W=(5,44,771,898) qBin_qxqy_FL,qBin_qxqy_FR,qBin_qxqy_FT,qBin_qxqy_FB,qBin_qxqy_FLR
1320	AppendToTable qBin_qxqy_FTB,qBin_qxqy_FLRTB
1321	ModifyTable format(Point)=1,width(qBin_qxqy_FLR)=136,width(qBin_qxqy_FLRTB)=120
1322	SetDataFolder fldrSav0
1323	EndMacro
1324
1325
1326
1327
1328
1329
1330
1331	//////NOTE///
1332	//// some chunks of the code here have been trimmed out for fun
1333	//Function V_CircularAverageTo1D(type)
1334	// String type
1335	//
1336	//////// get information about the detector (in type folder) that is needed for reduction
1337	//// pixel dimensions
1338	//// beam center
1339	//// distances
1340	//// wavelength
1341	//// wavelength spread
1342	////
1343	//
1344	// String destPath = "root:"
1345	//
1346	//// NVAR pixelsX = root:myGlobals:gNPixelsX
1347	//// NVAR pixelsY = root:myGlobals:gNPixelsY
1348	//
1349	//// pixelsX = 48
1350	//// pixelsY = 250
1351	//
1352	// // this is for non-linear corrections not applicable?
1353	//// xcenter = pixelsX/2 + 0.5 // == 64.5 for 128x128 Ordela
1354	//// ycenter = pixelsY/2 + 0.5 // == 64.5 for 128x128 Ordela
1355	//
1356	// // beam center, in pixels
1357	// x0 = 60 //reals[16]
1358	// y0 = 125 //reals[17]
1359	//// //detector calibration constants
1360	//// sx = reals[10] //mm/pixel (x)
1361	//// sx3 = reals[11] //nonlinear coeff
1362	//// sy = reals[13] //mm/pixel (y)
1363	//// sy3 = reals[14] //nonlinear coeff
1364	//
1365	//// dtsize = 10*reals[20] //det size in mm
1366	// dtdist = 1000 //1000*reals[18] // det distance in mm
1367	//
1368	//
1369	/////// decide how the binning width is to be determined
1370	//
1371	//// NVAR binWidth=root:Packages:NIST:gBinWidth
1372	//
1373	// dr = 1 //binWidth // *********annulus width set by user, default is one*********
1374	// ddr = 4 //dr*sx //step size, in mm (this value should be passed to the resolution calculation, not dr 18NOV03)
1375	//
1376	// Variable rcentr,large_num,small_num,dtdis2,nq,xoffst,dxbm,dybm,ii
1377	// Variable phi_rad,dphi_rad,phi_x,phi_y
1378	// Variable forward,mirror
1379	//
1380	////// do I need to pick sides?? probably, for consistency, but confusing nomenclature now
1381	//// String side = StringByKey("SIDE",keyListStr,"=",";")
1382	//
1383	///////// keep the sector calculations, I'll want to do this...
1384	//////
1385	//// if(!isCircular) //must be sector avg (rectangular not sent to this function)
1386	//// //convert from degrees to radians
1387	//// phi_rad = (Pi/180)*NumberByKey("PHI",keyListStr,"=",";")
1388	//// dphi_rad = (Pi/180)*NumberByKey("DPHI",keyListStr,"=",";")
1389	//// //create cartesian values for unit vector in phi direction
1390	//// phi_x = cos(phi_rad)
1391	//// phi_y = sin(phi_rad)
1392	//// Endif
1393	//
1394	// /// data wave is data in the current folder which was set at the top of the function
1395	//// WAVE data=$(destPath + ":data")
1396	// Make/O/N=(pixelsX,pixelsY) $(destPath + ":data")
1397	// Wave data = $(destPath + ":data")
1398	// data = 1
1399	//
1400	// //Check for the existence of the mask, if not, make one (local to this folder) that is null
1401	//
1402	// if(WaveExists($"root:Packages:NIST:MSK:data") == 0)
1403	// Print "There is no mask file loaded (WaveExists)- the data is not masked"
1404	// Make/O/N=(pixelsX,pixelsY) $(destPath + ":mask")
1405	// Wave mask = $(destPath + ":mask")
1406	// mask = 0
1407	// else
1408	// Wave mask=$"root:Packages:NIST:MSK:data"
1409	// Endif
1410	//
1411	// //
1412	// //pixels within rcentr of beam center are broken into 9 parts (units of mm)
1413	// rcentr = 100 //original
1414	//// rcentr = 0
1415	// // values for error if unable to estimate value
1416	// //large_num = 1e10
1417	// large_num = 1 //1e10 value (typically sig of last data point) plots poorly, arb set to 1
1418	// small_num = 1e-10
1419	//
1420	// // output wave are expected to exist (?) initialized to zero, what length?
1421	// // 200 points on VAX --- use 300 here, or more if SAXS data is used with 1024x1024 detector (1000 pts seems good)
1422	// Variable defWavePts=500
1423	// Make/O/N=(defWavePts) $(destPath + ":qval"),$(destPath + ":aveint")
1424	// Make/O/N=(defWavePts) $(destPath + ":ncells"),$(destPath + ":dsq"),$(destPath + ":sigave")
1425	// Make/O/N=(defWavePts) $(destPath + ":SigmaQ"),$(destPath + ":fSubS"),$(destPath + ":QBar")
1426	//
1427	// WAVE qval = $(destPath + ":qval")
1428	// WAVE aveint = $(destPath + ":aveint")
1429	// WAVE ncells = $(destPath + ":ncells")
1430	// WAVE dsq = $(destPath + ":dsq")
1431	// WAVE sigave = $(destPath + ":sigave")
1432	// WAVE qbar = $(destPath + ":QBar")
1433	// WAVE sigmaq = $(destPath + ":SigmaQ")
1434	// WAVE fsubs = $(destPath + ":fSubS")
1435	//
1436	// qval = 0
1437	// aveint = 0
1438	// ncells = 0
1439	// dsq = 0
1440	// sigave = 0
1441	// qbar = 0
1442	// sigmaq = 0
1443	// fsubs = 0
1444	//
1445	// dtdis2 = dtdist^2
1446	// nq = 1
1447	// xoffst=0
1448	// //distance of beam center from detector center
1449	////// // the linearity corrections for the 2D Ordela detectors are applied from the center of the detector,
1450	// // so figure where that is, relative to the beam center.
1451	//// dxbm = V_FX(x0,sx3,xcenter,sx)
1452	//// dybm = V_FY(y0,sy3,ycenter,sy)
1453	//
1454	// //BEGIN AVERAGE **********
1455	// Variable xi,dxi,dx,jj,data_pixel,yj,dyj,dy,mask_val=0.1
1456	// Variable dr2,nd,fd,nd2,ll,kk,dxx,dyy,ir,dphi_p
1457	//
1458	// // IGOR arrays are indexed from [0][0], FORTAN from (1,1) (and the detector too)
1459	// // loop index corresponds to FORTRAN (old code)
1460	// // and the IGOR array indices must be adjusted (-1) to the correct address
1461	//
1462	// /////
1463	// //// need to add in here a step that calculates the q-values, and bins based on q-value
1464	// //// since the 4 panels are not in the same plane (but relatively close)
1465	// //
1466	// // if we're always using pairs in the same plane, then binning in r is OK
1467	// ////
1468	//
1469	// ii=1
1470	// do
1471	// xi = ii
1472	// dxi = V_FX(xi,sx3,xcenter,sx)
1473	// dx = dxi-dxbm //dx and dy are in mm
1474	//
1475	// jj = 1
1476	// do
1477	// data_pixel = data[ii-1][jj-1] //assign to local variable
1478	// yj = jj
1479	// dyj = V_FY(yj,sy3,ycenter,sy)
1480	// dy = dyj - dybm
1481	// if(!(mask[ii-1][jj-1])) //masked pixels = 1, skip if masked (this way works...)
1482	// dr2 = (dx^2 + dy^2)^(0.5) //distance from beam center NOTE dr2 used here - dr used above
1483	// if(dr2>rcentr) //keep pixel whole
1484	// nd = 1
1485	// fd = 1
1486	// else //break pixel into 9 equal parts
1487	// nd = 3
1488	// fd = 2
1489	// endif
1490	// nd2 = nd^2
1491	// ll = 1 //"el-el" loop index
1492	// do
1493	// dxx = dx + (ll - fd)*sx/3
1494	// kk = 1
1495	// do
1496	// dyy = dy + (kk - fd)*sy/3
1497	// if(isCircular)
1498	// //circular average, use all pixels
1499	// //(increment)
1500	// nq = V_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1501	//// else
1502	//// //a sector average - determine azimuthal angle
1503	//// dphi_p = V_dphi_pixel(dxx,dyy,phi_x,phi_y)
1504	//// if(dphi_p < dphi_rad)
1505	//// forward = 1 //within forward sector
1506	//// else
1507	//// forward = 0
1508	//// Endif
1509	//// if((Pi - dphi_p) < dphi_rad)
1510	//// mirror = 1 //within mirror sector
1511	//// else
1512	//// mirror = 0
1513	//// Endif
1514	//// //check if pixel lies within allowed sector(s)
1515	//// if(cmpstr(side,"both")==0) //both sectors
1516	//// if ( mirror \|\| forward)
1517	//// //increment
1518	//// nq = V_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1519	//// Endif
1520	//// else
1521	//// if(cmpstr(side,"right")==0) //forward sector only
1522	//// if(forward)
1523	//// //increment
1524	//// nq = V_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1525	//// Endif
1526	//// else //mirror sector only
1527	//// if(mirror)
1528	//// //increment
1529	//// nq = V_IncrementPixel(data_pixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1530	//// Endif
1531	//// Endif
1532	//// Endif //allowable sectors
1533	// Endif //circular or sector check
1534	// kk+=1
1535	// while(kk<=nd)
1536	// ll += 1
1537	// while(ll<=nd)
1538	// Endif //masked pixel check
1539	// jj += 1
1540	// while (jj<=pixelsY)
1541	// ii += 1
1542	// while(ii<=pixelsX) //end of the averaging
1543	//
1544	// //compute q-values and errors
1545	// Variable ntotal,rr,theta,avesq,aveisq,var
1546	//
1547	// lambda = reals[26]
1548	// ntotal = 0
1549	// kk = 1
1550	// do
1551	// rr = (2kk-1)ddr/2
1552	// theta = 0.5*atan(rr/dtdist)
1553	// qval[kk-1] = (4Pi/lambda)sin(theta)
1554	// if(ncells[kk-1] == 0)
1555	// //no pixels in annuli, data unknown
1556	// aveint[kk-1] = 0
1557	// sigave[kk-1] = large_num
1558	// else
1559	// if(ncells[kk-1] <= 1)
1560	// //need more than one pixel to determine error
1561	// aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
1562	// sigave[kk-1] = large_num
1563	// else
1564	// //assume that the intensity in each pixel in annuli is normally
1565	// // distributed about mean...
1566	// aveint[kk-1] = aveint[kk-1]/ncells[kk-1]
1567	// avesq = aveint[kk-1]^2
1568	// aveisq = dsq[kk-1]/ncells[kk-1]
1569	// var = aveisq-avesq
1570	// if(var<=0)
1571	// sigave[kk-1] = small_num
1572	// else
1573	// sigave[kk-1] = sqrt(var/(ncells[kk-1] - 1))
1574	// endif
1575	// endif
1576	// endif
1577	// ntotal += ncells[kk-1]
1578	// kk+=1
1579	// while(kk<=nq)
1580	//
1581	// //Print "NQ = ",nq
1582	// // data waves were defined as 300 points (=defWavePts), but now have less than that (nq) points
1583	// // use DeletePoints to remove junk from end of waves
1584	// //WaveStats would be a more foolproof implementation, to get the # points in the wave
1585	// Variable startElement,numElements
1586	// startElement = nq
1587	// numElements = defWavePts - startElement
1588	// DeletePoints startElement,numElements, qval,aveint,ncells,dsq,sigave
1589	//
1590	// //////////////end of VAX sector_ave()
1591	//
1592	// //angle dependent transmission correction
1593	// Variable uval,arg,cos_th
1594	// lambda = reals[26]
1595	// trans = reals[4]
1596	//
1597	////
1598	//// The transmission correction is now done at the ADD step, in DetCorr()
1599	////
1600	//// ////this section is the trans_correct() VAX routine
1601	//// if(trans<0.1)
1602	//// Print "*transmission is less than 0.1* and is a significant correction"
1603	//// endif
1604	//// if(trans==0)
1605	//// Print "*transmission is ZERO* and has been reset to 1.0 for the averaging calculation"
1606	//// trans = 1
1607	//// endif
1608	//// //optical thickness
1609	//// uval = -ln(trans) //use natural logarithm
1610	//// //apply correction to aveint[]
1611	//// //index from zero here, since only working with IGOR waves
1612	//// ii=0
1613	//// do
1614	//// theta = 2asin(lambdaqval[ii]/(4*pi))
1615	//// cos_th = cos(theta)
1616	//// arg = (1-cos_th)/cos_th
1617	//// if((uval<0.01) \|\| (cos_th>0.99)) //OR
1618	//// //small arg, approx correction
1619	//// aveint[ii] /= 1-0.5uvalarg
1620	//// else
1621	//// //large arg, exact correction
1622	//// aveint[ii] /= (1-exp(-uvalarg))/(uvalarg)
1623	//// endif
1624	//// ii+=1
1625	//// while(ii<nq)
1626	//// //end of transmission/pathlength correction
1627	//
1628	//// ***************************************************************
1629	////
1630	//// Do the extra 3 columns of resolution calculations starting here.
1631	////
1632	//// ***************************************************************
1633	//
1634	// Variable L2 = reals[18]
1635	// Variable BS = reals[21]
1636	// Variable S1 = reals[23]
1637	// Variable S2 = reals[24]
1638	// Variable L1 = reals[25]
1639	// lambda = reals[26]
1640	// Variable lambdaWidth = reals[27]
1641	// String detStr=textRead[9]
1642	//
1643	// Variable usingLenses = reals[28] //new 2007
1644	//
1645	// //Two parameters DDET and APOFF are instrument dependent. Determine
1646	// //these from the instrument name in the header.
1647	// //From conversation with JB on 01.06.99 these are the current
1648	// //good values
1649	//
1650	// Variable DDet
1651	// NVAR apOff = root:myGlobals:apOff //in cm
1652	//
1653	//// DDet = DetectorPixelResolution(fileStr,detStr) //needs detector type and beamline
1654	// //note that reading the detector pixel size from the header ASSUMES SQUARE PIXELS! - Jan2008
1655	// DDet = reals[10]/10 // header value (X) is in mm, want cm here
1656	//
1657	//
1658	// //Width of annulus used for the average is gotten from the
1659	// //input dialog before. This also must be passed to the resolution
1660	// //calculator. Currently the default is dr=1 so just keeping that.
1661	//
1662	// //Go from 0 to nq doing the calc for all three values at
1663	// //every Q value
1664	//
1665	// ii=0
1666	//
1667	// Variable ret1,ret2,ret3
1668	// do
1669	// // commented out for compiler
1670	//// getResolution(qval[ii],lambda,lambdaWidth,DDet,apOff,S1,S2,L1,L2,BS,ddr,usingLenses,ret1,ret2,ret3)
1671	// sigmaq[ii] = ret1
1672	// qbar[ii] = ret2
1673	// fsubs[ii] = ret3
1674	// ii+=1
1675	// while(ii<nq)
1676	// DeletePoints startElement,numElements, sigmaq, qbar, fsubs
1677	//
1678	//// End of resolution calculations
1679	//// ***************************************************************
1680	//
1681	// //Plot the data in the Plot_1d window
1682	//// Avg_1D_Graph(aveint,qval,sigave) //commented out for compiler
1683	//
1684	// //get rid of the default mask, if one was created (it is in the current folder)
1685	// //don't just kill "mask" since it might be pointing to the one in the MSK folder
1686	// Killwaves/Z $(destPath+":mask")
1687	//
1688	// //return to root folder (redundant)
1689	// SetDataFolder root:
1690	//
1691	// Return 0
1692	//End
1693
1694	////returns nq, new number of q-values
1695	////arrays aveint,dsq,ncells are also changed by this function
1696	////
1697	//Function V_IncrementPixel(dataPixel,ddr,dxx,dyy,aveint,dsq,ncells,nq,nd2)
1698	// Variable dataPixel,ddr,dxx,dyy
1699	// Wave aveint,dsq,ncells
1700	// Variable nq,nd2
1701	//
1702	// Variable ir
1703	//
1704	// ir = trunc(sqrt(dxxdxx+dyydyy)/ddr)+1
1705	// if (ir>nq)
1706	// nq = ir //resets maximum number of q-values
1707	// endif
1708	// aveint[ir-1] += dataPixel/nd2 //ir-1 must be used, since ir is physical
1709	// dsq[ir-1] += dataPixel*dataPixel/nd2
1710	// ncells[ir-1] += 1/nd2
1711	//
1712	// Return nq
1713	//End
1714	//
1715	////function determines azimuthal angle dphi that a vector connecting
1716	////center of detector to pixel makes with respect to vector
1717	////at chosen azimuthal angle phi -> [cos(phi),sin(phi)] = [phi_x,phi_y]
1718	////dphi is always positive, varying from 0 to Pi
1719	////
1720	//Function V_dphi_pixel(dxx,dyy,phi_x,phi_y)
1721	// Variable dxx,dyy,phi_x,phi_y
1722	//
1723	// Variable val,rr,dot_prod
1724	//
1725	// rr = sqrt(dxx^2 + dyy^2)
1726	// dot_prod = (dxxphi_x + dyyphi_y)/rr
1727	// //? correct for roundoff error? - is this necessary in IGOR, w/ double precision?
1728	// if(dot_prod > 1)
1729	// dot_prod =1
1730	// Endif
1731	// if(dot_prod < -1)
1732	// dot_prod = -1
1733	// Endif
1734	//
1735	// val = acos(dot_prod)
1736	//
1737	// return val
1738	//
1739	//End
1740	//
1741	////calculates the x distance from the center of the detector, w/nonlinear corrections
1742	////
1743	//Function V_FX(xx,sx3,xcenter,sx)
1744	// Variable xx,sx3,xcenter,sx
1745	//
1746	// Variable retval
1747	//
1748	// retval = sx3tan((xx-xcenter)sx/sx3)
1749	// Return retval
1750	//End
1751	//
1752	////calculates the y distance from the center of the detector, w/nonlinear corrections
1753	////
1754	//Function V_FY(yy,sy3,ycenter,sy)
1755	// Variable yy,sy3,ycenter,sy
1756	//
1757	// Variable retval
1758	//
1759	// retval = sy3tan((yy-ycenter)sy/sy3)
1760	// Return retval
1761	//End

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

Download in other formats:

Original Format