#pragma rtGlobals=1		// Use modern global access method.
#pragma version = 2.5
// with version 2.4, I added the offset to bx in the calculation
// This is the new version that deals with multiple detecgtors, and has the beam along the z-axis
// The traditional beam line coordinates
#include  "Elements", version>=1.1
#include  "Neutron", version>=1.0
#include "vector-math"

// for the new motor used in September, depth = zi - surface - translator


Menu "SCD"
	"行行 General 行行"
	"tell about spot at cursor A",InfoAboutSpotOnDetector(NaN,NaN,0.13)
	"Channel to Wavelength",chan2lambda(NaN)
	"Wavelength to Channel",lambda2chan(Nan)
	"-"
	"行行 Reciprocal Lattice 行行"
	"find hkl of cursor using reip",print hkl_fromDetectorSpot_recip(NaN,NaN,$"")
	"angle between G vectors at cursors", print angleBetweenGofCursors()
	"angle between two vectors", angleBetweenHKLs(NaN,NaN,NaN,NaN,NaN,NaN)
	"detector xy from hkl and recip",printDetectorXY_of_hkl($"",NaN,NaN,NaN)
	"add predicted hkl spots",addPredictedHKL($"",NaN,NaN)
	"remove predicted hkl spots", SetDrawLayer /K UserFront		// clear the drawing layer
	"wavelength of hkl from recip",hkl2lambda($"",NaN,NaN,NaN)
	"hkl pointing at source",hkl_BeamAxis($"")
	"calc recip lattice from 2 cursors", spots2recip(NaN,NaN,NaN,NaN)
	"calc recip lattice from 2 cursors xx", xxspots2recip("",NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN)
	"-"
	"行行 Displaying 行行"
	"Show a Slice Step Graph", ShowSliceGraph("")
	"Run the Slice as a Movie", RunAsMovie()
	"Apply-Remove gamma Correction",ApplyGamma(NaN)
	"Make Movie of Top Slice Graph", makeMovieOfTopSliceGraph(NaN)
	"show SCD geometry",Panel_SCD_geometry()
	"-"
	"行行 Sorting and Reading In 行行"
	"Sort Out 1 Depth", SortDataOneDepth(NaN,NaN,"","",0)
	"Sum Up Over Wavelength", sumUpWavelength("","",NaN)
	"Make Sliced Sums at Multiple Depths", MakeSlicedSumsAtManyDepths("","",NaN,NaN,NaN,NaN)
	"Set One Geometry Parameters",DetectorGeometryWindow()
	"ReadRangeSCD", ReadRangeSCD("","",nan,"","")
	"Set Translator Distances", setTranslatorDistances("","",NaN,NaN)
	"-"
	"Initialize SCD package",SCDinitPackage(1)
End


Proc GraphImageStyle() : GraphStyle
	PauseUpdate; Silent 1		// modifying window...
	ModifyGraph/Z gfSize=14,width={Aspect,1}, tick=2, minor=1, lowTrip=0.001,standoff=0, mirror=1
	ModifyGraph/Z axOffset(left)=-2.3,axOffset(bottom)=-1.11111
	ModifyImage $StringFromList(0,ImageNameList("",";")) ctab= {*,*,Terrain,1}
//	ModifyImage $StringFromList(0,ImageNameList("",";")) ctab= {*,*,Rainbow,0}
EndMacro


// ====================================================================
// ====================================================================
// =========================  Start of Reciprocal Lattice   =======================
// ====================================================================

Function/T hkl_fromDetectorSpot_recip(detX,detY,recip)	// computes hkl of a pixel using recip
	Variable detX,detY
	Wave/D recip
	if (!WaveExists(recip) || DimSize(recip,0)!=3 || DimSize(recip,1)!=3 || DimSize(recip,2)!=0)
		String recipLattice, mat3x3List=WaveList_NxN(3,3)
		recipLattice = SelectString(exists("recip")==1,"","recip")
		Prompt recipLattice, "3x3 matrix with reciprocal lattice",popup,mat3x3List
		DoPrompt "reciprocal lattice",recipLattice
		if (V_flag)
			return ""
		endif
		Wave/D recip = $recipLattice
	endif

	detX = numtype(detX) ? hcsr(A) : detX
	detY = numtype(detY) ? vcsr(A) : detY
	String noteStr = note(ImageNameToWaveRef("",StringFromList(0,ImageNameList("",";"))))
	Variable deta = valFromWaveNoteOrGlobal(noteStr,"deta")
	Variable detd = valFromWaveNoteOrGlobal(noteStr,"detd")
	Variable bx = valFromWaveNoteOrGlobal(noteStr,"bx")
	Variable xp, zp, yp=detY						// y is vertical
	pixelXZfromDetector(deta,detd,detX,xp,zp)		// calculate pixel coordinates xp,zp

	Variable lambda = getLambdaFromPlot("")		// get lambda from a step slice plot
	Variable k_len = 2*PI/lambda
	Make/N=3/O kf, ki, Gvec
	ki={0, 0, k_len}
//	kf={xp,yp,zp}
	kf={xp-bx,yp,zp}
	normalize(kf)
	kf *= k_len										// kf and ki need to be k_len long
	Gvec = kf - ki

	// recip x H = G  -->  recipInv x recip x H = recipInv x G
	// hkl = recipInv x Gvec, so find inverse of recip
	Make/N=(3,3)/O/D recip_temp_,recipInv
	recipInv = (p==q)
	recip_temp_ = recip
	MatrixLinearSolve /O recip_temp_ recipInv		// recipInv now contains the inverse of hklN
	KillWaves/Z recip_temp_

	MatrixMultiply recipInv, Gvec					// result is in M_product
	Wave hkl = M_product
	String out
	sprintf out, "%.0f %.0f %.0f",hkl[0],hkl[1],hkl[2]

	KillWaves/Z ki,kf,Gvec, W_IPIV, recipInv, M_product
	return out
End



Function hkl2lambda(recip,h,k,l)			// find the wavelength for an hkl using the reciprocal lattice
	Wave/D recip							// 3x3 matrix with reciprocal lattice
	Variable h,k,l

	Variable printIt=0
	if (numtype(h) || numtype(k) || numtype(l) || !WaveExists(recip))
		h = numtype(h) ? 0 : h
		k = numtype(k) ? 2 : k
		l = numtype(l) ? 2 : l
		Prompt h, "h"
		Prompt k, "k"
		Prompt l, "l"
		String recipLattice, mat3x3List=WaveList_NxN(3,3)
		recipLattice = SelectString(exists("recip")==1,"","recip")
		Prompt recipLattice, "3x3 matrix with reciprocal lattice",popup,mat3x3List
		DoPrompt "hkl",recipLattice,h,k,l
		if (V_flag)
			return NaN
		endif
		Wave/D recip = $recipLattice
		printIt = 1
	endif
	if (numtype(h) || numtype(k) || numtype(l) || !WaveExists(recip))
		return NaN
	endif

	Make/N=3/D/O hkl__, kihat={0,0,1}
	hkl__ = {h,k,l}
	MatrixMultiply recip, hkl__	// result is in M_product
	KillWaves/Z Ghat
	Rename M_product Ghat			// G vector (from hkl__)
	Redimension/N=3 Ghat			// needed because dimSize(Ghat,1)==1, not 0
	Variable G = normalize(Ghat)	// d = 2�/G
		//   cos(alpha) = -MatrixDot(Ghat,kihat)
		//	theta = �/2-alpha
		//	sin(theta) = sin(�/2-alpha) = cos(alpha),  but I know cos(alpha) from the dot product
		//	sin(theta) = -MatrixDot(Ghat,kihat)
	Variable lambda = 4*PI/G * (-MatrixDot(Ghat,kihat))	// 4� sin(theta) / G
	lambda = (lambda<0) ? NaN :  lambda
	if (printIt)
		printf " [%s] x (%d %d %d)  -->  lambda=%g\r",NameOfWave(recip),h,k,l,lambda
	endif
	KillWaves/Z Ghat, kihat, hkl__
	return lambda
End



// find where (hkl) should land on the detector from reciprocal lattice, and print to history
Function/T printDetectorXY_of_hkl(recip,h,k,l)
	Wave/D recip							// 3x3 matrix with reciprocal lattice
	Variable h,k,l

	if (numtype(h) || numtype(k) || numtype(l) || !WaveExists(recip))
		h = numtype(h) ? 0 : h
		k = numtype(k) ? 2 : k
		l = numtype(l) ? 2 : l
		Prompt h, "h"
		Prompt k, "k"
		Prompt l, "l"
		String recipLattice, mat3x3List=WaveList_NxN(3,3)
		recipLattice = SelectString(exists("recip")==1,"","recip")
		Prompt recipLattice, "3x3 matrix with reciprocal lattice",popup,mat3x3List
		DoPrompt "hkl",recipLattice,h,k,l
		if (V_flag)
			return ""
		endif
		Wave/D recip = $recipLattice
	endif
	if (numtype(h) || numtype(k) || numtype(l) || !WaveExists(recip))
		return ""
	endif

	String graphName = StringFromList(0,WinList("*",";","WIN:"))
	if (strlen(graphName)<1)
		Abort "No Window for doing determining parameters"
	endif
	Wave image = ImageNameToWaveRef(graphName,StringFromList(0,ImageNameList(graphName,";")))
	if (!WaveExists(image))
		return ""
	endif

	Variable detX, detY									// position on detector (mm)
	Make/O hkl__temp__ = {h,k,l}
	Wave hkl = hkl__temp__
	calcXYofHKL(recip,hkl,detX,detY,note(image))		// find xy on the detector
	printf "for hkl=(%d %d %d) [%s] --> detector=(%g, %g)\r",hkl[0],hkl[1],hkl[2],NameOfWave(recip),detX,detY
	KillWaves hkl__temp__
	String result
	sprintf result "detX=%g;detY=%g;",detX,detY
	return result
End



Function/T hkl_BeamAxis(recip)	// computes hkl of a Gvector that is from sample back toward the source
	Wave/D recip
	Variable printIt=0
	if (!WaveExists(recip) || DimSize(recip,0)!=3 || DimSize(recip,1)!=3 || DimSize(recip,2)!=0)
		String recipLattice, mat3x3List=WaveList_NxN(3,3)
		recipLattice = SelectString(exists("recip")==1,"","recip")
		Prompt recipLattice, "3x3 matrix with reciprocal lattice",popup,mat3x3List
		DoPrompt "reciprocal lattice",recipLattice
		if (V_flag)
			return ""
		endif
		Wave/D recip = $recipLattice
		printIt = 1
	endif
	// recip x H = G  -->  recipInv x recip x H = recipInv x G
	// hkl = recipInv x Gvec, so find inverse of recip
	Make/N=3/O Gvec={0,0,-2*PI}				// back reflection for lambda=1�
	Make/N=(3,3)/O/D recip_temp_,recipInv
	recipInv = (p==q)
	recip_temp_ = recip
	MatrixLinearSolve /O recip_temp_ recipInv		// recipInv now contains the inverse of hklN
	KillWaves/Z recip_temp_

	MatrixMultiply recipInv, Gvec					// result is in M_product
	Wave hkl = M_product
	String out
	sprintf out, "%.2f, %.2f, %.2f",hkl[0],hkl[1],hkl[2]
	KillWaves/Z Gvec, W_IPIV, recipInv, M_product
	if (printIt)
		printf "[%s]x[hkl=(%s)] = G,  points at the source\r",NameOfWave(recip),out
	endif
	return out
End



Function addPredictedHKL(recip, hklMax,lamMin)
	Wave/D recip
	Variable hklMax
	Variable lamMin				// usually 0.42 �
	if (!WaveExists(recip) ||  numtype(hklMax) || hklMax<1 || lamMin<0 || numtype(lamMin))
		lamMin = (lamMin<0 || numtype(lamMin)) ? 0.42 : lamMin
		String mat3x3List=WaveList_NxN(3,3)
		String recipLattice = SelectString(exists("recip")==1,"","recip")
		hklMax = (numtype(hklMax) || hklMax<1)? 6 : hklMax
		Prompt hklMax, "max value for h,k, or l"
		Prompt recipLattice, "3x3 matrix with reciprocal lattice",popup,mat3x3List
		Prompt lamMin, "minimum allowed lambda (�)"
		DoPrompt "choose max", recipLattice, hklMax,lamMin
		Wave/D recip = $recipLattice
		if (V_flag || !WaveExists(recip))
			return NaN
		endif
		printf "   addPredictedHKL(%s,%d,%g)\r",recipLattice,hklMax,lamMin
	endif

	SVAR materialsList=root:Packages:SCD:materialsList
	String material = StringByKey("material",note(recip),"=")
	Variable i=WhichListItem(material,materialsList)
	if (i<0)
		Abort "unknown material '"+material+"'"
	endif
	SVAR BravaisList=root:Packages:SCD:BravaisList
	String funcName = "isLowestOrder"+StringFromList(i,BravaisList)
	FUNCREF isLowestOrderFCC isLowestOrder = $funcName			// function to use to determine lowest hkl

	Variable h,k,l
	Make/N=3/O/D hkl
	SetDrawLayer /K UserFront		// clear the drawing layer

	for (l=0;l<=hklMax;l+=(l>=0))
		l = -l
		for (k=0;k<=hklMax;k+=(k>=0))
			k = -k
			for (h=0;h<=hklMax;h+=(h>=0))
				h = -h
				hkl = {h,k,l}
				if (isLowestOrder(hkl))
					placeOneSpot(recip,hkl,lamMin)
				endif
			endfor
		endfor
	endfor
	KillWaves/Z hkl
End
//Function addPredictedHKL(recip, hklMax,lamMin)
//	Wave/D recip
//	Variable hklMax
//	Variable lamMin				// usually 0.42 �
//	if (!WaveExists(recip) ||  numtype(hklMax) || hklMax<1 || lamMin<0 || numtype(lamMin))
//		lamMin = (lamMin<0 || numtype(lamMin)) ? 0.42 : lamMin
//		String mat3x3List=WaveList_NxN(3,3)
//		String recipLattice = SelectString(exists("recip")==1,"","recip")
//		hklMax = (numtype(hklMax) || hklMax<1)? 6 : hklMax
//		Prompt hklMax, "max value for h,k, or l"
//		Prompt recipLattice, "3x3 matrix with reciprocal lattice",popup,mat3x3List
//		Prompt lamMin, "minimum allowed lambda (�)"
//		DoPrompt "choose max", recipLattice, hklMax,lamMin
//		Wave/D recip = $recipLattice
//		if (V_flag || !WaveExists(recip))
//			return NaN
//		endif
//		printf "   addPredictedHKL(%s,%d,%g)\r",recipLattice,hklMax,lamMin
//	endif
//
//	SVAR materialsList=root:Packages:SCD:materialsList
//	String material = StringByKey("material",note(recip),"=")
//	Variable i=WhichListItem(material,materialsList)
//	if (i<0)
//		Abort "unknown material '"+material+"'"
//	endif
//	SVAR BravaisList=root:Packages:SCD:BravaisList
//	String funcName = "isLowestOrder"+StringFromList(i,BravaisList)
//	FUNCREF isLowestOrderFCC isLowestOrder = $funcName			// function to use to determine lowest hkl
//
//	Variable h,k,l
//	Make/N=3/O/D hkl
//	SetDrawLayer /K UserFront		// clear the drawing layer
//
//	for (l=-hklMax;l<=hklMax;l+=1)
//		for (k=-hklMax;k<=hklMax;k+=1)
//			for (h=-hklMax;h<=hklMax;h+=1)
//				hkl = {h,k,l}
//				if (isLowestOrder(hkl))
//					placeOneSpot(recip,hkl,lamMin)
//				endif
//			endfor
//		endfor
//	endfor
//	KillWaves/Z hkl
//End


Function placeOneSpot(recip,hkl,lamMin)	// place hkl spots on top image according to reciprocal lattice vectors
	// coordinate system is the standard beamline coordinate system, z-along beam, y-up, x-sideways
	Wave/D recip				// a 3x3 reciprocal lattice a*=recip[][0], b*=recip[][1], c*=recip[][2]
	Wave/D hkl					// hkl vector
	Variable lamMin			// minimum allowed wavelength

	Variable lambda = hkl2lambda(recip,hkl[0],hkl[1],hkl[2])
	if (lambda<lamMin || numtype(lambda))
		return NaN				// wavelength too small
	endif

	String graphName = StringFromList(0,WinList("*",";","WIN:"))
	if (strlen(graphName)<1)
		Abort "No Window for doing a slice update"
	endif
	Wave image = ImageNameToWaveRef(graphName,StringFromList(0,ImageNameList(graphName,";")))
	if (!WaveExists(image))
		return NaN
	endif

	Variable detX, detY						// position on detector (mm)
	calcXYofHKL(recip,hkl,detX,detY,note(image))	// find xy on the detector
	if (numtype(detX) || numtype(detY))
		return NaN
	endif

	Variable xlo,xhi,ylo,yhi
	xlo = DimOffset(image,0)
	ylo = DimOffset(image,1)
	xhi = DimOffset(image,0) + (DimSize(image,0)-1)*DimDelta(image,0)
	yhi = DimOffset(image,1) + (DimSize(image,1)-1)*DimDelta(image,1)
	if (xlo>detX || detX>xhi || ylo>detY || detY>yhi)	// spot not on detector
		return NaN
	endif
	printf "drawing for  (%2d %2d %2d)  at %4.2f (�)   at  detector=(%.2f,%.2f)\r",hkl[0],hkl[1],hkl[2],lambda,detX,detY

	Variable boxSize=3
	SetDrawLayer UserFront
//	SetDrawEnv xcoord= bottom,ycoord= left,fillpat= 0
	SetDrawEnv xcoord= bottom,ycoord= left,linepat= 3,fillpat= 0
	DrawRect detX-boxSize,detY+boxSize, detX+boxSize,detY-boxSize		//	DrawRect left, top, right, bottom

	String str
	if (detY>-100)
		SetDrawEnv xcoord=bottom,ycoord=left,fsize=9,textxjust=1,textyjust=2
	else
		SetDrawEnv xcoord=bottom,ycoord=left,fsize=9,textxjust=1,textyjust=0
	endif

	sprintf str, "(%d %d %d)",hkl[0],hkl[1],hkl[2]
	DrawText detX,detY,str												// DrawText x0, y0, textStr
End




Function calcXYofHKL(recip,hkl,detX,detY,noteStr)	// calc place of hkl on detector from reciprocal lattice vectors
	// coordinate system is the standard beamline coordinate system, z-along beam, y-up, x-sideways
	Wave/D recip				// a 3x3 reciprocal lattice a*=recip[][0], b*=recip[][1], c*=recip[][2]
	Wave/D hkl					// hkl vector
	Variable &detX, &detY		// position on detector (mm), passed back
	String noteStr
	detX = NaN ; detY = NaN		// default is an error

	Variable G					// length of G vector
	MatrixMultiply recip, hkl	// result is in M_product
	KillWaves/Z Gvec
	Rename M_product Gvec		// G vector (from hkl)
	Wave Gvec=Gvec
	Redimension/N=3 Gvec
	Duplicate/O/D Gvec Ghat
	G = normalize(Ghat)		// Ghat is now unit vector in Gvec direction, and G = |Gvec|

	Make/O/D/N=3 kihat={0,0,1}, kfhat, detCentHat, inDetPlane, detXhat
	Variable k = -G/2 / MatrixDot(Ghat,kihat)	// = 2�/lambda,   -ki^ � G^ = sin(theta)

	// to solve this remember that kf = ki + Gvec
	kfhat = (k*kihat[p]) + Gvec[p]				// actually I wil only need kfhat, not kf
	normalize(kfhat)

	Variable deta = valFromWaveNoteOrGlobal(noteStr,"deta")
	Variable detd = valFromWaveNoteOrGlobal(noteStr,"detd")
	Variable cosDeta=cos(deta*PI/180), sinDeta=sin(deta*PI/180)
	detCentHat = {sinDeta, 0, cosDeta}			// unit vector from (000) to center of detector
	Variable dot = MatrixDot(detCentHat,kfhat)
	if (dot<0.01)
		return NaN								// kfhat is does not go towards detector (parallel or worse)
	endif
	Variable d = detd /dot						// distance from origin to detector of scattered ray (mm)
	// vector from detector center to pixel = (origin to pixel) - (origin to detector center)
	inDetPlane = d*kfhat[p] - detd*detCentHat[p]	// detd/dot is distance from origin pixel on detector (mm)
	detXhat = {-cosDeta,0,sinDeta}				// unit vector from detector center toward +x on detector
	detX = MatrixDot(inDetPlane,detXhat)
	detY = inDetPlane[1]
	KillWaves/Z Gvec, Ghat, kihat, kfhat, detCentHat, inDetPlane, detXhat
	return d
End



Function spots2recip(x0,y0,x1,y1)				// convert 2 spots on detector to reciprocal lattice
	Variable x0,y0,x1,y1
	x0 = numtype(x0) ? hcsr(A) : x0			// position on detector (mm), in detector frame
	y0 = numtype(y0) ? vcsr(A) : y0
	x1 = numtype(x1) ? hcsr(B) : x1
	y1 = numtype(y1) ? vcsr(B) : y1
	Make/N=3/O/D hkl0,hkl1, hkl2
	if (numtype(x0) || numtype(y0) || numtype(x1) || numtype(y1))
		Abort "need cursors on spots"
	endif

	SVAR materialsList=root:Packages:SCD:materialsList
	String material
	Prompt material, "material to index", popup,  materialsList
	DoPrompt "pick material", material
	if (V_Flag)
		return 1
	endif

	Variable i=WhichListItem(material,materialsList)
	if (i<0)
		Abort "unknown material '"+material+"'"
	endif
	SVAR aoList=root:Packages:SCD:aoList
	Variable ao = str2num(StringFromList(i,aoList))	// d-spaceing  (�)

	Variable h,k,l
	Prompt h, "h"
	Prompt k, "k"
	Prompt l, "l"
	h=0 ; k=2 ; l=2
	DoPrompt "hkl of first point (cursor A)",h,k,l
	if (V_flag)
		return NaN
	endif
	hkl0 = {h,k,l}

	h=1 ; k=3 ; l=3
	DoPrompt "hkl of second point (cursor B)",h,k,l
	if (V_flag)
		return NaN
	endif
	hkl1 = {h,k,l}

	String noteStr = note(ImageNameToWaveRef("",StringFromList(0,ImageNameList("",";"))))
	String recipName
	sprintf recipName, "recip_%s_%d", material, 10*NumberByKey("depth",noteStr,"=")

	// I have all of the info I need here, construct the reciprocal lattice.

	Make/O/D/N=3 kihat={0,0,1}, kfhat, G0,G1,G2
	// kihat points down stream, from source to sample
	Variable deta = valFromWaveNoteOrGlobal(noteStr,"deta")
	Variable detd = valFromWaveNoteOrGlobal(noteStr,"detd")
	Variable bx = valFromWaveNoteOrGlobal(noteStr,"bx")
	Variable xp, zp								// pixel position (in horizontal plane)

	pixelXZfromDetector(deta,detd,x0,xp,zp)	// calculate pixel coordinates xp,zp
//	kfhat = {xp,y0,zp}							// vertical height (y-direction) is y0
	kfhat = {xp-bx,y0,zp}						// vertical height (y-direction) is y0
	normalize(kfhat)
	G0 = kfhat-kihat							// points half way between kf and -ki,  normalize later

	pixelXZfromDetector(deta,detd,x1,xp,zp)	// calculate pixel coordinates xp,zp
//	kfhat = {xp,y1,zp}							// vertical height (y-direction) is y1
	kfhat = {xp-bx,y1,zp}						// vertical height (y-direction) is y1
	normalize(kfhat)
	G1 = kfhat-kihat

	cross(G0,G1,G2)							// vector cross product to get third direction
	cross(hkl0,hkl1,hkl2)
RemoveCommonFactors(hkl0)
RemoveCommonFactors(hkl1)
RemoveCommonFactors(hkl2)
printwave(hkl0)
printwave(hkl1)
printwave(hkl2)
printwave(G0)
printwave(G1)
printwave(G2)
printf "creating the reciprocal lattice in wave '%s'\r",recipName
print ""
	cross(G2,G0,G1)						// ensure an orthonormal system
	cross(hkl2,hkl0,hkl1)

doAlert 1, "do the funky 'G2 *= -1'?"
if(V_Flag==1)
	G2 *= -1								// this seems to be necessary, but I do not know why !!!!!!
	print "     did the funky  'G2 *= -1'"
else
	print "     took out the 'G2 *= -1'"
endif

	normalize(G0)
	normalize(G1)
	normalize(G2)
	normalize(hkl0)
	normalize(hkl1)
	normalize(hkl2)
printwave(hkl0)
printwave(hkl1)
printwave(hkl2)
printwave(G0)
printwave(G1)
printwave(G2)

	// now G0,G1,G2 and hkl0,hkl1,hkl2 are all defined as orthonormal systems
	// form a matrix of the hklN
	// recip x H = G  -->  recip = G x H^-1

	Make/N=(3,3)/D/O hklN, ident, $recipName
	Wave/D recip = $recipName
	noteStr = ReplaceStringByKey("material",noteStr,material,"=")
	Note/K recip
	Note recip noteStr
	ident = (p==q)							// identity matrix
	hklN[][0] = hkl0[p]
	hklN[][1] = hkl1[p]
	hklN[][2] = hkl2[p]

print "determinant(hklN)= ",MatrixDet(hklN)

	MatrixLinearSolve /O hklN  ident		// ident now contains the inverse of hklN
	KillWaves/Z W_IPIV

	recip[][0] = G0[p]
	recip[][1] = G1[p]
	recip[][2] = G2[p]
print "determinant(G012)= ",MatrixDet(recip)
	MatrixMultiply recip, ident				// multiplying G x H^-1
	Wave M_product=M_product
	recip = M_product * (2*PI/ao)

print "determinant(recip)= ",MatrixDet(recip)


hkl0 = {0,  1,  2}
hkl1 = {1,  2,  2}
hkl2 = {-2,  2,  -1}

print ""
MatrixMultiply recip, hkl0
printWave(hkl0)
normalize(M_product)
printWave(M_product)

MatrixMultiply recip, hkl1
printWave(hkl1)
normalize(M_product)
M_product *= 1.60552
printWave(M_product)

MatrixMultiply recip, hkl2
printWave(hkl2)
normalize(M_product)
printWave(M_product)
	KillWaves/Z G0,G1,G2, hkl0,hkl1,hkl2, M_product,kfhat,kihat, hklN, ident
End


// ====================================================================
// ==========================  End of Reciprocal Lattice   =======================
// ====================================================================
// ====================================================================




// ====================================================================
// ====================================================================
// =============================  Start of Movies   ==========================
// ====================================================================


Function RunAsMovie()
	GetWindow $StringFromList(0,WinList("*",";","WIN:")), note
	Variable oldVal = NumberByKey("value",S_value,"=")
	upDateSliceProc("",-inf,1)				// goto the start
	DoUpdate
	Sleep/T 30
	Variable i = 0
	do
		Sleep /T 4
		i = upDateSliceProc("up",1,1) ? i+1 : i
		DoUpdate
	while (i<2)
	Sleep /T 30
	upDateSliceProc("",oldVal,1)			// goto the original position
End


Function makeMovieOfTopSliceGraph(upDown)	//each frame is at one step in the slice graph
	Variable upDown							// direction of movie, true = up, false=down
	if (numtype(upDown))
		upDown = 1
		Prompt upDown, "direction", popup,"up;down"
		DoPrompt "direction of movie", upDown
		if (V_flag)
			return 1
		endif
		upDown = (upDown==1)
	endif

	String winNam=StringFromList(0,WinList("*",";","WIN:"))
	if (strlen(winNam)<1)
		DoAlert 0, "No Window to make movie of"
		return 1
	endif
	GetWindow $winNam, note
	String noteStr = S_value
	if (NumberByKey("sliceWindow",noteStr,"=")!=1)
		DoAlert 0, "There is no slice window as top graph"
		return 1
	endif

	if (upDown)
		upDateSliceProc("",-inf,1)				// goto the start
	else
		upDateSliceProc("",inf,1)				// goto the start
	endif
	String source = StringByKey("source",noteStr,"=")
	printf "Creating a movie named '%s.mov'\r",source
	NewMovie /F=10/P=home/Z as source+".mov"
	if (V_flag)
		DoAlert 0, "Could not create movie"
		return 1
	endif
	Variable i = 0
	do
		DoUpdate
		AddMovieFrame
		if (upDown)
			i = upDateSliceProc("up",1,1) ? i+1 : i
		else
			i = upDateSliceProc("down",1,1) ? i+1 : i
		endif
	while (i<2)
	CloseMovie
End


// ====================================================================
// ==============================  End of Movies   ==========================
// ====================================================================
// ====================================================================





// ====================================================================
// ====================================================================
// =============================  Start of hkl annotations   =====================
// ====================================================================


// this form is to be called from a marquee around a spot
Function EvaluateSpot(): GraphMarquee
	GetMarquee/K left, bottom
	if (!V_flag)
		return NaN
	endif
	Wave image = ImageNameToWaveRef(S_marqueeWin,StringFromList(0,ImageNameList(S_marqueeWin,";")))
	if (!WaveExists(image))
		return NaN
	endif
	//	String format = "flag: %g; left: %g; top: %g; right: %g; bottom: %g     window:'%s'\r"
	//	printf format, V_flag, V_left, V_top, V_right, V_bottom, S_marqueeWin
	Make/O/T/N=2 T_Constraints
	T_Constraints[0] = {"K6 > 0","K6 < 1"}
	Make/D/O W_coef={0,0,0,0,0,0,0.5}
	Variable V_FitError, V_fitOptions=4
	CurveFit/Q/N Gauss2D kwCWave=W_coef, image(V_left,V_right)(V_bottom,V_top) /C=T_Constraints 
	KillWaves/Z T_Constraints
	if (V_FitError)
		DoAlert 0, "could not fit peak"
		return NaN
	endif

	Variable tolerance=.13						// tolerated error in sqrt(h^2 + k^2 + l^2)
	Variable detX = w_coef[2]					// position on detector (mm)
	Variable detY = w_coef[4]
	Variable tth = calc2theta(detX,detY)
	Variable lambda = getLambdaFromPlot("")
	Variable Integral = integralGauss2d(V_left,V_right,V_bottom,V_top,w_coef[0],image)
	if (numtype(lambda) || tth>=PI || tth<0)	// no wavelength, just get 2theta
//		printf "at detector=(%.3f,%.3f),  2theta = %.4fr",detX,detY,tth*180/PI
		printf "at detector=(%.3f,%.3f),  2theta = %.4f�,    Integral = %g\r",detX,detY,tth*180/PI,Integral
		return 1
	endif
	KillWaves/Z W_coef,W_sigma,W_ParamConfidenceInterval

	SVAR materialsList=root:Packages:SCD:materialsList
	String material="all"
	if (ItemsInList(materialsList)>1)
		Prompt material, "material to index", popup,  materialsList+";-;all"
		DoPrompt "material", material
		if (V_flag)
			return 1
		endif
	endif
	String str,result
	Variable i
	if (!cmpstr("all",material))				// check all the materials in the materialsList
		Variable m,N = ItemsInList(materialsList)
		String mat
		for (m=0;m<N;m+=1)
			mat = StringFromList(m,materialsList)
			printf "for '%s' at %.3g �,   (Integral = %g):\r",mat,lambda,Integral
			result = calcHKL(detX,detY,mat,tolerance)
			i = 0
			do
				str = StringFromList(i,result)
				if (strlen(str)<1)
					break
				endif
				print str
				i += 1
			while(1)
		endfor
		return 0
	endif

	// just a single material, first print out the results
	printf "for '%s' at %.3g �,   (Integral = %g):\r",material,lambda,Integral
	result = calcHKL(detX,detY,material,tolerance)
	String list = ""
	Variable h,k,l,del,delMin=1e10
	i = 0
	do
		str = StringFromList(i,result)
		if (strlen(str)<1)
			break
		endif
		print str
		sscanf str, "(%d,%d,%d)", h,k,l
		del = str2num(str[strsearch(str,"� = ",0)+4,inf])
		delMin = min(delMin,del)
		sprintf str,"hkl=(%d %d %d),del=%g;", h,k,l,del
		list += str
		i += 1
	while(1)

	String item,text = ""
	for (i=0;i<ItemsInList(list); i+=1)
		item = StringFromList(i,list)
		if (delMin==NumberByKey("del",item,"=",","))
			if (strlen(text)>0)
				text += "\r"
			endif
			text += StringByKey("hkl",item,"=",",")
		endif
	endfor

	DoAlert 1, "add label to plot?"
	if (V_flag==1)
		i = WhichListItem(material,materialsList)
		str = "\\K("+SelectString(i-1,"65535,0,0","0,0,65535","0,65535,0")+")"
		addAnnotaionHKL(S_marqueeWin,detX,detY,str+text)
	endif
	if (exists("recip"))
		Wave/D recip=recip
		printf "using reciprocal lattice, hkl = (%s)\r",hkl_fromDetectorSpot_recip(detX,detY,recip)
	endif
End
Function integralGauss2d(V_left,V_right,V_bottom,V_top,K0,image)
	Variable V_left,V_right,V_bottom,V_top		// size of box in axis units (mm)
	Variable K0
	Wave image

	Variable ilo = (V_left-DimOffset(image,0)) / DimDelta(image,0)
	Variable ihi = (V_right-DimOffset(image,0)) / DimDelta(image,0)
	Variable jlo = (V_bottom-DimOffset(image,1)) / DimDelta(image,1)
	Variable jhi = (V_top-DimOffset(image,1)) / DimDelta(image,1)
	Variable i,j, N=0, Int=0
	for (j=jlo;j<=jhi;j+=1)
		for (i=ilo;i<=ihi;i+=1)
			Int += image[i][j]
			N += 1
		endfor
	endfor
	Int -= N*K0										// remove background
	Int *=  DimDelta(image,0)* DimDelta(image,1)		// scale by dx,dy to make it an area
	return Int
	//	K0+K1*exp((-1/(2*(1-K6^2)))*(((x-K2)/K3)^2 + ((y-K4)/K5)^2 - (2*K6*(x-K2)*(y-K4)/(K3*K5))))
End




// this form is to be called from a menu, or command line
Function InfoAboutSpotOnDetector(detX,detY,tolerance)
	Variable detX, detY						// postion on detector (mm), or NaN to use cursor A
	Variable tolerance=.1					// tolerated error in sqrt(h^2 + k^2 + l^2)
	Variable tth = calc2theta(detX,detY)
	printf "at detector=(%.3f,%.3f),  2theta = %.4fr",detX,detY,tth*180/PI
	if (numtype(getLambdaFromPlot("")))	// good wavelength, try for an hkl
		return tth
	endif

	SVAR materialsList=root:Packages:SCD:materialsList
	String material
	Prompt material, "material to index", popup,  materialsList
	DoPrompt "material", material
	if (V_flag)
		return 1
	endif
	String str = calcHKL(detX,detY,material,tolerance)
	if (strlen(str)<1)
		print "no hkl close enough to point on detector"
		return tth
	endif
	Variable i
	for (i=0;i<ItemsInList(str);i+=1)
		print StringFromList(i,str)
	endfor
	return tth
End



Function/T calcHKL(detX,detY,material,tolerance)	// get possible hkl's from the detector position and wavelength
	Variable detX,detY
	String material
	Variable tolerance=.1							// tolerated error in sqrt(h^2 + k^2 + l^2)

	SVAR materialsList=root:Packages:SCD:materialsList
	if (WhichListItem(material,materialsList)<0)	// material not in materialsList
		Prompt material, "material to index", popup,  materialsList
		DoPrompt "pick material", material
	endif

	Variable i=WhichListItem(material,materialsList)
	if (i<0)
		Abort "unknown material '"+material+"'"
	endif
	SVAR aoList=root:Packages:SCD:aoList
	SVAR BravaisList=root:Packages:SCD:BravaisList
	Variable ao = str2num(StringFromList(i,aoList))	// d-spaceing  (�)
	String funcName = StringFromList(i,BravaisList)+"allow"
	FUNCREF FCCallow allow = $funcName				// function to use to determine hkl is allowed

	Variable lambda = getLambdaFromPlot("")
	Variable tth = calc2theta(detX,detY)
	Variable d = lambda / sin(tth/2) / 2				// d-spaceing (�)
	if (numtype(d) || d<=0)
		return ""
	endif
	Variable hklMag = ao/d								// sqrt(h2+k^2+l^2)
	String str, answer=""
	Variable i1,i2,i3, del, count=0
	for (i3=0;i3<11;i3+=1)
		for (i2=i3;i2<11;i2+=1)
			for (i1=i2;i1<11;i1+=1)
				del = abs(hklMag-sqrt(i3^2+i2^2+i1^2))
				if (del<tolerance && allow(i3,i2,i1))
					sprintf str, "(%d,%d,%d) at [%.3f, %.3f, %.4f],    � = %.3f",i3,i2,i1,detX,detY,lambda,del
					answer += str+";"
					count += 1
				endif
			endfor
		endfor
	endfor
	if (count==0)
		printf "for %s, no hkl matches hkl^2 = %g to a tolerance of %g, 2theta = %g, d-spacing=%g\r",material,hklMag^2,tolerance,tth*180/PI,d
	endif
	return answer
End
Function FCCallow(h,k,l)					// all odd or all even
	Variable h,k,l
	Variable s = mod(abs(h),2)+mod(abs(k),2) + mod(abs(l),2)
	return (s==0 || s==3)
End
Function BCCallow(h,k,l)					// sum must be even
	Variable h,k,l
	Variable s = abs(h+k+l)
	return !(mod(s,2))
End
Function calc2theta(detX,detY)				// compute 2theta (rad) from the location on the neutron detector
	Variable &detX,&detY
	detX = numtype(detX) ? hcsr(A) : detX
	detY = numtype(detY) ? vcsr(A) : detY
	String noteStr = note(ImageNameToWaveRef("",StringFromList(0,ImageNameList("",";"))))
	Variable deta = valFromWaveNoteOrGlobal(noteStr,"deta")
	Variable detd = valFromWaveNoteOrGlobal(noteStr,"detd")
	Variable bx = valFromWaveNoteOrGlobal(noteStr,"bx")
	Variable xp, zp, yp=detY						// y is vertical
	pixelXZfromDetector(deta,detd,detX,xp,zp)		// calculate pixel coordinates xp,zp
	xp -= bx
	Variable tth = atan2(sqrt(xp*xp+yp*yp),zp)	// 2-theta (radians)
	return tth
End
//Function calc2theta(detX,detY)				// compute 2theta (rad) from the location on the neutron detector
//	Variable &detX,&detY
//	detX = numtype(detX) ? hcsr(A) : detX
//	detY = numtype(detY) ? vcsr(A) : detY
//	String noteStr = note(ImageNameToWaveRef("",StringFromList(0,ImageNameList("",";"))))
//	Variable deta = valFromWaveNoteOrGlobal(noteStr,"deta")
//	Variable detd = valFromWaveNoteOrGlobal(noteStr,"detd")
//	Variable xp, zp, yp=detY						// y is vertical
//	pixelXZfromDetector(deta,detd,detX,xp,zp)		// calculate pixel coordinates xp,zp
//	Variable tth = atan2(sqrt(xp*xp+yp*yp),zp)	// 2-theta (radians)
//	return tth
//End
Function getLambdaFromPlot(graphName)			// get lambda from a step slice plot
	String graphName								// if empty string, use top plot
	if (strlen(graphName)<1)
		graphName = StringFromList(0,WinList("*",";","WIN:"))
	endif
	if (strlen(graphName)<1)
		DoAlert 0, "No Window for doing a slice update"
		return 1
	endif
	GetWindow $graphName, note
	Variable lambda = NumberByKey("lambda",S_value,"=")
	return lambda
End



// Determine the angle between two vectors with coordinates (h1,k1,l1), (h2,k2,l2) (cubic!)
// This is the angle between the reciprocal lattice vectors, not the outgoing wave vectors
Function angleBetweenHKLs(h1,k1,l1,h2,k2,l2)
	Variable h1,k1,l1,h2,k2,l2
	Variable printIt = 0
	if (numtype(h1+k1+l1+h2+k2+l2))
		h1 = numtype(h1) ? 1 : h1
		k1 = numtype(k1) ? 1 : k1
		l1 = numtype(l1) ? 1 : l1
		h2 = numtype(h2) ? 1 : h2
		k2 = numtype(k2) ? 1 : k2
		l2 = numtype(l2) ? 1 : l2
		Prompt h1, "h 1"
		Prompt k1, "k 1"
		Prompt l1, "l 1"
		Prompt h2, "h 2"
		Prompt k2, "k 2"
		Prompt l2, "l 2"
		DoPrompt "two vectors", h1,h2,k1,k2,l1,l2
		if (V_flag)
			return NaN
		endif
		printIt = 1
	endif
	Make/N=3/O/D v1,v2
	v1 = {h1,k1,l1}
	v2 = {h2,k2,l2}
	Variable ang = angleBetweenVectors(v1,v2)*180/PI
	KillWaves/Z v1,v2
	if (printIt)
		printf "(%g, %g, %g) ^ (%g, %g, %g) = %gr",h1,k1,l1,h2,k2,l2,ang
	endif
	return ang
End



// Determine the angle between two cartesian vectors
Function angleBetweenVectors(v1,v2)		// retruns angle between v1 and v2  (radians)
	Wave v1
	Wave v2
	if (DimSize(v1,1)>0 || DimSize(v2,1)>0 || DimSize(v1,0)!=DimSize(v2,0))	// 1-d vectors must match
		return NaN
	endif
	Variable N = DimSize(v1,0)
	Variable i,m1=0,m2=0, dot=0
	for (i=0;i<N;i+=1)
		m1 += v1[i]^2
		m2 += v2[i]^2
		dot += v1[i]*v2[i]
	endfor
	dot /= sqrt(m1*m2)						// first nomalize
	dot = max(-1,min(1,dot))					// ensure that truncation does not |cosine| > 1
	return acos(dot)
End





// ===============================================================


Function addAnnotaionHKL(graphName,xx,yy,str)
	String graphName
	Variable xx,yy				// postion of annotaion
	String str					// string to add
	String list = AnnotationList(graphName)

	String name								// annotaion name
	Variable i=0
	do
		name = "hkl"+num2istr(i)
		i += 1
	while(WhichListItem(name,list)>0)	// this loop finds next name of type hklN

	Wave image = ImageNameToWaveRef(graphName,StringFromList(0,ImageNameList(graphName,";")))
	if (!WaveExists(image))
		return NaN
	endif

	Variable index		// index into image corresponding to xx,yy
	index = round((yy-DimOffset(image,1))/DimDelta(image,1))	// index in y direction
	index *= DimSize(image,0)									// ponits to correct row
	index += round((xx-DimOffset(image,0))/DimDelta(image,0))	// along row
	index = pnt2x(image,index)									// funny thing Igor needs
//	str = "\\Z09\\K(65535,0,0)"+str
	str = "\\Z09"+str
	Tag/W=$graphName/A=MC/B=1/F=0/X=0/Y=0/N=$name/I=1/L=0/Z=1 $NameOfWave(image),index,str
End


Function deleteAnnotationsHKL(): GraphMarquee
	GetMarquee/K left, bottom
	if (!V_flag)
		return NaN
	endif

	String list = AnnotationList(S_marqueeWin)
	String name
	Variable i, N=ItemsInList(list)
	for (i=0;i<N;i+=1)
		name = StringFromList(i,list)
		if (strsearch(name,"hkl",0)==0)
			TextBox/K/N=$name
		endif
	endfor
End
Function hideAnnotationsHKL(): GraphMarquee
	GetMarquee/K left, bottom
	if (!V_flag)
		return NaN
	endif

	String list = AnnotationList(S_marqueeWin)
	String name
	Variable i, N=ItemsInList(list)
	for (i=0;i<N;i+=1)
		name = StringFromList(i,list)
		if (strsearch(name,"hkl",0)==0)
			Tag/C/N=$name/V=0
		endif
	endfor
End
Function showAnnotationsHKL(): GraphMarquee
	GetMarquee/K left, bottom
	if (!V_flag)
		return NaN
	endif

	String list = AnnotationList(S_marqueeWin)
	String name
	Variable i, N=ItemsInList(list)
	for (i=0;i<N;i+=1)
		name = StringFromList(i,list)
		if (strsearch(name,"hkl",0)==0)
			Tag/C/N=$name/V=1
		endif
	endfor
End


// ====================================================================
// ==============================  End of hkl annotations   =====================
// ====================================================================
// ====================================================================




// ====================================================================
// ====================================================================
// ==============================  Sorting and Summing   =====================
// ====================================================================


// This routine makes a 3-d array, where each pixel position (x,y) is the thresholded sum over lambda 
// for one depth.  The z coordinate is the reconstructed depth.  So as you step through z, you see the 
// lambda summed laue pattern for each depth.
Function MakeSlicedSumsAtManyDepths(base,range,depth0,depthMax,delta,sumThresh)
	String base = "stackA"			// base name for raw data
	String range = "8523-8538"	// run range to include
	Variable depth0					// first depth (mm)
	Variable depthMax				// last depth (mm)
	Variable delta					// step size in depth (mm)
	Variable sumThresh=5			// for summing in angle, only include pixels � sumThresh

	if (numtype(depth0+depthMax+delta+sumThresh) || sumThresh<0 || strlen(base)<1 || strlen(range)<1)
		depth0 = numtype(depth0) ? 0 : depth0
		depthMax = numtype(depthMax) ? 17 : depthMax
		delta = numtype(delta) ? 1 : delta
		sumThresh = numtype(sumThresh) ? 5 : sumThresh
		base = SelectString(strlen(base) , "stackA", base)
		SVAR rangeList = root:Packages:SCD:rangeList
		range = SelectString(strlen(range),StringFromList(ItemsInList(rangeList)-1,rangeList),range)
		Prompt depth0, "first depth to process (mm)"
		Prompt depthMax, "last depth to process (mm)"
		Prompt delta, "step size in depth (mm)"
		Prompt base, "base name of raw data wave name, e.g. base='stackA'"
		Prompt range, "range of run numbers in original data set, e.g. '8523-8538'",popup,rangeList+";other"
		Prompt sumThresh, "threshold for summing each depth in time"
		DoPrompt "depth summing",base,range,depth0,depthMax,delta,sumThresh
		if (V_flag)
			return 1
		endif
		if (!cmpstr(range,"other"))
			Prompt range, "range of run numbers, e.g. '3-10,12'"
			DoPrompt "new range", range
			if (V_flag)
				return 1
			endif
			rangeList = AddListItem(range,rangeList,";",Inf)
		endif
		printf "MakeSlicedSumsAtManyDepths(\"%s\",\"%s\",%g,%g,%g,%d)\r",base,range,depth0,depthMax,delta,sumThresh
	endif
	if (numtype(depth0+depthMax+delta+sumThresh) || sumThresh<0 || strlen(base)<1 || strlen(range)<1)
		Abort "illegal values on input"
	endif
	Variable start1, start = ticks

	Variable depth
	Variable tol = delta/2
	String allSums = base+"AllSums"
	String sliceName, sumName="", noteStr
	Variable i, Ndepths=floor((depthMax-depth0)/delta +1)
	printf "storing all the sum slices in '%s'\r",allSums
	Variable once=1

	for (i=0,depth=depth0;depth<=depthMax;depth+=delta,i+=1)
		start1 = ticks
		sliceName = SortDataOneDepth(depth,tol,base,range,sumThresh)
//		printf "  sorted '%s'  at depth=%g,  in %.2f sec\r",sliceName,NumberByKey("depth",note($sliceName),"="),(ticks-start1)/60.15
		NVAR V_nonZero=V_nonZero
		if (V_nonZero)
			sumName = sumUpWavelength(sliceName,"0-"+num2istr(DimSize($sliceName,2)-1),sumThresh)
//			sumName = sumUpWavelength(sliceName,"0-225,227-271",sumThresh)
			Wave asum = $sumName
			if (once)
				once = 0
				Duplicate/O $sumName $allSums
				Redimension/N=(DimSize(asum,0),DimSize(asum,1),Ndepths) $allSums
				SetScale/P z, depth0, delta, "mm", $allSums
				Wave all = $allSums
				noteStr = note(all)
				noteStr = RemoveByKey("depth",noteStr,"=")
				Note /K all
				Note all, noteStr
			endif
			all[][][i] = asum[p][q]
		endif
		printf "  sorted and summed '%s'  at depth=%g,  in %.2f sec\r",sliceName,NumberByKey("depth",note($sliceName),"="),(ticks-start1)/60.15
		KillWaves/Z $sliceName, $sumName
	endfor
	printf "sorting and summing %d depths took %s\r",i,num2sexigesmal((ticks-start)/60.15,0)
	beep
	return 0
End



Function/T SortDataOneDepth(depth,tol,base,range,thresh)	// sort data into 1 depth array
	// sort data and save only one depth at the given depth (the result is 3-d, x y & lambda)
	// returns the full path name of the new 3-d wave
	//	Note that +z is along beam, +y is up, and +x is horizontally perp to beam (points to where you stand).
	Variable depth						// target depth to get (mm)
	Variable tol							// tolerance on the depth (mm), typically 0.5
	String base							// for file name, e.g. base="CuSCD", or "stackA"
	String range						// range of run numbers in original data set
	Variable thresh						// only include pixels that are >= thresh, 0 includes all
//	Variable res=1.036/2				// resolution half-width (mm), assumed flat topped

	Variable printIt=0
	if (numtype(depth+tol+thresh) || tol<=0 || thresh<0 || strlen(base)<1 ||  strlen(range)<1)
		depth = numtype(depth) ? 0 : depth
		tol = (numtype(tol) || tol<=0) ? 0.5 : tol
		thresh = numtype(thresh) ? 0 : abs(thresh)
		base = SelectString(strlen(base) , "stackA", base)
		SVAR rangeList = root:Packages:SCD:rangeList
		range = SelectString(strlen(range),StringFromList(ItemsInList(rangeList)-1,rangeList),range)
		Prompt depth, "desired depth in the sample (mm)"
		Prompt tol, "tolerance on the depth (depth眛ol) (mm)"
		Prompt thresh, "threshold for including a pixel, 0 includes all"
		Prompt base, "for wave name, e.g. base='stackA'"
		Prompt range, "range of run numbers in original data set, e.g. '8523-8538'",popup,rangeList+";other"
		DoPrompt "the depth to sort out", depth, tol,base,range,thresh
		if (V_flag)
			return ""
		endif
		if (!cmpstr(range,"other"))
			Prompt range, "range of run numbers, e.g. '3-10,12'"
			DoPrompt "new range", range
			if (V_flag)
				return ""
			endif
			rangeList = AddListItem(range,rangeList,";",Inf)
		endif
		printIt = 1
	endif
	thresh = round(thresh)
	if (numtype(depth+tol+thresh) || thresh<0 || tol<=0 || strlen(base)<1 || strlen(range)<1)
		print "Illegal depth or tolerance"
		return ""
	endif
	String sliceName
	sprintf sliceName "%s_%s%02d_%02d",base, SelectString(depth<0,"","m"),abs(depth*10),10*tol
	if (printIt)
		printf "SortDataOneDepth(%g,%g,\"%s\",\"%s\",%g)\t\t\t// into: '%s'\r",depth,tol,base,range,thresh,sliceName
	endif

	range = expandRange(range,";")
	String frameName
	sprintf frameName, "root:raw:%s%04d",base,str2num(StringFromList(0,range))
	Wave frame = $frameName
	Variable dimX=DimSize(frame,0),dimY=DimSize(frame,1),dimZ=DimSize(frame,2)
	String noteStr = note(frame)

	Make/N=(dimX,dimY,dimZ)/O/U/W $sliceName	// create 1 depth slice for x, y, and lambda
	Wave slice = $sliceName
	slice = 0
	SetXYtoANGER(slice)
	noteStr = ReplaceNumberByKey("depth",noteStr,depth,"=")
	noteStr = ReplaceNumberByKey("tol",noteStr,tol,"=")
	noteStr = ReplaceNumberByKey("threshSort",noteStr,thresh,"=")
	noteStr = RemoveByKey("translate_mm",noteStr,"=")
	Note slice, noteStr

	Variable deta = valFromWaveNoteOrGlobal(noteStr,"deta")
	Variable detd = valFromWaveNoteOrGlobal(noteStr,"detd")
	Variable x2mm = valFromWaveNoteOrGlobal(noteStr,"x2mm")
	Variable y2mm = valFromWaveNoteOrGlobal(noteStr,"y2mm")
	Variable xLeftmm = valFromWaveNoteOrGlobal(noteStr,"xLeftmm")
	Variable slitx = valFromWaveNoteOrGlobal(noteStr,"slitx")
	Variable slitz = valFromWaveNoteOrGlobal(noteStr,"slitz")
	Variable bx = valFromWaveNoteOrGlobal(noteStr,"bx")
	Variable surface = valFromWaveNoteOrGlobal(noteStr,"surface")

	Variable k								// index over each of the runs in range
	Variable i								// index over each pixels in the x-direcion in the detector
	Variable trans							// translation position for each SCD run
	Variable xi								// x position on detector from its center (mm)
	Variable xp, zp							// x,z coodinates of the pixel (mm), in beam line system
	Variable zintercept						// position on beam axis where beam hits (mm)
//	Variable x0=DimOffset(slice,0)
//	Variable dx=DimDelta(slice,0)
//	Variable dmin, dmax					// range of the overlap
	Variable/G V_nonZero=0				// flags if this depth is empty
	Variable angle							// angle between detector normal and exit ray
	Variable detaR = deta*PI/180			// radian version of deta
	Variable iorigin							// index where beam would land if this volume were at (0,0)
	Variable start = ticks
	do										// look at each run in the range
		if (strlen(StringFromList(k,range))<1)
			break
		endif
		sprintf frameName, "root:raw:%s%04d",base,str2num(StringFromList(k,range))
		Wave frame = $frameName
		if (!WaveExists(frame))
			Abort "the input wave '"+frameName+"' cannot be found"
		endif
		trans = NumberByKey("translate_mm",note(frame),"=")
		for (i=0;i<dimX;i+=1)
			xi = i*x2mm + xLeftmm					// horiz distance on detector from center
			pixelXZfromDetector(deta,detd,xi,xp,zp)	// calculate pixel coordinates xp,zp
//			zintercept = slitz - slitx*(zp-slitz)/(xp-slitx)
//			zintercept = (slitz*xp-slitx*zp) / (xp-slitx)	// this line is the same as the one before it!!!
			zintercept = ((bx-slitx)*zp - (bx-xp)*slitz) / (xp-slitx)	// added an offset to incident beam
			// surface is z coordinate of surface when trans=0
			// depth = zintercept - trans - surface,  a bigger zintercept is deeper
			//  old			// surface at center, when trans=surface
			//  old			// depth = zintercept + trans - surface,  a bigger zintercept is deeper
			if (abs(zintercept-trans-surface-depth)<=tol)
				V_nonZero = 1									// found something, not empty
				// add to sorted image at the position on CCD for sample at (0,0)
				// ideal Xccd = detd * tan(deta - atan2(kx,kz))
				// for atan2(kx,kz) = atan2((xp-bx),(zp-zi)),  a ray starting at (zi,bx) and ending at (zp,xp)
				angle = detaR - atan2((xp-bx),(zp-zintercept)) // correct to position z==0
				iorigin = (detd*tan(angle)-xLeftmm)/x2mm	// index for diffraction from (0,0)
				if (angle<-PI/2 || angle>PI/2)
					printf "angle = %g�-%g� = %g�,    i = %d,   iorigin = %d,       %s\r",deta,atan2((xp-bx),(zp-zintercept))*180/PI,angle*180/PI,i,iorigin,NameOfWave(frame)
				elseif (thresh<0.5)
					slice[iorigin][][] += frame[i][q][r]
				else
					slice[iorigin][][] += (frame[i][q][r]>thresh) ? frame[i][q][r] : 0
				endif
			endif
//			dmin = max(zmm-depth-res,depth-tol)		// bottom of overlap range
//			dmax = min(zmm-depth+res,depth+tol)	// top of overlap range
//			if (dmax>dmin)
//				slice[i][][] += round((dmax-dmin)/res/2*frame[i][q][r])
//			endif
		endfor
		k += 1
	while(1)
	if (printIt)
		printf "sorting took %.2f sec\r",(ticks-start)/60.15
		DoAlert 1, "put up a plot of "+sliceName
		if (V_flag)
			ShowSliceGraph(sliceName)
		endif
	endif
	return GetWavesDataFolder(slice,2)				// return the full path
End
Function/T SortDataOneDepthOLD(depth,tol,base,range,thresh)	// sort data into 1 depth array
	// sort data and save only one depth at the given depth (the result is 3-d, x y & lambda)
	// returns the full path name of the new 3-d wave
	//	Note that +z is along beam, +y is up, and +x is horizontally perp to beam (points to where you stand).
	Variable depth						// target depth to get (mm)
	Variable tol							// tolerance on the depth (mm), typically 0.5
	String base							// for file name, e.g. base="CuSCD", or "stackA"
	String range						// range of run numbers in original data set
	Variable thresh						// only include pixels that are >= thresh, 0 includes all
//	Variable res=1.036/2				// resolution half-width (mm), assumed flat topped

	Variable printIt=0
	if (numtype(depth+tol+thresh) || tol<=0 || thresh<0 || strlen(base)<1 ||  strlen(range)<1)
		depth = numtype(depth) ? 0 : depth
		tol = (numtype(tol) || tol<=0) ? 0.5 : tol
		thresh = numtype(thresh) ? 0 : abs(thresh)
		base = SelectString(strlen(base) , "stackA", base)
		SVAR rangeList = root:Packages:SCD:rangeList
		range = SelectString(strlen(range),StringFromList(ItemsInList(rangeList)-1,rangeList),range)
		Prompt depth, "desired depth in the sample (mm)"
		Prompt tol, "tolerance on the depth (depth眛ol) (mm)"
		Prompt thresh, "threshold for including a pixel, 0 includes all"
		Prompt base, "for wave name, e.g. base='stackA'"
		Prompt range, "range of run numbers in original data set, e.g. '8523-8538'",popup,rangeList+";other"
		DoPrompt "the depth to sort out", depth, tol,base,range,thresh
		if (V_flag)
			return ""
		endif
		if (!cmpstr(range,"other"))
			Prompt range, "range of run numbers, e.g. '3-10,12'"
			DoPrompt "new range", range
			if (V_flag)
				return ""
			endif
			rangeList = AddListItem(range,rangeList,";",Inf)
		endif
		printIt = 1
	endif
	thresh = round(thresh)
	if (numtype(depth+tol+thresh) || thresh<0 || tol<=0 || strlen(base)<1 || strlen(range)<1)
		print "Illegal depth or tolerance"
		return ""
	endif
	String sliceName
	sprintf sliceName "%s_%s%02d_%02d",base, SelectString(depth<0,"","m"),abs(depth*10),10*tol
	if (printIt)
		printf "SortDataOneDepth(%g,%g,\"%s\",\"%s\",%g)\t\t\t// into: '%s'\r",depth,tol,base,range,thresh,sliceName
	endif

	range = expandRange(range,";")
	String frameName
	sprintf frameName, "root:raw:%s%04d",base,str2num(StringFromList(0,range))
	Wave frame = $frameName
	Variable dimX=DimSize(frame,0),dimY=DimSize(frame,1),dimZ=DimSize(frame,2)
	String noteStr = note(frame)

	Make/N=(dimX,dimY,dimZ)/O/U/W $sliceName	// create 1 depth slice for x, y, and lambda
	Wave slice = $sliceName
	slice = 0
	SetXYtoANGER(slice)
	noteStr = ReplaceNumberByKey("depth",noteStr,depth,"=")
	noteStr = ReplaceNumberByKey("tol",noteStr,tol,"=")
	noteStr = ReplaceNumberByKey("threshSort",noteStr,thresh,"=")
	noteStr = RemoveByKey("translate_mm",noteStr,"=")
	Note slice, noteStr

	Variable deta = valFromWaveNoteOrGlobal(noteStr,"deta")
	Variable detd = valFromWaveNoteOrGlobal(noteStr,"detd")
	Variable x2mm = valFromWaveNoteOrGlobal(noteStr,"x2mm")
	Variable y2mm = valFromWaveNoteOrGlobal(noteStr,"y2mm")
	Variable xLeftmm = valFromWaveNoteOrGlobal(noteStr,"xLeftmm")
	Variable slitx = valFromWaveNoteOrGlobal(noteStr,"slitx")
	Variable slitz = valFromWaveNoteOrGlobal(noteStr,"slitz")
	Variable bx = valFromWaveNoteOrGlobal(noteStr,"bx")
	Variable surface = valFromWaveNoteOrGlobal(noteStr,"surface")

	Variable k								// index over each of the runs in range
	Variable i								// index over each pixels in the x-direcion in the detector
	Variable trans							// translation position for each SCD run
	Variable xi								// x position on detector from its center (mm)
	Variable xp, zp							// x,z coodinates of the pixel (mm), in beam line system
	Variable zintercept						// position on beam axis where beam hits (mm)
//	Variable x0=DimOffset(slice,0)
//	Variable dx=DimDelta(slice,0)
//	Variable dmin, dmax					// range of the overlap
	Variable/G V_nonZero=0				// flags if this depth is empty
	Variable start = ticks
	do										// look at each run in the range
		if (strlen(StringFromList(k,range))<1)
			break
		endif
		sprintf frameName, "root:raw:%s%04d",base,str2num(StringFromList(k,range))
		Wave frame = $frameName
		if (!WaveExists(frame))
			Abort "the input wave '"+frameName+"' cannot be found"
		endif
		trans = NumberByKey("translate_mm",note(frame),"=")
		for (i=0;i<dimX;i+=1)
			xi = i*x2mm + xLeftmm					// horiz distance on detector from center
			pixelXZfromDetector(deta,detd,xi,xp,zp)	// calculate pixel coordinates xp,zp
//			zintercept = slitz - slitx*(zp-slitz)/(xp-slitx)
//			zintercept = (slitz*xp-slitx*zp) / (xp-slitx)	// this line is the same as the one before it!!!
			zintercept = ((bx-slitx)*zp - (bx-xp)*slitz) / (xp-slitx)	// added an offset to incident beam
			// surface is z coordinate of surface when trans=0
			// depth = zintercept - trans - surface,  a bigger zintercept is deeper
			//  old			// surface at center, when trans=surface
			//  old			// depth = zintercept + trans - surface,  a bigger zintercept is deeper
			if (abs(zintercept-trans-surface-depth)<=tol)
				V_nonZero = 1				// found something, not empty
				if (thresh<.5)
					slice[i][][] += frame[i][q][r]
				else
					slice[i][][] += (frame[i][q][r]>thresh) ? frame[i][q][r] : 0
				endif
			endif
//			dmin = max(zmm-depth-res,depth-tol)		// bottom of overlap range
//			dmax = min(zmm-depth+res,depth+tol)	// top of overlap range
//			if (dmax>dmin)
//				slice[i][][] += round((dmax-dmin)/res/2*frame[i][q][r])
//			endif
		endfor
		k += 1
	while(1)
	if (printIt)
		printf "sorting took %.2f sec\r",(ticks-start)/60.15
		DoAlert 1, "put up a plot of "+sliceName
		if (V_flag)
			ShowSliceGraph(sliceName)
		endif
	endif
	return GetWavesDataFolder(slice,2)				// return the full path
End
Function/T SortDataOneDepth_RealOLD(depth,tol,base,range,thresh)	// sort data into 1 depth array
	// sort data and save only one depth at the given depth (the result is 3-d, x y & lambda)
	// returns the full path name of the new 3-d wave
	//	Note that +z is along beam, +y is up, and +x is horizontally perp to beam (points to where you stand).
	Variable depth						// target depth to get (mm)
	Variable tol							// tolerance on the depth (mm), typically 0.5
	String base							// for file name, e.g. base="CuSCD", or "stackA"
	String range						// range of run numbers in original data set
	Variable thresh						// only include pixels that are >= thresh, 0 includes all
	Variable res=1.036/2				// resolution half-width (mm), assumed flat topped

	Variable printIt=0
	if (numtype(depth+tol+thresh) || tol<=0 || thresh<0 || strlen(base)<1 ||  strlen(range)<1)
		depth = numtype(depth) ? 0 : depth
		tol = (numtype(tol) || tol<=0) ? 0.5 : tol
		thresh = numtype(thresh) ? 0 : abs(thresh)
		base = SelectString(strlen(base) , "stackA", base)
		SVAR rangeList = root:Packages:SCD:rangeList
		range = SelectString(strlen(range),StringFromList(ItemsInList(rangeList)-1,rangeList),range)
		Prompt depth, "desired depth in the sample (mm)"
		Prompt tol, "tolerance on the depth (depth眛ol) (mm)"
		Prompt thresh, "threshold for including a pixel, 0 includes all"
		Prompt base, "for wave name, e.g. base='stackA'"
		Prompt range, "range of run numbers in original data set, e.g. '8523-8538'",popup,rangeList+";other"
		DoPrompt "the depth to sort out", depth, tol,base,range,thresh
		if (V_flag)
			return ""
		endif
		if (!cmpstr(range,"other"))
			Prompt range, "range of run numbers, e.g. '3-10,12'"
			DoPrompt "new range", range
			if (V_flag)
				return ""
			endif
			rangeList = AddListItem(range,rangeList,";",Inf)
		endif
		printIt = 1
	endif
	thresh = round(thresh)
	if (numtype(depth+tol+thresh) || thresh<0 || tol<=0 || strlen(base)<1 || strlen(range)<1)
		print "Illegal depth or tolerance"
		return ""
	endif
	String sliceName
	sprintf sliceName "%s_%s%02d_%02d",base, SelectString(depth<0,"","m"),abs(depth*10),10*tol
	if (printIt)
		printf "SortDataOneDepth(%g,%g,\"%s\",\"%s\",%g)\t\t\t// into: '%s'\r",depth,tol,base,range,thresh,sliceName
	endif

	range = expandRange(range,";")
	String frameName
	sprintf frameName, "root:raw:%s%04d",base,str2num(StringFromList(0,range))
	Wave frame = $frameName
	Variable dimX=DimSize(frame,0),dimY=DimSize(frame,1),dimZ=DimSize(frame,2)
	String noteStr = note(frame)

	Make/N=(dimX,dimY,dimZ)/O/U/W $sliceName	// create 1 depth slice for x, y, and lambda
	Wave slice = $sliceName
	slice = 0
	SetXYtoANGER(slice)
	noteStr = ReplaceNumberByKey("depth",noteStr,depth,"=")
	noteStr = ReplaceNumberByKey("tol",noteStr,tol,"=")
	noteStr = ReplaceNumberByKey("threshSort",noteStr,tol,"=")
	Note slice, noteStr

	Variable deta = valFromWaveNoteOrGlobal(noteStr,"deta")
	Variable detd = valFromWaveNoteOrGlobal(noteStr,"detd")
	Variable x2mm = valFromWaveNoteOrGlobal(noteStr,"x2mm")
	Variable y2mm = valFromWaveNoteOrGlobal(noteStr,"y2mm")
	Variable xLeftmm = valFromWaveNoteOrGlobal(noteStr,"xLeftmm")
	Variable slitx = valFromWaveNoteOrGlobal(noteStr,"slitx")
	Variable slitz = valFromWaveNoteOrGlobal(noteStr,"slitz")
	Variable surface = valFromWaveNoteOrGlobal(noteStr,"surface")

	Variable k								// index over each of the runs in range
	Variable i								// index over each pixels in the x-direcion in the detector
	Variable trans							// translation position for each SCD run
	Variable xi								// x position on detector from its center (mm)
	Variable xp, zp							// x,z coodinates of the pixel (mm), in beam line system
	Variable zintercept						// position on beam axis where beam hits (mm)
//	Variable x0=DimOffset(slice,0)
//	Variable dx=DimDelta(slice,0)
//	Variable dmin, dmax					// range of the overlap
	Variable start = ticks
	do										// look at each run in the range
		if (strlen(StringFromList(k,range))<1)
			break
		endif
		sprintf frameName, "root:raw:%s%04d",base,str2num(StringFromList(k,range))
		Wave frame = $frameName
		if (!WaveExists(frame))
			Abort "the input wave '"+frameName+"' does cannot be found"
		endif
		trans = NumberByKey("translate_mm",note(frame),"=")
		for (i=0;i<dimX;i+=1)
			xi = i*x2mm + xLeftmm					// horiz distance on detector from center
			pixelXZfromDetector(deta,detd,xi,xp,zp)	// calculate pixel coordinates xp,zp
			zintercept = slitz - slitx*(zp-slitz)/(xp-slitx)
			// surface at center, when trans=surface
			// depth = zintercept + trans - surface,  a bigger zintercept is deeper
			if (abs(zintercept+trans-surface+depth)<=tol)
				if (thresh<.5)
					slice[i][][] += frame[i][q][r]
				else
					slice[i][][] += (frame[i][q][r]>thresh) ? frame[i][q][r] : 0
				endif
			endif
//			dmin = max(zmm-depth-res,depth-tol)		// bottom of overlap range
//			dmax = min(zmm-depth+res,depth+tol)	// top of overlap range
//			if (dmax>dmin)
//				slice[i][][] += round((dmax-dmin)/res/2*frame[i][q][r])
//			endif
		endfor
		k += 1
	while(1)
	printf "sorting took %.2f sec\r",(ticks-start)/60.15
	if (printIt)
		DoAlert 1, "put up a plot of "+sliceName
		if (V_flag)
			ShowSliceGraph(sliceName)
		endif
	endif
	return GetWavesDataFolder(slice,2)				// return the full path
End
// this routine extracts the value of varName from noteStr if it is there, otherwise, it tries to 
// use the global the global of same name in root:Packages:SCD.  If all fail, it returns NaN.
Function valFromWaveNoteOrGlobal(noteStr,varName)
	String noteStr
	String varName

	Variable val = NumberByKey(varName,noteStr,"=")
	if (numtype(val))
		String gName = "root:Packages:SCD:"+varName
		val = NumVarOrDefault(gName,NaN)
	endif
	return val
End


// ====================================================================
// ===========================  End Sorting and Summing   =====================
// ====================================================================
// ====================================================================



// ====================================================================
// ====================================================================
// ==============================  Start Displaying   =========================
// ====================================================================

Function ShowSliceGraph(source)
	String source
	if (exists(source)!=1)
		source = ""
	elseif (DimSize($source,2)<1)
		source = ""
	endif
	if (strlen(source)<1)
		source = StringFromList(ItemsInList(WaveList("*",";","DIMS:3"))-1, WaveList("*",";","DIMS:3"))
 		Prompt source, "3D wave to show", popup, WaveList("*",";","DIMS:3")
		DoPrompt "pick wave", source
		if (V_flag)
			return 1
		endif
	endif
	if (exists(source)!=1 || DimSize($source,2)<1)
		if (strlen(source)<1)
			print "no 3-d waves here"
		else
			printf "wave '%s' is un-usable\r",source
		endif
		return 1
	endif
	Display /W=(86,104,456,430)/K=1
	TextBox/N=text1/F=0/S=3/A=LB/X=-31/Y=-12 source
	String slice = UniqueName("slice_temp",1,0)
	String binName = GetDataFolder(3)+UniqueName("slice_binning",3,0)		// variable to hold binning
	Variable nx=DimSize($source,0)
	Variable ny=DimSize($source,1)
	Variable nz=DimSize($source,2)
	Variable z0=DimOffset($source,2)
	Variable dz=DimDelta($source,2)
	Make/N=(nx,ny)/O/U/I $slice
	Note $slice note($source)
	if (SetXYtoANGER($slice))			// not an Anger detector, so copy the scales.
		SetScale/P x DimOffset($source,0),DimDelta($source,0),WaveUnits($source,0), $slice
		SetScale/P y DimOffset($source,1),DimDelta($source,1),WaveUnits($source,1), $slice
	endif
	Variable/G $binName=round(nz/12)
	NVAR binning=$binName
	binning = (nz<=25) ? 1 : binning
	AppendImage $slice
	ModifyImage $slice ctab= {*,*,Terrain,1}
	ModifyGraph margin(right)=7,gfSize=14,width={Aspect,1}
	ModifyGraph tick=2,mirror=1,minor=1,standoff=0
	ModifyGraph lblMargin(left)=63, lblLatPos(left)=0
	ModifyGraph axOffset(left)=4.86,axOffset(bottom)=-0.75
	TextBox/N=text0/F=0/S=3/A=LB/X=-32/Y=-1.5 "\\Z160"
	Slider sliceSlider,pos={9,5},size={46,199},proc=upDateSliceProc
	Slider sliceSlider,limits={z0,(nz-1)*dz+z0,dz},value= z0,live= 0,tkLblRot= 90
	Button down,pos={12,207},size={17,20},proc=PushButtonProc,title="<"
	Button up,pos={36,207},size={17,20},proc=PushButtonProc,title=">"
	SetVariable setBin,pos={3,232},size={61,15},title="bin",format="%d"
	SetVariable setBin,limits={1,nz,1},value= $binName
	SetWindow kwTopWin, hook=onWindowClose
	SetWindow kwTopWin, note="sliceWindow=1;"
	SetWindow kwTopWin, note+="source="+source+";"
	SetWindow kwTopWin, note+="sourceFull="+GetDataFolder(3)+source+";"
	SetWindow kwTopWin, note+="binName="+binName+";"
	SetWindow kwTopWin, note+="value="+num2str(z0)+";"
	SetWindow kwTopWin, note+="killVariable="+binName+";"
	SetWindow kwTopWin, note+="killWave="+GetDataFolder(3)+slice+";"
	DoWindow /T $WinName(0,1), WinName(0,1)+":"+source
	upDateSliceProc("",z0,1)
End
//Function ShowSliceGraph(source)
//	String source
//	if (exists(source)!=1)
//		source = ""
//	elseif (DimSize($source,2)<1)
//		source = ""
//	endif
//	if (strlen(source)<1)
//		source = StringFromList(ItemsInList(WaveList("*",";","DIMS:3"))-1, WaveList("*",";","DIMS:3"))
 //		Prompt source, "3D wave to show", popup, WaveList("*",";","DIMS:3")
//		DoPrompt "pick wave", source
//		if (V_flag)
//			return 1
//		endif
//	endif
//	if (exists(source)!=1 || DimSize($source,2)<1)
//		if (strlen(source)<1)
//			print "no 3-d waves here"
//		else
//			printf "wave '%s' is un-usable\r",source
//		endif
//		return 1
//	endif
//	Display /W=(86,104,456,430)/K=1
//	TextBox/N=text1/F=0/S=3/A=RB/X=-12.00/Y=-13.00 source
//	String slice = UniqueName("slice_temp",1,0)
//	String binName = GetDataFolder(3)+UniqueName("slice_binning",3,0)		// variable to hold binning
//	Variable nx=DimSize($source,0)
//	Variable ny=DimSize($source,1)
//	Variable nz=DimSize($source,2)
//	Variable z0=DimOffset($source,2)
//	Variable dz=DimDelta($source,2)
//	Make/N=(nx,ny)/O/U/I $slice
//	Note $slice note($source)
//	SetXYtoANGER($slice)
//	Variable/G $binName=round(nz/12)
//	NVAR binning=$binName
//	binning = (nz<=25) ? 1 : binning
//	AppendImage $slice
////	ModifyImage $slice ctab= {*,*,Rainbow,0}
//	ModifyImage $slice ctab= {*,*,Terrain,1}
//	ModifyGraph margin(right)=60,gfSize=14,width={Aspect,1}
//	ModifyGraph tick=2,mirror=1,minor=1,standoff=0
//	ModifyGraph axOffset(left)=-4.1,axOffset(bottom)=-0.75
//	TextBox/N=text0/F=0/S=3/A=RB/X=2.55/Y=3.64 "\\Z160"
//	Slider sliceSlider,pos={322,10},size={46,199},proc=upDateSliceProc
//	Slider sliceSlider,limits={z0,(nz-1)*dz+z0,dz},value= z0,live= 0,tkLblRot= 90
//	Button down,pos={325,212},size={17,20},proc=PushButtonProc,title="<"
//	Button up,pos={349,212},size={17,20},proc=PushButtonProc,title=">"
//	SetVariable setBin,pos={316,237},size={61,15},title="bin",format="%d"
//	SetVariable setBin,limits={1,nz,1},value= $binName
//	SetWindow kwTopWin, hook=onWindowClose
//	SetWindow kwTopWin, note="sliceWindow=1;"
//	SetWindow kwTopWin, note+="source="+source+";"
//	SetWindow kwTopWin, note+="sourceFull="+GetDataFolder(3)+source+";"
//	SetWindow kwTopWin, note+="binName="+binName+";"
//	SetWindow kwTopWin, note+="value="+num2str(z0)+";"
//	SetWindow kwTopWin, note+="killVariable="+binName+";"
//	SetWindow kwTopWin, note+="killWave="+GetDataFolder(3)+slice+";"
//	upDateSliceProc("",z0,1)
//End
Function onWindowClose(infoStr)
	String infoStr
	if (stringmatch(StringByKey("EVENT",infoStr),"kill"))
		String graphName = StringByKey("WINDOW",infoStr)
		GetWindow $graphName, note
		String noteStr = S_value
		Wave image = ImageNameToWaveRef(graphName,StringFromList(0,ImageNameList(graphName,";")))
		String wName = StringFromList(0,ImageNameList(graphName,";"))
		String binName = StringByKey("binName",noteStr,"=")
		RemoveImage $wName
		wName=StringByKey("killWave",noteStr,"=")
		KillWaves $wName
		String varName=StringByKey("killVariable",noteStr,"=")
		KillVariables $varName
	endif
	return 0
End
Function upDateSliceProc(name,value,event)
	String name			// name of this slider control
	Variable value		// value of the slider (if value==-1 increase, value=-2 decrease)
	Variable event		// bit field:   1=valueset, 2=mousedown, 4=mouse up,8=mouse moved

	if (!(event & 1))
		return 0
	endif
	String graphName = StringFromList(0,WinList("*",";","WIN:"))
	if (strlen(graphName)<1)
		DoAlert 0, "No Window for doing a slice update"
		return 1
	endif
	GetWindow $graphName, note
	String noteStr = S_value
	String str, source
	Wave image = ImageNameToWaveRef(graphName,StringFromList(0,ImageNameList(graphName,";")))
	source = StringByKey("sourceFull",noteStr,"=")
	Wave ws=$source
	if (!WaveExists(ws))
		print "wave '"+source+"' not found"
		return 1
	endif
	Variable nz=DimSize($source,2)
	Variable z0=DimOffset($source,2)
	Variable dz=DimDelta($source,2)
	Variable vmax = (nz-1)*dz+z0

	NVAR binning = $StringByKey("binName",noteStr,"=")
	binning = max(1,min(binning,nz))
	Variable upDownButton=0				// 0=normal, -1=down, +1=up
	if (stringmatch(name,"up"))
		upDownButton = 1
	elseif (stringmatch(name,"down"))
		upDownButton = -1
	endif
	if (upDownButton)
		Variable oldVal = NumberByKey("value",noteStr,"=")
		if (numtype(oldVal))
			print "old value not found in note"
			return 1
		endif
		value = oldVal + sign(upDownButton)*binning*dz
		if (value<z0 || value>vmax)
			beep
		endif
	endif
	value = max(z0,min(vmax,value))
	Slider sliceSlider, value=value

	Variable ivalue = round((value-z0)/dz)
	Variable i,im=min(nz,ivalue+binning)
	image = 0
	for (i=ivalue;i<im;i+=1)
		image += ws[p][q][i]
	endfor
	Wave delays = root:raw:delays
	Variable Ndelays = WaveExists(delays) ? numpnts(delays) : NaN
	Variable lambda = (nz==Ndelays-1) ? chan2lambda(value) : NaN

	noteStr = ReplaceNumberByKey("value",noteStr,value,"=")
	if (numtype(lambda))
		sprintf str, "\\Z16 %g %s " ,value, WaveUnits($source,2)
	else
		noteStr = ReplaceNumberByKey("lambda",noteStr,lambda,"=")
		sprintf str, "\\Z16 %.2f � " ,lambda
	endif
	SetWindow $graphName note=noteStr
	TextBox/C/N=text0 str
	if (value>=vmax || value<=z0)
		return -1
	endif
	return 0
End
Function PushButtonProc(ctrlName) : ButtonControl
	String ctrlName
	if (!cmpstr(ctrlName,"up"))
		upDateSliceProc(ctrlName,1,1)
	endif
	if (!cmpstr(ctrlName,"down"))
		upDateSliceProc(ctrlName,-1,1)
	endif
End



Function/T sumUpWavelength(wName,range,thresh)
	String wName				// name of wave to sum up
	String range				// specifies channel ranges to sum up
	Variable thresh				// only include pixels that are >= thresh, 0 includes all
//	if (strlen(range)<1)
//		range = "0-225,227-271"
//	endif

	Variable printIt=0
	if (strlen(wName)<1 || strlen(range)<1 || numtype(thresh) || thresh<0)
		String wList = WaveList("*",";","DIMS:3")
		String fldrSav= GetDataFolder(1)
		SetDataFolder root:raw:
		wList += WaveList("*",";","DIMS:3")
		SetDataFolder fldrSav
		if (strlen(range)<1)
			sprintf range, "0-%d",DimSize($wName,2)-1
		endif
		thresh = numtype(thresh) ? 0 : abs(thresh)
		Prompt wName, "name of wave to sum", popup, wList
		Prompt range, "range of channels (or wavelengths) to sum"
		Prompt thresh, "threshold for including a pixel, 0 includes all"
		DoPrompt "pick a wave", wName, range, thresh
		if (V_flag)
			return ""
		endif
		printIt = 1
		printf "  sumUpWavelength(\"%s\",\"%s\",%d)\r",wName,range,round(thresh)
	endif
	thresh = round(thresh)
	if (strsearch(range,".",0)>=0)			// if a decimal point found, assume input is wavelength, convert
		printf "converted lambda range = '%s'", range
		range = LambdaRange2Channels(range)
		printf "  to channel range '%s'\r", range
	endif
	String outName = wName+"sum"
	if (thresh>0)
		outName += "_"+num2istr(thresh)
	endif
	if (exists(wName)!=1)
		wName = "root:raw:"+wName
	endif
	Variable nx=DimSize($wName,0)
	Variable ny=DimSize($wName,1)
	Variable nz=DimSize($wName,2)

	String command, switches=makeWaveSwitches (NumberByKey("NUMTYPE", WaveInfo($wName,0)))
	sprintf command, "Make/O/N=(%d,%d)%s %s",nx,ny,switches,outName
	Execute command
	Wave out=$outName
	Wave in=$wName
	SetXYtoANGER(out)

	String noteStr = note(in)
	noteStr = ReplaceStringByKey("sumSource",noteStr,GetWavesDataFolder(in,2),"=")
	noteStr = ReplaceNumberByKey("threshSum",noteStr,thresh,"=")
	Note /K out
	Note out, noteStr
	range = expandRange(range,";")

	Duplicate /O in in_temp_
	in_temp_ = (in_temp_[p][q][r]>=thresh) ? in_temp_[p][q][r] : 0		// zero all pixels below thresh
	out = 0
	Variable k,m
	for (k=0;k<ItemsInList(range);k+=1)				// loop over the range of runs
		m = str2num(StringFromList(k,range))
		out += in_temp_[p][q][m]						// accumulate all pixels in mth time slice
	endfor
	KillWaves/Z in_temp_
	if (printIt)
		Display;AppendImage $outName
		Execute "GraphImageStyle()"
	endif
	return GetWavesDataFolder(out,2)				// return the full path
End




Function ApplyGamma(yesNo)				// apply or remove the gamma correction to the top most image
	Variable yesNo							// true means apply gamma, false means remove it
	if (numtype(yesNo))
		Prompt yesNo, "apply or remove gamma correction", popup, "Apply; Remove"
		DoPrompt "gamma Correction", yesNo
		if (V_flag)
			return NaN
		endif
		yesNo = (yesNo==1)
	endif
	String wName = StringFromList(0,ImageNameList("",";" ))
	if (strlen(wName)<1)
		DoAlert 0, "No image on top graph"
		return 1
	endif

	if (!yesNo)
		ModifyImage $wName lookup=$""		// remove gamma correction
		return 0
	else											// apply gamma correction
		Wave gamWave=$SelectString(exists("gamma")!=1,"gamma","root:gamma")
		if (!WaveExists(gamWave))
			DoAlert 0, "gamma wave cannot be found"
			return 1
		endif
		ModifyImage $wName lookup=gamWave
	endif
	return 0
End


Window Graph_gamma() : Graph
	PauseUpdate; Silent 1		// building window...
	Display /W=(571,202,966,410) gamma
	ModifyGraph gfSize=14
	ModifyGraph lSize=2
	ModifyGraph rgb=(0,3,65280)
	ModifyGraph tick=2
	ModifyGraph mirror=1
	ModifyGraph minor=1
	ModifyGraph lowTrip=0.001
	ModifyGraph standoff=0
	SetAxis left 0,1
	SetAxis bottom 0,1
EndMacro

// ====================================================================
// ===========================  End Displaying   =============================
// ====================================================================
// ====================================================================



// ====================================================================
// ====================================================================
// ===========================  Start Reading IN   ============================
// ====================================================================

Function setTranslatorDistances(base,range,start,delta)		// set the translator distance for each run
	String base								// for file name, e.g. base="CuSCD", or "stackA"
	String range							// range of run numbers to set
	Variable start							// position of first point in the range
	Variable delta							// step size for the range

	Variable printIt=0
	if (numtype(start) || numtype(delta) || strlen(range)<1 || strlen(base)<1)
		base = SelectString(strlen(base) , "stackA", base)
		range = SelectString(strlen(range) , "8523-8538", range)
		start = numtype(start) ? -13 : start
		delta = numtype(delta) ? 1 : delta
		Prompt range, "range of run numbers to use, e.g. '8523-8538'"
		Prompt start, "translator position of first run in range (mm)"
		Prompt delta, "translator step size between runs (mm)"
		Prompt base, "for wave name, e.g. base='stackA'"
		DoPrompt "set translator positions", range,base,start,delta
		if (V_flag)
			return 1
		endif
		printIt = 1
	endif
	printf "    setTranslatorDistances(\"%s\",\"%s\",%g,%g)\r",base,range,start,delta
	range = expandRange(range,";")
	if (numtype(start) || numtype(delta) || numtype(str2num(StringFromList(0,range))) || strlen(base)<1)
		printf "invalid inputs to setTranslatorDistances()\r"
		return 1
	endif

	String wName, str
	Variable i=0, j, j0, translate_mm
	j0 = str2num(StringFromList(0,range))	// j0 is run number of first run in range
	do
		j = str2num(StringFromList(i,range))
		if (numtype(j))
			break
		endif
		sprintf wName, "%s%04d",base,j
		translate_mm = start + (j-j0)*delta
		str = ReplaceNumberByKey("translate_mm",note($wName),translate_mm,"=")
		Note/K $wname
		Note $wName, str
		i += 1
	while(1)
End



Function ReadRangeSCD(range,suffix,num1,base,wNote)
	String range				// range of run numbers to load
	String suffix				// file name suffix (not extension), note suffix may be empty
	Variable num1				// run num to subtract to get index for wave name, 7409 for Cu"
	String base					// for file name, e.g. base="CuSCD"
	String wNote				// wave note with parameters

	if (strlen(range)<1 || numtype(num1) || strlen(base)<1)
		num1 = numtype(num1) ? 0 : num1
		base = SelectString(strlen(base) , "stackA", base)
		SVAR rangeList = root:Packages:SCD:rangeList
		range = SelectString(strlen(range),StringFromList(ItemsInList(rangeList)-1,rangeList),range)
		Prompt range, "range of run numbers in a depth scan, e.g. '3-10,12'",popup,rangeList+";other"
		Prompt suffix, "suffix for data file (probably 'a', 'b' or '')"
		Prompt num1, "run num to subtract to get index for wave name, 7409 for Cu"
		Prompt base, "for wave name, e.g. base='stackA'"
		DoPrompt "files to load", range, suffix, num1, base
		if (V_flag)
			return 0
		endif
		if (!cmpstr(range,"other"))
			Prompt range, "range of run numbers, e.g. '3-10,12'"
			DoPrompt "new range", range
			if (V_flag)
				return 0
			endif
			rangeList = AddListItem(range,rangeList,";",Inf)
		endif
		printf "   ReadRangeSCD(\"%s\",\"%s\",%d,\"%s\",\"%s\")\r",range, suffix, num1, base, wNote
	endif
	if (strlen(range)<1 || numtype(num1) || strlen(base)<1)
		Abort "invalid inputs to ReadRangeSCD()"
	endif

	PathInfo raw							// put full path into S_path
	if (!V_Flag)							// Path 'raw' does not exist, create it
		NewPath /M="Location of raw data files"/Z raw
		PathInfo raw						// put full path into S_path
	endif
	if (!V_Flag)							// path still does not exist, give up
		DoAlert 0, "cannot create path 'raw'"
		return  1
	endif
	if (strlen(wNote)<1)						// if no parameters passed, set them with dialongs
		wNote = ParametersInNote(NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN)
	endif
	String range0 = range
	range = expandRange(range,";")			// expand the range
	String fileName							// full path name to file to open
	Variable rnum, i=0
	do
		rnum = str2num(StringFromList(i,range))
		if (numtype(rnum))
			break
		endif
		sprintf fileName, "%s%04d%s.txt", S_path,rnum,suffix
		if (readSCDasciiData(fileName,base,rnum-num1,wNote))
			break
		endif
		i += 1
	while(1)
	SetAllANGERXYscale()					// set the x-y scales

	DoAlert 1,"Set the translator distances for these scans"
	if (V_Flag==1)
		setTranslatorDistances(base,range0,NaN,1)	// set the translator distance for each run
	endif
End
Function ReadRangeSCD_OLD(range,num1,base,wNote)
	String range				// range of run numbers to load
	Variable num1				// run num to subtract to get index for wave name, 7409 for Cu"
	String base					// for file name, e.g. base="CuSCD"
	String wNote				// wave note with parameters

	if (strlen(range)<1 || numtype(num1) || strlen(base)<1)
		num1 = numtype(num1) ? 0 : num1
		base = SelectString(strlen(base) , "stack", base)
		SVAR rangeList = root:Packages:SCD:rangeList
		range = SelectString(strlen(range),StringFromList(ItemsInList(rangeList)-1,rangeList),range)
		Prompt range, "range of run numbers in a depth scan, e.g. '3-10,12'",popup,rangeList+";other"
		Prompt num1, "run num to subtract to get index for wave name, 7409 for Cu"
		Prompt base, "for wave name, e.g. base='stack'"
		DoPrompt "files to load", range, num1, base
		if (V_flag)
			return 0
		endif
		if (!cmpstr(range,"other"))
			Prompt range, "range of run numbers, e.g. '3-10,12'"
			DoPrompt "new range", range
			if (V_flag)
				return 0
			endif
			rangeList = AddListItem(range,rangeList,";",Inf)
		endif
		printf "   ReadRangeSCD(\"%s\",%d,\"%s\",\"%s\")\r",range, num1, base, wNote
	endif
	if (strlen(range)<1 || numtype(num1) || strlen(base)<1)
		Abort "invalid inputs to ReadRangeSCD()"
	endif

	PathInfo raw							// put full path into S_path
	if (!V_Flag)							// Path 'raw' does not exist, create it
		NewPath /M="Location of raw data files"/Z raw
		PathInfo raw						// put full path into S_path
	endif
	if (!V_Flag)							// path still does not exist, give up
		DoAlert 0, "cannot create path 'raw'"
		return  1
	endif
	if (strlen(wNote)<1)						// if no parameters passed, set them with dialongs
		wNote = ParametersInNote(NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN)
	endif
	String range0 = range
	range = expandRange(range,";")			// expand the range
	String fileName							// full path name to file to open
	Variable rnum, i=0
	do
		rnum = str2num(StringFromList(i,range))
		if (numtype(rnum))
			break
		endif
		sprintf fileName, "%s%04da.txt", S_path,rnum
		if (readSCDasciiData(fileName,base,rnum-num1,wNote))
			break
		endif
		i += 1
	while(1)
	SetAllANGERXYscale()					// set the x-y scales

	DoAlert 1,"Set the translator distances for these scans"
	if (V_Flag==1)
		setTranslatorDistances(base,range0,NaN,-1)	// set the translator distance for each run
	endif
End


Static Function SetAllANGERXYscale()			// sets the x-y scale for (xSize x ySize) arrays to correct mm scaling
	String list = WaveList("*",";","")
	String wName
	NVAR xSize = root:Packages:SCD:xSize		// x dimension of detector (pixels)
	NVAR ySize = root:Packages:SCD:ySize		// y dimension of detector (pixels)
	Variable i=0
	do
		wName = StringFromList(i,list)
		if (strlen(wName)<1)
			break
		endif
		if (DimSize($wName,0)==xSize && DimSize($wName,1==ySize))
			SetXYtoANGER($wName)
		endif
		i += 1
	while(1)
End

Function SetXYtoANGER(w)				// sets the x-y scale of this wave to ANGER values (if OK)
	Wave w
	NVAR xSize = root:Packages:SCD:xSize		// x dimension of detector (pixels)
	NVAR ySize = root:Packages:SCD:ySize		// y dimension of detector (pixels)

	NVAR x2mm = root:Packages:SCD:x2mm		// distance between pixels in x direction (mm)
	NVAR y2mm = root:Packages:SCD:y2mm		// distance between pixels in y direction (mm)
	NVAR xLeftmm = root:Packages:SCD:xLeftmm	// distance of first (left most) from center (mm)
	NVAR yLowermm = root:Packages:SCD:yLowermm	// distance of first (bottom most) from center (mm)

	if (!WaveExists(w) || DimSize(w,0)!=xSize || DimSize(w,1)!=ySize)
		return 1								// not the right size, exit with error
	endif
	Variable xlo,dx, ylo,dy
//	SetScale/P x (xLeftmm-x2mm),x2mm,"mm", w	//	first channel is 1 not 0
//	SetScale/P y (yLowermm-y2mm),y2mm,"mm", w	//	first channel is 1 not 0
	SetScale/P x xLeftmm,x2mm,"mm", w				//	first channel is 1 not 0
	SetScale/P y yLowermm,y2mm,"mm", w			//	first channel is 1 not 0
	return 0
End



Static Function readSCDasciiData(fileName,base,waveNum,wNote)
	String fileName						// full path name to file to open
	String base							// base name of the wave to create, empty defaults to 'scd'
	Variable waveNum					// wave name = base+waveNum
	String wNote						// wave note with parameters

	if (numtype(waveNum) || strlen(base)<1)
		return 1
	endif
	if (exists("root:Packages:SCD:xSize")!=2 || exists("root:Packages:SCD:ySize")!=2 || exists("root:Packages:SCD:zSize")!=2)
		SCDinitPackage(1)
	endif
	NVAR xSize = root:Packages:SCD:xSize	// x dimension of detector (pixels)
	NVAR ySize = root:Packages:SCD:ySize	// y dimension of detector (pixels)
	NVAR zSize = root:Packages:SCD:zSize	// number of time slices

	Variable fid=0

	Open /R/Z=1 fid  as fileName
	if (!fid)
		print "unable to open file ",fileName
		return 1
	endif
	if (strlen(wNote)<1)						// if no parameters passed, set them with dialongs
		wNote = ParametersInNote(NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN,NaN)
	endif

	base = SelectString(strlen(base) , "scd", base)
	String wName=base+num2istr(waveNum)

	printf "about to read from file '%s'\r",fileName
	Variable c1, skip=0							// find the start of the data
	String line
	do
		FReadLine fid, line
		skip +=1
		c1 = char2num(line)
	while (c1==char2num("#") || c1>64)		// skip the first few lines starting with '#'
	skip -= 1
	Close fid

	// the lines in the file contain the following: 'time delay (祍)' followed by (xSize*ySize) pixel values
	// there is one line for each time slice (zSize of them)
	// and the final line is end time of the last slice

	Variable start = ticks
	LoadWave/J/L={0,skip,0,0,0 }/M/V={"\t "," ",0, 0 } fileName
	Variable readTimeSec = (ticks-start)/60.15

	Variable i = ItemsInList(S_waveNames)
	if (i!=1)
		Abort "the LoadWave command loaded "+num2str(i)+" waves"
	endif
	String wave0Name = StringFromList(0,S_waveNames)
	Wave wav = $wave0Name

	Variable N0,N1
	N0 = DimSize(wav,0)
	N1 = DimSize(wav,1)

	Make/N=(N0)/O delays_temp_
	delays_temp_ = wav[p][0]

	Redimension/n=(N0*N1) wav
	DeletePoints 0,N0, wav
	N1 -= 1
	Redimension/n=(N0,N1) wav

	Variable Nx,Ny,Nz, wrongSize=0					// make sure image dimensions match globals
	if (xSize*ySize == N1)
		Nx = xSize
		Ny = ySize
	else
		Nx = round(sqrt(N1))
		Ny = round(sqrt(N1))
		wrongSize = 1
	endif
	Nz = N0-1
	wrongSize = wrongSize || (Nz !=zSize) 

	if (wrongSize)
		Prompt Nx, "number of pixels along x direction,  global="+num2str(xSize)
		Prompt Ny, "number of pixels along y direction,  global="+num2str(ySize)
		Prompt Nz, "number of time slices in detector,  global="+num2str(zSize)
		DoPrompt "Size of image in file does not match globals",Nx,Ny,Nz
		if (V_flag)
			DoAlert 0, "Leaving original values as they are"
		else
			xSize = Nx
			ySize = Ny
			zSize = Nz
			printf "updating global images sizes to xSize=%d, ySize=%d, zSize=%d\r",xSize,ySize,zSize
		endif
	endif

	Make/N=(Nx,Ny,Nz)/U/W/O $wName				// unsigned 2-byte int
	Wave scd = $wName
	SetXYtoANGER(scd)
	scd = wav[r][Nx*q+p]								// this does not copy the last line that is NaN's
	Note /K scd
	Note scd wNote
	printf "reading took %g sec,   whole load of this file took %g sec\r",readTimeSec,(ticks-start)/60.15
	KillWaves/Z $wave0Name
	UpdateGlobalsFromList(wNote)						// update globals to reflect wave note values

//	check that delays matches the delays_temp_
	Variable updateDelay = (exists("delays")!=1)		// need to update (or make) delays
	if (!updateDelay)
		if (!compareWaves(delays,delays_temp_))
			DoAlert 1, "delays do not match, update?"
			updateDelay = (V_flag==1)
		endif
	endif
	if (updateDelay)
		print "updating the 'delays' wave"
		Duplicate/O delays_temp_ delays
		SetScale d 0,0,"祍", delays
	endif
	KillWaves/Z delays_temp_
End



Static Function/T ParametersInNote(x2mm,y2mm,xLeftmm,yLowermm,deta,detd,L1mm,Tzero,slitx,slitz,bx,slitWidth,surface)
	Variable x2mm			// distance between pixels in x direction (mm)
	Variable y2mm			// distance between pixels in y direction (mm)
	Variable xLeftmm		// distance of first (left most) from center (mm)
	Variable yLowermm		// distance of first (bottom most) from center (mm)
	Variable deta			// angle from incident beam  ot center of detector (clockwise deg.)
	Variable detd			// distance from sample to detector center (mm)
	Variable L1mm			// distance from sample to source (mm)
	Variable Tzero			// time offset (祍ec), actual time is (time in file) -Tzero
	Variable slitx			// x position of slit (mm)
	Variable slitz			// z position of slit (mm)
	Variable bx				// x-offset of neutron beam horizontally from correct axis (mm)
	Variable slitWidth		// width of slit (mm)
	Variable surface		// z-position of surface, when translation=0, surface is z-coordinate of the surface (mm)

	if (numtype(x2mm+y2mm+xLeftmm+yLowermm+deta+detd+L1mm+Tzero+slitx+slitz+slitWidth+surface))	// need to ask for values
		x2mm = numtype(x2mm) ? NumVarOrDefault("root:Packages:SCD:x2mm",1.52215) : x2mm
		y2mm = numtype(y2mm) ? NumVarOrDefault("root:Packages:SCD:y2mm",1.52535) : y2mm
		xLeftmm = numtype(xLeftmm) ? NumVarOrDefault("root:Packages:SCD:xLeftmm",-66.58399) : xLeftmm
		yLowermm = numtype(yLowermm) ? NumVarOrDefault("root:Packages:SCD:yLowermm",-85.97354) : yLowermm
		deta = numtype(deta) ? NumVarOrDefault("root:Packages:SCD:deta",-75) : deta
		detd = numtype(detd) ? NumVarOrDefault("root:Packages:SCD:detd",230) : detd
		L1mm = numtype(L1mm) ? NumVarOrDefault("root:Packages:SCD:L1mm",9371.1) : L1mm	// 9378
		Tzero = numtype(Tzero) ? NumVarOrDefault("root:Packages:SCD:Tzero",21.4) : Tzero
		slitx = numtype(slitx) ? NumVarOrDefault("root:Packages:SCD:slitx",-6.9) : slitx
		slitz = numtype(slitz) ? NumVarOrDefault("root:Packages:SCD:slitz",0) : slitz
		bx = numtype(bx) ? NumVarOrDefault("root:Packages:SCD:bx",0) : bx
		slitWidth = numtype(slitWidth) ? NumVarOrDefault("root:Packages:SCD:slitWidth",1) : slitWidth
		surface = numtype(surface) ? NumVarOrDefault("root:Packages:SCD:surface",0) : surface
		Prompt x2mm, "x distance between pixels (mm)"
		Prompt y2mm, "y distance between pixels (mm)"
		Prompt xLeftmm, "dist: left most pixel to center (mm)"
		Prompt yLowermm, " dist: lowest pixel to center (mm)"
		Prompt deta, "angle from incident beam to center of detector (clockwise degrees)"
		Prompt detd "sample to detector center distance (mm)"
		Prompt L1mm, "sample to source distance (mm)"
		Prompt Tzero, "time offset (祍ec), actual time is (time in file) -This"
		Prompt slitx, "x position of slit (mm)"
		Prompt slitz, "z position of slit (mm)"
		Prompt bx, "x-offset of neutron beam horizontally from correct axis (mm)"
		Prompt slitWidth, "width of slit (mm)"
		Prompt surface, "z-coord of surface (@ trans=0) (mm)"
		DoPrompt "Detector Size (mm not cm)",x2mm,x2mm,xLeftmm,yLowermm,deta,detd,L1mm,Tzero
		if (V_flag)
			DoAlert 0, "Default values have been used, if you don't like thisyou shouldn't have pushed cancel"
		endif
		DoPrompt "slit (mm not cm)",slitx,slitz,bx,slitWidth,surface
		if (V_flag)
			DoAlert 0, "Default values have been used, if you don't like thisyou shouldn't have pushed cancel"
		endif
	endif

	String list=""
	list = ReplaceNumberByKey("x2mm",list,x2mm,"=")
	list = ReplaceNumberByKey("y2mm",list,y2mm,"=")
	list = ReplaceNumberByKey("xLeftmm",list,xLeftmm,"=")
	list = ReplaceNumberByKey("yLowermm",list,yLowermm,"=")
	list = ReplaceNumberByKey("deta",list,deta,"=")
	list = ReplaceNumberByKey("detd",list,detd,"=")
	list = ReplaceNumberByKey("L1mm",list,L1mm,"=")
	list = ReplaceNumberByKey("Tzero",list,Tzero,"=")
	list = ReplaceNumberByKey("slitx",list,slitx,"=")
	list = ReplaceNumberByKey("slitz",list,slitz,"=")
	list = ReplaceNumberByKey("bx",list,bx,"=")
	list = ReplaceNumberByKey("slitWidth",list,slitWidth,"=")
	list = ReplaceNumberByKey("surface",list,surface,"=")

	return list
End



Static Function UpdateGlobalsFromList(values)
	String values				// keyed list of values for resetting global values
	String var,varList = "x2mm;y2mm;xLeftmm;yLowermm;deta;detd;L1mm;Tzero;slitx;slitz;bx;slitWidth;surface"
	String wholeVar
	Variable val,i=-1
	do
		i += 1
		var = StringFromList(i,varList)
		if (strlen(var)<1)
			break
		endif
		val = NumberByKey(var,values,"=")
		if (numtype(val))
			continue
		endif
		var = "root:Packages:SCD:"+var
		if (exists(var)!=2)
			continue
		endif
		NVAR global=$var
		global = val
	while(1)
End
//Static Function UpdateGlobalsFromList(values)
//	String values				// keyed list of values for resetting global values
//	String var,varList = "x2mm;y2mm;xLeftmm;yLowermm;deta;detd;L1mm;Tzero;slitx;slitz;slitWidth;surface"
//	Variable val,i=0
//	do
//		var = StringFromList(i,varList)
//		if (strlen(var)<1)
//			break
//		endif
//		var = "root:Packages:SCD:"+var
//		NVAR global=$(var+"root:Packages:SCD:")
//		if (!NVAR_Exists(global))
//			continue
//		endif
//		val = NumberByKey(var,values,"=")
//		if (numtype(val))
//			continue
//		endif
//		global = val
//		i += 1
//	while(1)
//End


// ====================================================================
// ====================================================================
// ===========================  End Reading IN   =============================
// ====================================================================



// ====================================================================
// ====================================================================
// ===========================  Start Panel Parameters   =======================
// ====================================================================


Function DetectorGeometryWindow()
	if (strlen(WinList("DetectorGeometry", ";","WIN:64"))>1)
		DoWindow /F DetectorGeometry
		return 1
	endif
	NewPanel/K=1/W=(23,49,459,367)
	DoWindow /C DetectorGeometry

	Variable deta_deg = NumVarOrDefault("root:Packages:SCD:deta",-90)
	Variable detd_mm = NumVarOrDefault("root:Packages:SCD:detd",250)
	Variable detWidth = NumVarOrDefault("root:Packages:SCD:xSize",100)*NumVarOrDefault("root:Packages:SCD:x2mm",1.5)

	// find current detector description
	Variable current = 1e20
	String icurrent="set me"
	SVAR DetectorList = root:Packages:SCD:DetectorList
	String dectectorNames="",desc,item
	Variable i
	do
		item = StringFromList(i,DetectorList)
		desc = StringByKey("desc", item,"=",",")
		if (strlen(desc)>0)
			dectectorNames += desc+";"
			if (abs(NumberByKey("deta",item,"=",",")-deta_deg)+abs(NumberByKey("detd",item,"=",",")-detd_mm) < current)
				current = abs(NumberByKey("deta",item,"=",",")-deta_deg)+abs(NumberByKey("detd",item,"=",",")-detd_mm)
				if (current <.01)
					icurrent = desc
				endif
			endif
		endif
		i += 1
	while(strlen(desc)>0)

	icurrent="set me"
	PopupMenu DetectorListPopup,pos={220,2},size={180,20},proc=PopMenuSetFromDetectorList,title="Presets"
	String cmd
	sprintf cmd, "PopupMenu DetectorListPopup,mode=1,popvalue=\"%s\",value=\"%s\"",icurrent,dectectorNames
	Execute cmd

	SetDrawLayer/K UserBack
	SetDrawLayer/K ProgBack
	RedrawDetectorGeometryInPanel(deta_deg,detd_mm,detWidth,1)

	SetDrawLayer ProgBack
	SetDrawEnv fsize= 16
	DrawText 10,18,"Set Detector Geometry"
	SetDrawLayer UserBack

	SetVariable detaDisp,pos={22,28},size={86,17},title="deta",fSize=12,limits={-180,0,0},value=root:Packages:SCD:deta
	SetVariable detdDisp,pos={22,52},size={86,17},title="detd",fSize=12,format="%g",limits={0,500,0},value=root:Packages:SCD:detd
	SetVariable xSizeDisp,pos={149,28},size={86,17},title="xSize",fSize=12,format="%d",limits={0,4096,0},value= root:Packages:SCD:xSize
	SetVariable ySizeDisp,pos={147,52},size={86,17},title="ySize",fSize=12,format="%d",limits={0,4096,0},value=root:Packages:SCD:ySize
	SetVariable zSizeDisp,pos={147,76},size={86,17},title="zSize",fSize=12,format="%d",limits={0,4096,0},barmisc={0,1000},value=root:Packages:SCD:zSize
	SetVariable x2mmDisp,pos={259,28},size={118,17},title="x2mm",fSize=12,format="%g",limits={-5,5,0},value=root:Packages:SCD:x2mm
	SetVariable y2mmDisp,pos={261,52},size={118,17},title="y2mm",fSize=12,format="%g",limits={-5,5,0},value=root:Packages:SCD:y2mm
	SetVariable xLeftmmDisp,pos={249,75},size={126,17},title="xLeftmm",fSize=12,format="%g",limits={-Inf,0,0},value=root:Packages:SCD:xLeftmm
	SetVariable yLowermmDisp,pos={239,99},size={136,17},title="yLowermm",fSize=12,format="%g",limits={-Inf,0,0},value=root:Packages:SCD:yLowermm
	SetVariable L1mmDisp,pos={298,217},size={118,17},title="L1mm",fSize=12,format="%g",limits={10,Inf,0},value= root:Packages:SCD:L1mm
	SetVariable TzeroDisp,pos={300,241},size={118,17},title="Tzero",fSize=12,format="%g",limits={-Inf,Inf,0},value=root:Packages:SCD:Tzero
	SetVariable slitxDisp,pos={22,262},size={118,17},title="slitx",fSize=12,format="%g",limits={-Inf,0,0},value=root:Packages:SCD:slitx
	SetVariable slitzDisp,pos={22,286},size={118,17},title="slitz",fSize=12,format="%g",limits={-Inf,Inf,0},value=root:Packages:SCD:slitz
	SetVariable bxDisp,pos={170,262},size={78,17},title="bx",fSize=12,format="%g",limits={-Inf,Inf,0},value=root:Packages:SCD:bx
	SetVariable slitWidthDisp,pos={171,286},size={103,17},title="slitWidth",fSize=12,format="%g",limits={0,Inf,0},value=root:Packages:SCD:slitWidth
	SetVariable surfaceDisp,pos={294,286},size={103,17},title="surface",fSize=12,format="%g",limits={-Inf,Inf,0},value=root:Packages:SCD:surface
	return 0
End


Function RedrawDetectorGeometryInPanel(deta,detd,detWidth,coords)
	Variable deta=-75				// detector angle (degrees)
	Variable detd=230				// detector distance (mm)
	Variable detWidth=150			// full width of detector (mm)
	Variable coords					// if true, draw coordinates

	if (strlen(WinList("DetectorGeometry", ";","WIN:64"))<2)
		return 1
	elseif (numtype(deta+detd+detWidth+coords))
		DoAlert 0,"invalid values in 'RedrawDetectorGeometryInPanel'"
		return 1
	endif
	DoWindow /F DetectorGeometry

	Variable xc=130,yc=250											// origin is at top left
	Variable dist = detd/2												// distance of detector in screen units
	Variable detHalf=detWidth/4										// half width of detector in screen units

	SetDrawLayer/K UserBack
	SetDrawLayer UserBack

	Variable cosdeta=cos(deta*PI/180), sindeta=sin(deta*PI/180)
	Variable detx = xc-dist*cosdeta
	Variable dety = yc+dist*sindeta
	SetDrawEnv arrow=1 ; DrawLine xc,yc,detx,dety					// arrow to detector center

	SetDrawEnv linethick=5											// This is the detector
	DrawLine detx+detHalf*sindeta,dety+detHalf*cosdeta,detx-detHalf*sindeta,dety-detHalf*cosdeta
	SetDrawEnv textxjust= (deta>-90) ? 2 : 0
	DrawText detx,dety,"detector"
	DrawText xc+(detx-xc)/2+4, yc+(dety-yc)/2+2,"detd="+num2str(detd)+"mm"

	// make a polygon of an arc from 0� to -deta�, radius is 1/2 detector distance
	Variable angle,rad=dist/2,N=10,i
	String arc
	sprintf arc, "DrawPoly %g,%g,1,1,{",xc-rad,yc
	for (i=0;i<=N;i+=1)
		angle = i*deta/N*PI/180
		arc += num2str(xc-rad*cos(angle))+","+num2str(yc+rad*sin(angle))
		if (i<N)
			arc += ","
		else
			arc += "}"
		endif

	endfor
	SetDrawEnv arrow=1,fillpat=0
	Execute arc
	SetDrawEnv textxjust= 2
	DrawText xc-rad*cos( deta/2*PI/180)-5,yc+rad*sin( deta/2*PI/180)+1,"deta="+num2str(deta)+"�"

	if (coords)
		SetDrawLayer/K ProgBack
		SetDrawLayer ProgBack
		SetDrawEnv linethick=2,arrow=2 ; DrawLine xc+6,yc,xc+158,yc// incident beam
		SetDrawEnv dash=3 ; DrawLine xc,yc,xc-dist,yc				// dashed line forward from sample
		DrawText xc+37,yc-3,"source to sample is L1"
		SetDrawEnv fillfgc=(1,52428,26586) ; DrawOval xc-5,yc-5,xc+5,yc+5	// green circle at the origin

//		xc +=1.5*dist  ;  yc -= 90									// draw the xyz coordinate frame at (xc,yc)
		xc +=172  ;  yc -= 90										// draw the xyz coordinate frame at (xc,yc)
		SetDrawEnv linethick=2,arrow=1
		DrawLine xc-3,yc+3,xc-3,yc+50
		SetDrawEnv linethick=2,arrow=1
		DrawLine xc-3,yc+3,xc-50,yc+3
		DrawText xc-24,yc+50,"+x"
		DrawText xc-66,yc+9,"+z"
		DrawText xc+9,yc+7,"+y"
		SetDrawEnv linethick=2,fillpat=0,fillfgc=(1,52428,26586)
		DrawOval xc-5,yc-5,xc+5,yc+5
		SetDrawEnv fillfgc=(0,0,0)
		DrawOval xc-2,yc-2,xc+2,yc+2
	endif
	return 0
End

Function PopMenuSetFromDetectorList(ctrlName,popNum,popStr) : PopupMenuControl
	String ctrlName
	Variable popNum
	String popStr
	if (cmpstr("DetectorListPopup",ctrlName))
		return 0
	endif

	SVAR DetectorList=root:Packages:SCD:DetectorList
	String item=""
	Variable detWidth,i
	i = 0
	do
		item = StringFromList(i,DetectorList)
		if (!cmpstr(popStr,StringByKey("desc", item,"=",",")))
			break
		endif
		i += 1
	while(strlen(item)>0)

	if (strlen(item)<5)
		return 1
	endif

	NVAR deta = root:Packages:SCD:deta
	NVAR detd = root:Packages:SCD:detd
	NVAR L1 = root:Packages:SCD:L1mm
	NVAR Tzero = root:Packages:SCD:Tzero
	NVAR xSize = root:Packages:SCD:xSize
	NVAR ySize = root:Packages:SCD:ySize
	NVAR zSize = root:Packages:SCD:zSize
	NVAR x2mm = root:Packages:SCD:x2mm
	NVAR y2mm = root:Packages:SCD:y2mm
	NVAR xLeft = root:Packages:SCD:xLeftmm
	NVAR yLower = root:Packages:SCD:yLowermm

	deta = NumberByKey("deta",item,"=",",")
	detd = NumberByKey("detd",item,"=",",")
	detWidth = NumberByKey("x2mm",item,"=",",")*NumberByKey("xSize",item,"=",",")
	RedrawDetectorGeometryInPanel(deta,detd,detWidth,0)

	L1  = NumberByKey("L1",item,"=",",")
	Tzero = NumberByKey("Tzero",item,"=",",")
	xSize = NumberByKey("xSize",item,"=",",")
	ySize = NumberByKey("ySize",item,"=",",")
	zSize = NumberByKey("zSize",item,"=",",")
	x2mm = NumberByKey("x2mm",item,"=",",")
	y2mm = NumberByKey("y2mm",item,"=",",")
	xLeft = NumberByKey("xLeft",item,"=",",")
	yLower = NumberByKey("yLower",item,"=",",")
End


Function/T FindMatchToDetectorList(wav,DetectorList)
	Wave wav
	String DetectorList			// SVAR DetectorList = root:Packages:SCD:DetectorList
	String noteStr = note(wav)
	Variable deta = NumberByKey("deta",noteStr,"=")
	Variable detd = NumberByKey("detd",noteStr,"=")
	Variable L1 = NumberByKey("L1mm",noteStr,"=")
	Variable Tzero = NumberByKey("Tzero",noteStr,"=")
//	Variable xSize = NumberByKey("xSize",noteStr,"=")
//	Variable ySize = NumberByKey("ySize",noteStr,"=")
//	Variable zSize = NumberByKey("zSize",noteStr,"=")
	Variable x2mm = NumberByKey("x2mm",noteStr,"=")
	Variable y2mm = NumberByKey("y2mm",noteStr,"=")
	Variable xLeft = NumberByKey("xLeftmm",noteStr,"=")
	Variable yLower = NumberByKey("yLowermm",noteStr,"=")
//	Variable slitx = NumberByKey("slitx",noteStr,"=")
//	Variable slitz = NumberByKey("slitz",noteStr,"=")
//	Variable bx = NumberByKey("bx",noteStr,"=")
//	Variable slitWidth = NumberByKey("slitWidth",noteStr,"=")

	Variable diff, i
	String item
	for (i=0;i<ItemsInList(DetectorList);i+=1)
		item = StringFromList(i,DetectorList)
		diff = 0
		diff += abs(deta-NumberByKey("deta",item,"=",","))
		diff += abs(detd-NumberByKey("detd",item,"=",","))
		diff += abs(L1-NumberByKey("L1",item,"=",","))/10
		diff += abs(Tzero-NumberByKey("Tzero",item,"=",","))*10
		diff += abs(x2mm-NumberByKey("x2mm",item,"=",","))*100
		diff += abs(y2mm-NumberByKey("y2mm",item,"=",","))*100
		diff += abs(xLeft-NumberByKey("xLeft",item,"=",","))
		diff += abs(yLower-NumberByKey("yLower",item,"=",","))
		if (!numtype(diff) && diff<.01)
			return StringByKey("desc",item,"=",",")
		endif
	endfor
	return ""
End
//Function ttt()
//	String DetectorList = StrVarOrDefault("root:Packages:SCD:DetectorList","")
//	Wave ww = root:raw:dent1a_20_05KspaceXYZ
//	printf "FindMatchToDetectorList(%s)=%s\r",NameOfWave(ww),FindMatchToDetectorList(ww,DetectorList)
//	Wave ww = root:raw:dent1b_20_05KspaceXYZ
//	printf "FindMatchToDetectorList(%s)=%s\r",NameOfWave(ww),FindMatchToDetectorList(ww,DetectorList)
//End

// ====================================================================
// ====================================================================
// ===========================  End Panel Parameters  ========================
// ====================================================================



// ====================================================================
// ====================================================================
// ===========================  Start Utility   ==============================
// ====================================================================

Static Function/T makeWaveSwitches (type)		// get make switches from the NUMTYPE from a WaveInfo command
	Variable type
	switch(type)
		case 8:
			return "/W"
		case 16:
		case 32:
			return "/I"
		case 72:			// 8+64
			return "/U/W"
		case 80:			//16+64
		case 96:			//32+64
			return "/U/I"
		case 4:
			return "/D"
		default:
			return ""
	endswitch
End



Function/T LambdaRange2Channels(Lrange)		// convert a range in wavelength to a range in channels
	String Lrange		// = ".7-.75,.81-.98,1.12-1.27,1.65-1.93"
	String Crange=""	// channel range = "22-27, 32-45,54-62,80-90"
	String char, num=""
	Variable dashOn=0						// flags whether to use floor() or ceil()
	Variable nchar, i
	for (i=0;i<strlen(Lrange);i+=1)
		char = Lrange[i]
		nchar=char2num(char)
		if(nchar>=48 && nchar<=57 || nchar==46)	// a digit, or '.'
			num += char					// part of a number, accumulate
		else
			if (strlen(num)>0)				// done accumulating number, so include it
				if (dashOn)
					Crange += num2str(ceil(lambda2chan(str2num(num))))
				else
					Crange += num2str(floor(lambda2chan(str2num(num))))
				endif
				dashOn = (nchar==45)
				num = ""
			endif
			Crange += char
		endif
	endfor
	if (strlen(num)>0)						// include any trailing numbers
		Crange += num2istr(lambda2chan(str2num(num)))
		num = ""
	endif
	return Crange
End



Function chan2lambda(index)		// convert time channel to wavelength (�)
	Variable index					// index into the time part of the data
	Variable printIt=0
	NVAR detd = root:Packages:SCD:detd						// distance from sample to detector center (mm)
	NVAR L1mm = root:Packages:SCD:L1mm					// distance from sample to source (mm)
	NVAR Tzero = root:Packages:SCD:Tzero					// time offset (祍ec), actual time is (time in file) -Tzero
	NVAR zSize = root:Packages:SCD:zSize					// number of time slices
	if (numtype(index))
		index = 0
		Prompt index, "channel number (range 0-"+num2istr(zSize)+")"
		DoPrompt "pick a channel", index
		if (V_flag)
			return NaN
		endif
		printIt = 1
	endif

	Variable distance = (L1mm + detd)/1000	// total path length of neutrons (m)
	Variable t0=Tzero							// delay in 祍ec
	Variable lam1, lam2, lam					// wavelengths at either end of the channel
	Wave delays = delays						// prefer local values
	if (!WaveExists(delays))
		Wave delays = root:raw:delays
	endif
	if (!WaveExists(delays))
		Abort "cannot find the 'delays' wave in chan2lambda("+num2str(index)+")"
	endif
	if (numtype(index))
		return NaN
	endif
	index = max(index,0)
	index = min(index,DimSize(delays,0)-2)
	lam1 = velocity2lambda(distance / (delays[index]+t0) * 1e6)
	lam2 = velocity2lambda(distance / (delays[index+1] +t0)* 1e6)
	lam = (lam1+lam2)/2
	if (printIt)
		printf "for channel %d, the wavelength is %g (�)\r", index, lam
	endif
	return lam
End

Function lambda2chan(lambda)			// convert wavelength (�) to time channel, returns NaN on error
	Variable lambda						// wavelength (�)
	Variable printIt=0
	if (numtype(lambda) || lambda<0)
		lambda = 1.8
		Prompt lambda, "wavelength (�)"
		DoPrompt "pick a wavelength", lambda
		if (V_flag)
			return NaN
		endif
		printIt = 1
	endif
	NVAR detd = root:Packages:SCD:detd						// distance from sample to detector center (mm)
	NVAR L1mm = root:Packages:SCD:L1mm					// distance from sample to source (mm)
	NVAR Tzero = root:Packages:SCD:Tzero					// time offset (祍ec), actual time is (time in file) -Tzero

//	Variable delay = distance / lambda2velocity(lambda)	// neutron velocity from the wavelength
	Variable distance = (L1mm + detd)/1000				// path length of neutrons (m)
	Variable delay = distance/lambda2velocity(lambda)		// neutron velocity from the wavelength
	delay = delay*1e6 - Tzero								// true delay (祍ec)
	Wave delays = delays									// prefer local values
	if (!WaveExists(delays))
		Wave delays = root:raw:delays
	endif
	Variable index = BinarySearchInterp(delays,delay)// index into the time part of the data
	if (!printIt)
		return index										// just return value, no matter what
	elseif (delay<delays[0])
			print "wavelength is too small"
			DoAlert 0, "wavelength is too small"
	elseif (delay>delays[inf])
			print "wavelength is too long"
			DoAlert 0, "wavelength is too long"
	elseif (numtype(index))
			print "something is very wrong computing wavelength"
			Abort "something is very wrong computing wavelength"
	elseif (printIt)
		printf "for a wavelength of %g (�), the nearest channel is %.1f\r", lambda,index
	endif
	return index
End



Static Function/T WaveList_NxN(n0,n1)
	Variable n0,n1
	String wList=WaveList("*",";","DIMS:2"), NxNList="", item
	Variable i
	for (i=0;i<ItemsInList(wList);i+=1)
		item = StringFromList(i,wList)
		if (DimSize($item,0)==n0 && DimSize($item,1)==n1)
			NxNList = AddListItem(item,NxNList)
		endif
	endfor
	return NxNList
End



Function isLowestOrderFCC(hkl)				// returns true if hkl is lowest order allowed FCC
	Wave hkl									// 220 is true, 044 is false
	if (!allowFCC(hkl))
		return 0
	endif

	// to have gotten this far, it is allowed FCC
	Duplicate/O hkl h_temp_
	Variable factor = RemoveCommonFactors(h_temp_)
	if (factor<1.1)								// a lowest order all odd type
		KillWaves/Z h_temp_
		return 1
	endif
	if (!allowFCC(h_temp_) && factor<2.1)	// a lowest order all even type
		KillWaves/Z h_temp_
		return 1
	endif
	KillWaves/Z h_temp_
	return 0
End

Function isLowestOrderBCC(hkl)				// returns true if hkl is lowest order allowed BCC
	Wave hkl
	if (!allowBCC(hkl))
		return 0
	endif

	// to have gotten this far, it is allowed BCC
	Duplicate/O hkl h_temp_
	Variable factor = RemoveCommonFactors(h_temp_)
	if (factor>2.1)								// common factors of 3 or greater, definitely not lowest order
		KillWaves/Z h_temp_
		return 0
	elseif (factor <1.1)							// no common factors
		KillWaves/Z h_temp_
		return 1
	endif

	// at this point, it is allowed, and factor is 2, so check if factor of 2 is needed
	h_temp_ /= 2
	Variable i = !allowBCC(h_temp_)
	KillWaves/Z h_temp_
	return i
End


// Determine the angle between two sets of G vectors whose diffracted spots are at cursors A and B
// Note, that this is not the angle between the cursors themselves, but between the G vectors
//	Note that +z is along beam, +y is up, and +x is horizontally perp to beam (points to where you stand).
Function angleBetweenGofCursors()				// determine angle of G vector between cursor A and B
	Variable x0 = hcsr(A)						// pixel poxitions on detector (mm from center)
	Variable y0 = vcsr(A)						// 		not beamline coords
	Variable x1 = hcsr(B)
	Variable y1 = vcsr(B)

	String noteStr = note(ImageNameToWaveRef("",StringFromList(0,ImageNameList("",";"))))
	Variable deta = valFromWaveNoteOrGlobal(noteStr,"deta")
	Variable detd = valFromWaveNoteOrGlobal(noteStr,"detd")
	Variable bx = valFromWaveNoteOrGlobal(noteStr,"bx")
	Variable xp, zp								// pixel position (in horizontal plane)
	Make/O/D/N=3 kihat={0,0,1}, kfhat, G0,G1
	// kihat points downstream, from the sampe away from the source
	pixelXZfromDetector(deta,detd,x0,xp,zp)	// calculate pixel coordinates xp,zp
	kfhat = {xp-bx,y0,zp}						// vertical height (y-direction) is y0
	normalize(kfhat)
	G0 = (kfhat-kihat)							// points half way between kf and -ki

	pixelXZfromDetector(deta,detd,x1,xp,zp)	// calculate pixel coordinates xp,zp
	kfhat = {xp-bx,y1,zp}						// vertical height (z-direction) is y1
	normalize(kfhat)
	G1 = (kfhat-kihat)

	Variable angle = angleBetweenVectors(G0,G1)*180/PI
	KillWaves G0,G1, kihat,kfhat
	return angle
End
//Function angleBetweenGofCursors()				// determine angle of G vector between cursor A and B
//	Variable x0 = hcsr(A)						// pixel poxitions on detector (mm from center)
//	Variable y0 = vcsr(A)						// 		not beamline coords
//	Variable x1 = hcsr(B)
//	Variable y1 = vcsr(B)
//
//	String noteStr = note(ImageNameToWaveRef("",StringFromList(0,ImageNameList("",";"))))
//	Variable deta = valFromWaveNoteOrGlobal(noteStr,"deta")
//	Variable detd = valFromWaveNoteOrGlobal(noteStr,"detd")
//	Variable xp, zp								// pixel position (in horizontal plane)
//	Make/O/D/N=3 kihat={0,0,1}, kfhat, G0,G1
//	// kihat points downstream, from the sampe away from the source
//	pixelXZfromDetector(deta,detd,x0,xp,zp)	// calculate pixel coordinates xp,zp
//	kfhat = {xp,y0,zp}							// vertical height (y-direction) is y0
//	normalize(kfhat)
//	G0 = (kfhat-kihat)							// points half way between kf and -ki
//
//	pixelXZfromDetector(deta,detd,x1,xp,zp)	// calculate pixel coordinates xp,zp
//	kfhat = {xp,y1,zp}							// vertical height (z-direction) is y1
//	normalize(kfhat)
//	G1 = (kfhat-kihat)
//
//	Variable angle = angleBetweenVectors(G0,G1)*180/PI
//	KillWaves G0,G1, kihat,kfhat
//	return angle
//End


Function pixelXZfromDetector(deta,detd,xDet,xp,zp)
	Variable deta			// angle from exit beam to detector center (measured counter clock wise) (deg)
	Variable detd			// distance from sample to detector center (mm)
	Variable xDet			// horizontal distance along the detector from detector center to pixel (mm)
	Variable &xp, &zp		// pixel position in beamline coordinates (in horizontal plane) (mm)

	Variable cosDeta = cos(deta*PI/180)
	Variable sinDeta = sin(deta*PI/180)
	xp = detd*sinDeta - xDet*cosDeta			// pixel position (in horizontal plane)
	zp = detd*cosDeta + xDet*sinDeta
	return NaN
End


Function RemoveCommonFactors(vec)			// if elements of vec are integres, remove all common factors
	Wave vec
	Variable N=numpnts(vec)
	WaveStats/Q vec
	V_max = abs(V_max)
	if (V_max>2000)
		Abort "elements of vector way too big for factoring"
	endif

	Variable i
	for (i=V_max; i>1; i-=1)
		if (mod(vec[0],i)==0 && mod(vec[1],i)==0 && mod(vec[2],i)==0)
			vec /= i
		endif
	endfor

	Variable vm=V_max
	WaveStats/Q vec
	return vm/V_max
End


//Function MatrixRotate(index,angle,mat)		// rotates matrix about x,y,z for index = 0,1,2
//	Variable index
//	Variable angle		// (radians)
//	Wave mat
//
//	Variable N = DimSize(mat,0)
//	if (DimSize(mat,1)!=N && DimSize(mat,2)!=0)
//		Abort "can only rotate square 2-d matricies"
//	endif
//
//	Variable cosa = cos(angle)
//	Variable sina = sin(angle)
//
//	Duplicate/O mat rot_temp_
//	rot_temp_ = 0
//	rot_temp_[index][index] = 1
//
//	Variable i0 = mod(index+1,N)
//	Variable i1 = mod(index+2,N)
//	rot_temp_[i0][i0] = cosa
//	rot_temp_[i1][i1] = cosa
//	rot_temp_[i0][i1] = sina
//	rot_temp_[i1][i0] = -sina
//
//	MatrixMultiply mat, rot_temp_
//	Wave M_product=M_product
//	mat = M_product
//	KillWaves/Z rot_temp_
//End


// ====================================================================
// ============================  End Utility   ==============================
// ====================================================================
// ====================================================================



// ====================================================================
// ====================================================================
// ===========================  Start Initialize   =============================
// ====================================================================


Window Panel_SCD_geometry() : Panel
	PauseUpdate; Silent 1		// building window...
	if (strlen(WinList("SCDgeometry", ";","WIN:64"))>1)
		DoWindow /F SCDgeometry
		return
	endif
	NewPanel /K=1 /W=(451,52,760,198)
	DoWindow /C SCDgeometry
	SetDrawLayer UserBack
	SetDrawEnv fillfgc= (1,52428,26586)
	DrawOval 112,121,123,132
	SetDrawEnv linethick= 2,arrow= 2
	DrawLine 123,126,275,126
	DrawText 154,123,"source to sample is L1"
	SetDrawEnv linethick= 5
	DrawLine 121,7,19,62
	SetDrawEnv arrow= 1
	DrawLine 115,122,69,36
	SetDrawEnv dash= 3
	DrawLine 110,126,10,126
	DrawText 90,73,"detd"
	DrawText 30,106,"裠eta"
	SetDrawEnv arrow= 1,fillpat= 0
	DrawPoly 59,127,1,1,{85,148,89,130,97,116,111,104,120,100,120,100}
	DrawText 43,25,"detector"
	SetDrawEnv linethick= 2,arrow= 1
	DrawLine 233,40,233,87
	SetDrawEnv linethick= 2,arrow= 1
	DrawLine 234,40,185,40
	DrawText 213,87,"+x"
	DrawText 171,46,"+z"
	DrawText 246,44,"+y"
	SetDrawEnv linethick= 2,fillpat= 0,fillfgc= (1,52428,26586)
	DrawOval 232,31,243,42
	SetDrawEnv fillfgc= (0,0,0)
	DrawOval 235,35,240,39
EndMacro



Function SCDinitPackage(doDialog)
	Variable doDialog					// if true, do the dialog even if all values look OK
	// this routine sets up global variables that describe the size and position of the detector
	if (exists("NeutronDataInitPackage")!=6 || exists("root:Packages:Neutron:Isotope")!=1 || exists("root:Packages:Neutron:Abs_xs")!=1)
		Execute/P "INSERTINCLUDE  \"Neutron\", version>=1.0"
		Execute/P "COMPILEPROCEDURES"
		NeutronDataInitPackage()
	endif

	// fwhm is only used in the reolution stuff
	//	Execute "FWHM",Execute/P "INSERTINCLUDE  \"FWHM\"";Execute/P "COMPILEPROCEDURES"

	NewDataFolder /O root:Packages
	NewDataFolder /O root:Packages:SCD

	//	Fri Sep 12, 2003  10:01:35  AM
	//	Jon,
	//	Here are the current calibration parameters.
	//
	//	DN     DETA     DETD     L1        Tzero        x2cm         y2cm         xLeft  		   yLower		     Description
	//	17  -75.00   23.00  937.11  21.4 	0.152215  0.152535  -6.658399  -8.597354	5-SCDcal 2003-08-05
	//	19 -120.00 18.00   937.11  21.4 	0.145137  0.144757  -6.551019  -8.185886	5-SCDcal 2003-08-05
	//
	//	Where DETA the detector angle in degrees, DETD is the sample-to-detector
	//	distance in centimeters, L1 is the sample-to-moderator distance in
	//	centimeters, x2cm and y2cm are channel widths in cm, xLeft and yLower
	//	are the distance of the left and lower edges from the "center" of the
	//	detector, which has coordinates of 0.0,0.0 in centimeter units.
	//
	//   the coordinate of the [0,0] pixel is (xLeft,yLower)
	//
	//	Art
	//	IPNS, Bldg. 360               Ph: 630-252-3465


		// Much   OLDER   values
		//
		// from instprm.dat
		//	phi=75,chi=90,omega=45
		//
		//	DN   DETA   DETD    L1         Tzero   x2cm     y2cm     xLeft   yLower  
		//	17 -90.00  25.00  933.68  0.0       0.1      0.1      -7.5     -7.5


	if (exists("root:Packages:SCD:DetectorList")!=2)
		String/G root:Packages:SCD:DetectorList = ""	// list of detector parameters
		SVAR DetectorList=root:Packages:SCD:DetectorList
		DetectorList = AddListItem("deta=-90,detd=320,L1=9437.1,Tzero=5.5,xSize=85,ySize=85,zSize=120,x2mm=2.9267,y2mm=3.0238,xLeft=-126.157,yLower=-124.4681,desc=first single 90� detector",DetectorList)
		DetectorList = AddListItem("deta=-120,detd=180,L1=9378,Tzero=0,xSize=100,ySize=100,zSize=272,x2mm=1.5,y2mm=1.5,xLeft=-75,yLower=-75,desc=120� nominal detector",DetectorList)
		DetectorList = AddListItem("deta=-75,detd=230,L1=9378,Tzero=0,xSize=100,ySize=100,zSize=272,x2mm=1.5,y2mm=1.5,xLeft=-75,yLower=-75,desc=75� nominal detector",DetectorList)
		DetectorList = AddListItem("deta=-120,detd=180,L1=9371.1,Tzero=21.4,xSize=100,ySize=100,zSize=272,x2mm=1.45137,y2mm=1.44757,xLeft=-65.51019,yLower=-81.85886,desc=120� 2003-08-05",DetectorList)
		DetectorList = AddListItem("deta=-75,detd=230,L1=9371.1,Tzero=21.4,xSize=100,ySize=100,zSize=272,x2mm=1.52215,y2mm=1.52535,xLeft=-66.58399,yLower=-85.97354,desc=75� 2003-08-05",DetectorList)
	endif

	if (exists("root:Packages:SCD:xSize")!=2)
		Variable/G root:Packages:SCD:xSize = NaN	// x dimension of detector (pixels)
	endif
	if (exists("root:Packages:SCD:ySize")!=2)
		Variable/G root:Packages:SCD:ySize = NaN	// y dimension of detector (pixels)
	endif
	if (exists("root:Packages:SCD:zSize")!=2)
		Variable/G root:Packages:SCD:zSize = NaN	// number of time slices
	endif
	if (exists("root:Packages:SCD:x2mm")!=2)
		Variable/G root:Packages:SCD:x2mm = NaN	// distance between pixels in x direction (mm)
	endif
	if (exists("root:Packages:SCD:y2mm")!=2)
		Variable/G root:Packages:SCD:y2mm = NaN	// distance between pixels in y direction (mm)
	endif
	if (exists("root:Packages:SCD:xLeftmm")!=2)
		Variable/G root:Packages:SCD:xLeftmm = NaN	// distance of first (left most) from center (mm)
	endif
	if (exists("root:Packages:SCD:yLowermm")!=2)
		Variable/G root:Packages:SCD:yLowermm = NaN	// distance of first (bottom most) from center (mm)
	endif

	if (exists("root:Packages:SCD:phi")!=2)
		Variable/G root:Packages:SCD:deta = NaN		// angle from incident beam to center of detector (clockwise deg.)
	endif
	if (exists("root:Packages:SCD:detd")!=2)
		Variable/G root:Packages:SCD:detd = NaN		// distance from sample to detector center (mm)
	endif
	if (exists("root:Packages:SCD:L1mm")!=2)
		Variable/G root:Packages:SCD:L1mm = NaN	// distance from sample to source (mm)
	endif
	if (exists("root:Packages:SCD:Tzero")!=2)
		Variable/G root:Packages:SCD:Tzero = NaN	// time offset (祍ec), actual time is (time in file) -Tzero
	endif

	if (exists("root:Packages:SCD:slitx")!=2)
		Variable/G root:Packages:SCD:slitx = NaN	// x position of slit (mm)
	endif
	if (exists("root:Packages:SCD:slitz")!=2)
		Variable/G root:Packages:SCD:slitz = NaN	// z position of slit (mm)
	endif
	if (exists("root:Packages:SCD:bx")!=2)
		Variable/G root:Packages:SCD:bx = NaN		// x-offset of neutron beam horizontally from correct axis (mm)
	endif
	if (exists("root:Packages:SCD:slitWidth")!=2)
		Variable/G root:Packages:SCD:slitWidth = NaN // width of slit (mm)
	endif

	if (exists("root:Packages:SCD:surface")!=2)
		Variable/G root:Packages:SCD:surface = NaN	// z-coord of surface (when trans=0) (mm)
	endif

	if (exists("root:Packages:SCD:rangeList")!=2)
		String/G root:Packages:SCD:rangeList = "8461-8491;8492-8522;8523-8538;"
	endif
	if (exists("root:Packages:SCD:materialsList")!=2)
		String/G root:Packages:SCD:materialsList = "Cu;Nb;"
	endif
	if (exists("root:Packages:SCD:aoList")!=2)
		String/G root:Packages:SCD:aoList = ""
	endif
	if (exists("root:Packages:SCD:BravaisList")!=2)
		String/G root:Packages:SCD:BravaisList = ""
	endif

	NVAR xSize = root:Packages:SCD:xSize		// x dimension of detector (pixels)
	NVAR ySize = root:Packages:SCD:ySize		// y dimension of detector (pixels)
	NVAR zSize = root:Packages:SCD:zSize		// number of time slices
	NVAR x2mm = root:Packages:SCD:x2mm		// distance between pixels in x direction (mm)
	NVAR y2mm = root:Packages:SCD:y2mm		// distance between pixels in y direction (mm)
	NVAR xLeftmm = root:Packages:SCD:xLeftmm	// distance of first (left most) from center (mm)
	NVAR yLowermm = root:Packages:SCD:yLowermm// distance of first (bottom most) from center (mm)
	NVAR deta = root:Packages:SCD:deta			// angle from incident beam ot center of detector (clockwise deg.)
	NVAR detd = root:Packages:SCD:detd			// distance from sample to detector center (mm)
	NVAR L1mm = root:Packages:SCD:L1mm		// distance from sample to source (mm)
	NVAR Tzero = root:Packages:SCD:Tzero		// time offset (祍ec), actual time is (time in file) -Tzero
	NVAR slitx = root:Packages:SCD:slitx			// x position of slit (mm)
	NVAR slitz = root:Packages:SCD:slitz			// z position of slit (mm)
	NVAR bx = root:Packages:SCD:bx				// x position of neutron beam (beam not on z-axis) (mm)
	NVAR slitWidth = root:Packages:SCD:slitWidth// width of slit (mm)
	NVAR surface = root:Packages:SCD:surface	// z-coord of surface (when trans=0) (mm)
	SVAR materialsList = root:Packages:SCD:materialsList	// list of all materials in the sample

	doDialog += numtype(xSize)|| numtype(ySize)|| numtype(zSize)
	doDialog += numtype(x2mm)|| numtype(y2mm)|| numtype(xLeftmm)|| numtype(yLowermm)
	doDialog += numtype(deta)|| numtype(detd)|| numtype(L1mm)|| numtype(Tzero)
	doDialog += numtype(slitx+slitz+slitWidth+bx)
	doDialog += (strlen(materialsList)<1)
	if (doDialog)
		xSize = numtype(xSize) ? 100 : xSize
		ySize = numtype(ySize) ? 100 : ySize
		zSize = numtype(zSize) ? 272 : zSize
		x2mm = numtype(x2mm) ? 1.52215 : x2mm
		y2mm = numtype(y2mm) ? 1.52535 : y2mm
		xLeftmm = numtype(xLeftmm) ? -66.58399 : xLeftmm
		yLowermm = numtype(yLowermm) ? -85.97354 : yLowermm
		deta = numtype(deta) ? -75 : deta
		detd = numtype(detd) ? 230 : detd
		L1mm = numtype(L1mm) ? 9371.1 : L1mm		// 9336.8
		Tzero = numtype(Tzero) ? 21.4 : Tzero
		slitx = numtype(slitx) ? -6.9 : slitx
		slitz = numtype(slitz) ? 0.00 : slitz
		bx = numtype(bx) ? 0.0 : bx
		slitWidth = numtype(slitWidth) ? 1.00 : slitWidth
		surface = numtype(surface) ? 0 : surface
		Variable x_Size=xSize, y_Size=ySize, z_Size=zSize
		Variable dxpixel_mm=x2mm, dypixel_mm=y2mm
		Variable xLeft_mm=xLeftmm, yLower_mm=yLowermm
		Variable deta_deg=deta, detd_mm=detd, L1_mm=L1mm, Tzero_usec=Tzero
		Variable slitx_mm=slitx, slitz_mm=slitz, bx_mm=bx, slitWidth_mm=slitWidth,surface_mm=surface
		String materials_List=materialsList
		Prompt x_Size, "number of pixels along x direction"
		Prompt y_Size, "number of pixels along y direction"
		Prompt z_Size, "number of time slices in detector"
		Prompt dxpixel_mm, "distance between pixels in x (mm)"
		Prompt dypixel_mm, "distance between pixels in y (mm)"
		Prompt xLeft_mm, "dist: left most pixel to center (mm)"
		Prompt yLower_mm, " dist: lowest pixel to center (mm)"
		Prompt deta_deg, "angle from incident beam to center of detector (clockwise degrees)"
		Prompt detd_mm, "distance from sample to detector center (mm)"
		Prompt L1_mm, "distance from sample to source (mm)"
		Prompt Tzero_usec, "time offset (祍ec), actual time is (time in file) -This"
		Prompt slitx_mm, "x position of slit (mm)"
		Prompt slitz_mm, "z position of slit (mm)"
		Prompt bx_mm, "offset of neutron beam in x-direcion from correct axis (mm)"
		Prompt slitWidth_mm, "width of slit (mm)"
		Prompt surface_mm, "z-coord of surface (@ trans=0) (mm)"
		Prompt materials_List, "list of all materials, semi-colon separated"
		DoPrompt "Detector Size (mm not cm)",x_Size,y_Size,z_Size,dxpixel_mm,dypixel_mm,xLeft_mm,yLower_mm
		if (V_flag)
			DoAlert 0, "Default values have been set, if you don't like this, rerun 'SCDinitPackage(1)'"
		endif
		xSize = x_Size
		ySize = y_Size
		zSize = z_Size
		x2mm = dxpixel_mm
		y2mm = dypixel_mm
		xLeftmm = xLeft_mm
		yLowermm = yLower_mm
		if (V_flag)
			return 1
		endif
		DoPrompt "Detector Placement (mm not cm)",deta_deg,detd_mm,L1_mm,Tzero_usec
		if (V_flag)
			DoAlert 0, "Default values have been set, if you don't like this, rerun 'SCDinitPackage(1)'"
		endif
		deta = deta_deg
		detd = detd_mm
		L1mm = L1_mm
		Tzero = Tzero_usec
		DoPrompt "Slit Placement (mm not cm)",slitx_mm,slitz_mm,bx_mm,slitWidth_mm,surface_mm
		if (V_flag)
			DoAlert 0, "Default values have been set, if you don't like this, rerun 'SCDinitPackage(1)'"
		endif
		slitx = slitx_mm
		slitz = slitz_mm
		bx = bx_mm
		slitWidth = slitWidth_mm
		surface = surface_mm
		DoPrompt "materials",materials_List
		if (V_flag)
			DoAlert 0, "Default values have been set, if you don't like this, rerun 'SCDinitPackage(1)'"
		endif
		materialsList = materials_List
	endif

	// now make sure that lattice constants and types are set properly
	SVAR aoList = root:Packages:SCD:aoList				// lattice constants list of all materials in materialsList
	SVAR BravaisList = root:Packages:SCD:BravaisList	// lattice type list of all materials in materialsList
	Variable i, N=ItemsInList(materialsList)
	if (!doDialog && N==ItemsInList(aoList) && N==ItemsInList(BravaisList))		// should be OK
		return 0
	endif
	aoList=""
	BravaisList=""
	for (i=0;i<N;i+=1)
		strswitch(StringFromList(i,materialsList))
			case "Cu":				// for Copper
				aoList = AddListItem("3.615", aoList,";",Inf)
				BravaisList = AddListItem("FCC", BravaisList,";",Inf)
				break
			case "Nb":				// for Niobium
				aoList = AddListItem("3.306", aoList,";",Inf)
				BravaisList = AddListItem("BCC", BravaisList,";",Inf)
				break
			case "Si":				// for Silicon
				aoList = AddListItem("5.43102088", aoList,";",Inf)
				BravaisList = AddListItem("FCC", BravaisList,";",Inf)
				break
			case "Al":				// for Aluminum
				aoList = AddListItem("4.05", aoList,";",Inf)
				BravaisList = AddListItem("FCC", BravaisList,";",Inf)
				break
			default:
				DoAlert 0, "Do not know about lattice of "+StringFromList(i,materialsList)
		endswitch
	endfor
End