#pragma rtGlobals=1		// Use modern global access method.
#pragma ModuleName=microGeo
#pragma version = 2.49


// version 2.4, changed pixel2XYZ() to return position in microns (not mm)
// version 2.45, added ability to set geometry from web server, and the wire axis

// with 2.1, I can  reverse the sense of pixel X using the constant CCD_X_REVERSE
Constant CCD_X_REVERSE = 0								// 0 is no reverse, 1 is reverse
//Constant USE_DISTORTION_DEFAULT = 0					// default is to NOT USE distortion
Constant USE_DISTORTION_DEFAULT = 1					// default is TO USE distortion
//Constant cosThetaWire=0.707106781186547, sinThetaWire=0.707106781186548		//  angle wire moves, usually 45 degrees
Strconstant EPICS_PREFIX="yum:geometry:"				// prefix for the geometry PVs

Menu "Help"
	"Micro Diffraction",DisplayHelpTopic/K=1 "3D-Xray Diffraction"
End
Menu "micro", dynamic
	"-"
	"ÐÐ> micro Panel", /Q, MakeMicroPanel(-1)
	"-"
	"Set Geometry Parameters...",MakeGeometryParametersPanel("")
	help={"Show panel used to set (or view) all of the geometry parameters"}
	Submenu DistortionCorrectionMenuItem(-1)
		DistortionCorrectionMenuItem(0) ,/Q, setUseDistortion(0)		// 0 is OFF
		DistortionCorrectionMenuItem(1) ,/Q, setUseDistortion(1)		// 1 is ON
	End
End
//
Function/S DistortionCorrectionMenuItem(m)				// note, "!\022(" causes a menu item to be checked and disabled
	Variable m
	Init_microGeo()
	NVAR useDistortion=root:Packages:geometry:useDistortion
	if (m==-1)											// -1 is special for the menu header
		return "Distortion Correction is "+SelectString(useDistortion,"OFF","ON")
	elseif (m %^ useDistortion)								// only one is true
		return SelectString(useDistortion,"Start","Stop")+" Using Distortion Correction"
	else
		return "!\022("+SelectString(useDistortion,"Not ","")+"Using Distortion Correction"
	endif
End


Function MakeMicroPanel(tab)						// makes the main microPanel
	Variable tab
	if (WinType("micropanel")==7)
		DoWindow/F microPanel
		if (!(tab>=0 && tab<=4))
			ControlInfo/W=microPanel tabMicro
			tab = V_Value
		endif
	else
		NewPanel /W=(672,60,991,683)/K=1/N=microPanel
		TabControl tabMicro,pos={0,0},size={317,18},proc=microGeneralTabProc
		TabControl tabMicro,help={"The main micro-beam functions"},userdata(tabNum)="3"
//		TabControl tabMicro,tabLabel(0)="Index",tabLabel(1)="E-W scan"
		TabControl tabMicro,tabLabel(0)="Index",tabLabel(1)="E scans"
		TabControl tabMicro,tabLabel(2)="3D-arrays",tabLabel(3)="Geo"
		TabControl tabMicro,tabLabel(4)="Xtal"
		MoveWindowToCorner("microPanel","tr")	// re-position to the top right corner
	endif
	tab = (tab>=0 && tab<=4) ? tab : 3				// default to geometry
	TabControl tabMicro,value=tab					// set control to tab
	STRUCT WMTabControlAction  tca				// open panel to the geometry tab
	tca.ctrlName="tabMicro"
	tca.win="microPanel"
	tca.eventCode=2
	tca.tab = tab
	tca.ctrlRect.bottom = 18
	microGeneralTabProc(tca)
End
//
Function microGeneralTabProc(tca) : TabControl		// changes the panel for each tab when a tab is selected
	STRUCT WMTabControlAction &tca
	switch( tca.eventCode )
		case -1:									// control killed
			return 1
		case 2:										// mouse up
			Variable tab = tca.tab
			break
	endswitch
	if (NumberByKey("tabNum",GetUserData(tca.win,"tabMicro","tabNum"),"=")==tab)
		return 0
	endif
	TabControl $tca.ctrlName, win=$tca.win,userdata(tabNum)=num2istr(tab)
	String win = GetUserData(tca.win,"tabMicro","panelName")
	if (WinType(win)==7 && strlen(win))
		KillWindow $win
	endif

	String fillFunctionList = "Index;EWscan;Arrays3d;Geometry;Lattice"
	String funcName = "Fill"+StringFromList(tab,fillFunctionList)+"ParametersPanel"
	if (exists(funcName)!=6)
		funcName = "Load"+funcName[4,Inf]
	endif
	if (exists(funcName)==6)
		FUNCREF FillGeometryParametersPanel FillParametersPanel = $funcName
		win = tca.win+FillParametersPanel("","microPanel",45,2+tca.ctrlRect.bottom)
	else
		win = ""
	endif
	TabControl $tca.ctrlName, win=$tca.win,userdata(panelName)=win
	SetActiveSubwindow microPanel
	return 0
End
//
Function/T LoadLatticeParametersPanel(strStruct,hostWin,left,top)
	String strStruct									// optional passed value of xtal structure, this is used if passed
	String hostWin									// name of home window
	Variable left, top									// offsets from the left and top

	SetWindow kwTopWin,userdata(LatticePanelName)=hostWin+"#LatticePanel"
	NewPanel/K=1/W=(left,top,left+221,top+365)/HOST=$hostWin
	ModifyPanel frameStyle=0, frameInset=0
	RenameWindow #,LatticePanel
	Button buttonInitLattice,pos={33,31},size={150,50},title="Init Lattice\rpackage",proc=LoadPackageButtonProc
	Button buttonInitLattice,help={"Load lattice symmetry procedures"}
	return "#LatticePanel"
End
//
Function/T LoadIndexParametersPanel(strStruct,hostWin,left,top)
	String strStruct									// optional passed value of xtal structure, this is used if passed
	String hostWin									// name of home window
	Variable left, top									// offsets from the left and top

	SetWindow kwTopWin,userdata(IndexPanelName)=hostWin+"#IndexPanel"
	NewPanel/K=1/W=(left,top,left+221,top+365)/HOST=$hostWin
	ModifyPanel frameStyle=0, frameInset=0
	RenameWindow #,IndexPanel
	Button buttonInitIndex,pos={33,31},size={150,50},title="Init Indexing\rpackage",proc=LoadPackageButtonProc
	Button buttonInitIndex,help={"Load Package for Fitting Peaks and Indexing patterns"}
	return "#IndexPanel"
End
//
Function/T LoadEWscanParametersPanel(strStruct,hostWin,left,top)
	String strStruct									// optional passed value of xtal structure, this is used if passed
	String hostWin									// name of home window
	Variable left, top									// offsets from the left and top

	SetWindow kwTopWin,userdata(EWscanPanelName)=hostWin+"#EWscanPanel"
	NewPanel/K=1/W=(left,top,left+221,top+365)/HOST=$hostWin
	ModifyPanel frameStyle=0, frameInset=0
	RenameWindow #,EwscanPanel
	Button buttonInitEWscan,pos={33,31},size={150,50},title="Init E-W scan\rpackage",proc=LoadPackageButtonProc
	Button buttonInitEWscan,help={"Load Package for processing E-W scans"}
	return "#EWscanPanel"
End
//Function/T LoadEWscanParametersPanel(strStruct,hostWin,left,top)
//	String strStruct									// optional passed value of xtal structure, this is used if passed
//	String hostWin									// name of home window
//	Variable left, top									// offsets from the left and top
//
//	SetWindow kwTopWin,userdata(EWscanPanelName)=hostWin+"#EWscanPanel"
//	NewPanel/K=1/W=(left,top,left+221,top+365)/HOST=$hostWin
//	ModifyPanel frameStyle=0, frameInset=0
//	RenameWindow #,EwscanPanel
//	if (strlen(WinList("EnergyWireScans.*", "", "WIN:128")))
//		TitleBox titleLoadedEnergyWire,pos={18,37},size={180,40},title="\\JCmacros already loaded,\rsee the menubar"
//		TitleBox titleLoadedEnergyWire,fSize=16,frame=0
//	else
//		Button buttonInitEWscan,pos={33,31},size={150,50},title="Init E-W scan\rpackage",proc=LoadPackageButtonProc
//		Button buttonInitEWscan,help={"Load Package for processing E-W scans"}
//	endif
//	return "#EWscanPanel"
//End
//
Function/T LoadArrays3dParametersPanel(strStruct,hostWin,left,top)
	String strStruct									// optional passed value of xtal structure, this is used if passed
	String hostWin									// name of home window
	Variable left, top									// offsets from the left and top

	SetWindow kwTopWin,userdata(Arrays3dPanelName)=hostWin+"#Arrays3dPanel"
	NewPanel/K=1/W=(left,top,left+221,top+365)/HOST=$hostWin
	ModifyPanel frameStyle=0, frameInset=0
	RenameWindow #,Arrays3dPanel
	if (strlen(WinList("ArrayOf3dOrients.*", "", "WIN:128")))
		TitleBox titleLoadedArray3d,pos={18,37},size={180,40},title="\\JCmacros already loaded,\rsee the menubar"
		TitleBox titleLoadedArray3d,fSize=16,frame=0
	else
		Button buttonInitArrays3d,pos={33,31},size={150,50},title="Init 3d Arrays\rpackage",proc=LoadPackageButtonProc
		Button buttonInitArrays3d,help={"Processes 3D array of matricies into viewable of rotataion angles, dislocation tensors, GND, etc."}
	endif
	return "#Arrays3dPanel"
End
//
Function LoadPackageButtonProc(ctrlName) : ButtonControl
	String ctrlName

	if (stringmatch(ctrlName,"buttonInitLattice"))
		Execute/P "INSERTINCLUDE  \"LatticeSym\", version>=3.32";Execute/P "COMPILEPROCEDURES ";Execute/P "InitLatticeSymPackage()"
	elseif (stringmatch(ctrlName,"buttonInitIndex"))
		Execute/P "INSERTINCLUDE  \"Indexing\", version>=2.21";Execute/P "COMPILEPROCEDURES ";Execute/P "InitIndexingPackage()"
	elseif (stringmatch(ctrlName,"buttonInitIndexLots"))
		Execute/P "INSERTINCLUDE  \"IndexLots\", version>=2.06";Execute/P "COMPILEPROCEDURES ";Execute/P "initSymmetryOperations()"
	elseif (stringmatch(ctrlName,"buttonInitEWscan"))
		Execute/P "INSERTINCLUDE  \"EnergyWireScans\", version>=0.981";Execute/P "COMPILEPROCEDURES "
	elseif (stringmatch(ctrlName,"buttonInitArrays3d"))
		Execute/P "INSERTINCLUDE  \"ArrayOf3dOrients\", version>=1.8";Execute/P "COMPILEPROCEDURES "
	else
		return 1
	endif

	TabControl tabMicro,value=3				// set control to tab
	STRUCT WMTabControlAction  tca			// open panel to the geometry tab
	tca.ctrlName="tabMicro"
	tca.win="microPanel"
	tca.eventCode=2
	tca.ctrlRect.bottom = 18
	tca.tab = 3
	microGeneralTabProc(tca)
	return 0
End
//
Static Function MoveWindowToCorner(win,corner)
	String win
	String corner			//  one of "TL", "TR", "BL", "BR"
	if (strlen(win)<1)
		win = StringFromList(0,WinList("*",";","WIN:71"))	// graphs, tables, layouts, panels
	endif
	if (strlen(win)<1)
		return 1
	endif
	String list = StringByKey("SCREEN1",IgorInfo(0))
	Variable left,top,right,bottom, dx,dy, i=strsearch(list,"RECT=",0)
	sscanf list[i+5,Inf], "%g,%g,%g,%g", left,top,right,bottom
	GetWindow $win,wsize
	if (stringmatch(corner,"TL"))
		dx = left - V_left
		dy = (top-V_top)+44
	elseif (stringmatch(corner,"TR"))
		dx = right-V_right
		dy = (top-V_top)+44
	elseif (stringmatch(corner,"BL"))
		dx = left - V_left
		dy = bottom-V_bottom
	elseif (stringmatch(corner,"BR"))
		dx = right-V_right
		dy = bottom-V_bottom
	else
		return 1												// invalid corner
	endif
	MoveWindow/W=$win V_left+dx,V_top+dy,V_right+dx,V_bottom+dy
	return 0
End














// to convert from a ROI pixel (zero based) to full unbinned chip use:
//
//	px = (startx-1) + px*groupx + (groupy-1)/2	// pixel is zero based here
//	py = (starty-1) + py*groupy + (groupy-1)/2




//  Note about the origin.	March 22, 2006
//
// Igor uses (0,0) as the first point in the image.   Wenge, uses (1,1) as the first point
// So:
// 
// changed peakcorrection2() to correct for this, nothing else.  Be Careful of meaning of geo.xcent and geo.ycent
// Wenge's program computes them with origin at (1,1), I use the origin (0,0).  So I should take his values minus 1.



// Beam line corodinate system (X=out door, Y=up, Z=downstream)
// Wenge corodinate system (X=out door, Y=upstream, Z=up)
//
// Xbeam = Xwenge
// Ybeam = Zwenge
// Zbeam = -Ywenge
//
// Xwenge = Xbeam
// Ywenge = -Zbeam
// Zwenge = Ybeam
//
Function YZ2F(y,z)		// F = -Y*cos(angle) + Z*sin(angle)
	Variable Y,Z
	Variable cosTheta = NumVarOrDefault("root:Packages:geometry:cosThetaWire",cos(PI/4))
	Variable sinTheta = NumVarOrDefault("root:Packages:geometry:sinThetaWire",sin(PI/4))
	return -Y*cosTheta + Z*sinTheta
End
Function YZ2H(Y,Z)		// H =  Y*sin(angle) + Z*cos(angle)
	Variable Y,Z
	Variable cosTheta = NumVarOrDefault("root:Packages:geometry:cosThetaWire",cos(PI/4))
	Variable sinTheta = NumVarOrDefault("root:Packages:geometry:sinThetaWire",sin(PI/4))
	return Y*sinTheta + Z*cosTheta
End
Function HF2Y(H,F)	// Y =  H*sin(angle) - F*cos(angle)
	Variable H,F
	Variable cosTheta = NumVarOrDefault("root:Packages:geometry:cosThetaWire",cos(PI/4))
	Variable sinTheta = NumVarOrDefault("root:Packages:geometry:sinThetaWire",sin(PI/4))
	return H*sinTheta - F*cosTheta
End
Function HF2Z(H,F)	// Z =  H*cos(angle) + F*sin(angle)
	Variable H,F
	Variable cosTheta = NumVarOrDefault("root:Packages:geometry:cosThetaWire",cos(PI/4))
	Variable sinTheta = NumVarOrDefault("root:Packages:geometry:sinThetaWire",sin(PI/4))
	return H*cosTheta + F*sinTheta
End




Structure microGeometry	// structure definition
	int32 NxCCD, NyCCD		// # of un-binned pixels in CCD, also called xdim, ydim
	double dpsx, dpsy			// size of CCD (mm)
	double xcent, ycent		// center of CCD in un-binned 1 based pixels directly over the Si calibration std.
	double ddOffset			// dd = CCDy + ddOffset, CCDy is CCD height reported by Heidenhain (µm)
	double dd					// height of CCD (µm)
	double xbet, xgam			// incident beam tilt (degrees)
	double xalfd, xbetd		// CCD tilt (degrees)
	double ki[3]				// unit vector in inpt direction
	double rded00, rded01, rded02	// rotation matrix for detector
	double rded10, rded11, rded12
	double rded20, rded21, rded22
	double rde00, rde01, rde02		// rotatation matrix from Wenge coord to Ideal BeamLine system
	double rde10, rde11, rde12
	double rde20, rde21, rde22
	STRUCT wireGeometry wire
EndStructure
//
Structure wireGeometry		// structure definition
	double F					// F of the wire for the Si calibration wire scan (µm)
	double H0					// H of wire centered on the direct beam (µm)
	double Hyc					// H of wire centered on the vertical over the Si (µm)
	double X					// X of wire for Si (µm)
	double dia					// diameter of wire (µm)
	double xyzSi[3]			// wire beam line coordinates that would put wire center at the Si position, computed (µm)
	double xyz0[3]			// wire beam line coordinates when wire is on the incident beam & Fwire is F, computed (µm)
	double axis[3]				// unit vector in direction of wire axis
EndStructure
//
Function printGeometry(geo)
	STRUCT microGeometry &geo
	printf "using geomery parameters:\r"
	printf "	NxCCD=%d, NyCCD=%d		// number of un-binned pixels in CCD, also called xdim, ydim\r",geo.NxCCD,geo.NyCCD
	printf "	dpsx=%g, dpsy=%g			// size of CCD (mm)\r",geo.dpsx,geo.dpsy
	printf "	xcent=%.3f, ycent=%.3f	// center of CCD in un-binned 1 based pixels\r",geo.xcent,geo.ycent
	printf "	dd=%g						// height of CCD (µm)\r",geo.dd
	printf "	ddOffset=%g					// dd = CCDy + ddOffset (µm)\r",geo.ddOffset
	printf "	xbet=%.5f, xgam=%.5f		// incident beam tilt (degrees)\r",geo.xbet,geo.xgam
	printf "	xalfd=%.5f, xbetd=%.5f	// CCD tilt (degrees)\r",geo.xalfd,geo.xbetd
	printf "	ki = {%.4f, %.4f, %.4f}	// incident beam direction (Wenge coordinates)\r",geo.ki[0],geo.ki[1],geo.ki[2]
	if (NumVarOrDefault("root:Packages:geometry:printVerbose",0))
		printf "			{%+.4f, %+.4f, %+.4f}	// rotation matrix for detector (Wenge coordinates)\r",geo.rded00, geo.rded01, geo.rded02
		printf "	rded =	{%+.4f, %+.4f, %+.4f}\r",geo.rded10, geo.rded11, geo.rded12
		printf "			{%+.4f, %+.4f, %+.4f}\r",geo.rded20, geo.rded21, geo.rded22
		printf "			{%+.4f, %+.4f, %+.4f}	// rotation matrix for incident beam (Wenge coordinates)\r",geo.rde00, geo.rde01, geo.rde02
		printf "	rde =	{%+.4f, %+.4f, %+.4f}\r",geo.rde10, geo.rde11, geo.rde12
		printf "			{%+.4f, %+.4f, %+.4f}\r",geo.rde20, geo.rde21, geo.rde22
	endif
	printf "		for the wire:\r"
	printf "	F=%.2f							// F of the wire for the Si calibration wire scan (µm)\r",geo.wire.F
	printf "	H0=%.2f						// H of wire centered on the direct beam (µm)\r",geo.wire.H0
	printf "	Hyc=%.2f						// H of wire centered on the vertical over the Si (µm)\r",geo.wire.Hyc
	printf "	X=%g								// X of wire for Si (µm)\r",geo.wire.X
	printf "	diameter=%.2f						// diameter of wire (µm)\r",geo.wire.dia
	printf "	wire axis direction = {%.6f, %.6f, %.6f}	// wire (PM500 wire coordinates) on the incident beam (µm)\r",geo.wire.axis[0],geo.wire.axis[1],geo.wire.axis[2]
	printf "	Si = {%.2f, %.2f, %.2f}	// wire (PM500 wire coordinates) of the Si position (µm)\r",geo.wire.xyzSi[0],geo.wire.xyzSi[1],geo.wire.xyzSi[2]
	printf "	@beam = {%.2f, %.2f, %.2f}	// wire (PM500 wire coordinates) on the incident beam (µm)\r",geo.wire.xyz0[0],geo.wire.xyz0[1],geo.wire.xyz0[2]
End


// These two routines convert betweeen Wenge and Ideal Beam Line coordinate systems.  They are used to get to/from Indexing program (Euler)
Function wenge2IdealBeamline(geo,qW,qBL)
	STRUCT microGeometry &geo
	Wave qW, qBL
	Variable q0,q1,q2
	q0 = qW[0]*geo.rde00 + qW[1]*geo.rde01 + qW[2]*geo.rde02	// rotate qW from Wenge to Ideal Beam-line
	q1 = qW[0]*geo.rde10 + qW[1]*geo.rde11 + qW[2]*geo.rde12	//   Ideal Beam LIne has incident beam along {0,0,1}
	q2 = qW[0]*geo.rde20 + qW[1]*geo.rde21 + qW[2]*geo.rde22	// 	Wenge has incident near {0,-1,0}, but tilted a bit
	qBL = {q0,q1,q2}
End
//
Function IdealBeamLine2wenge(geo,qBL,qW)
	STRUCT microGeometry &geo
	Wave qBL, qW
	Variable q0,q1,q2
	q0 = qBL[0]*geo.rde00 + qBL[1]*geo.rde10 + qBL[2]*geo.rde20	// rotate qBL from Ideal Beam-line to Wenge
	q1 = qBL[0]*geo.rde01 + qBL[1]*geo.rde11 + qBL[2]*geo.rde21	//   Ideal Beam LIne has incident beam along {0,0,1}
	q2 = qBL[0]*geo.rde02 + qBL[1]*geo.rde12 + qBL[2]*geo.rde22	// 	Wenge has incident near {0,-1,0}, but tilted a bit
	qW = {q0,q1,q2}
End

Function EulerMatrix(alpha,bet,gam)	// Make the rotation matrix M_Euler from Euler angles (degree)
	Variable alpha,bet,gam				// Euler angles, passed as degrees
	alpha *= PI/180					// need radians internally
	bet *= PI/180
	gam *= PI/180
	Make/N=(3,3)/O/D M_Euler
	if (alpha<=-2*PI || alpha>=2*PI || bet<=-2*PI || bet>=2*PI || gam<=-2*PI || gam>=2*PI)
		M_Euler = NaN
		return 1
	endif
	Variable ca=cos(alpha), cb=cos(bet), sa=sin(alpha), sb=sin(bet)
	Variable cg = cos(gam), sg = sin(gam)
	M_Euler[0][0] = cg*cb*ca - sg*sa
	M_Euler[1][0] = -sg*cb*ca - cg*sa
	M_Euler[2][0] = sb*ca
	M_Euler[0][1] = cg*cb*sa + sg*ca
	M_Euler[1][1] = -sg*cb*sa + cg*ca
	M_Euler[2][1] = sb*sa
	M_Euler[0][2] = -cg*sb
	M_Euler[1][2] = sg*sb
	M_Euler[2][2] = cb
	return 0
End

Function MovingSiOrigin(geo)			// this calculates Wenge's moving origin
	// see the procedure "calctotalrange" in "ThreeDimX_RayMicroscopy.pro" for the origin of this calculation
	STRUCT microGeometry &geo
	Variable SiHstart = -7030							// this one parameter is not in the geometry strucutre

	Variable cosTheta = NumVarOrDefault("root:Packages:geometry:cosThetaWire",cos(PI/4))
	Variable sinTheta = NumVarOrDefault("root:Packages:geometry:sinThetaWire",sin(PI/4))
	Variable ddc = geo.dd*1e3							// 40220
	Variable wirepitchy = geo.dpsy/geo.NxCCD * 1000
	Variable yoff = 0									// full chip image
	Variable CCDypitch=24.0, CCDysize=24*2084		// un-binned images, yes, these are hard coded

	Variable startwirevertical, startwirehorizontal, sourcemin
	startwirevertical=(SiHstart-geo.wire.H0)*cosTheta
	startwirehorizontal=(SiHstart-geo.wire.Hyc)*sinTheta
	sourcemin=CCDYsize-(1+geo.ycent-yoff)*wirepitchy-yoff*CCDypitch-ddc*(CCDYsize-(1+geo.ycent-yoff)*wirepitchy-yoff*CCDypitch-startwirehorizontal)/(ddc-startwirevertical)
	return -sourcemin
End
//Function TestMovingSiOrigin()
//	STRUCT microGeometry geo
//	FillGeometryStructDefault(geo)
//	geo.dd = 40.22
//	geo.ycent = 923.22
//	geo.wire.H0 = -7277
//	geo.wire.Hyc = -6858.27
//	Variable depth = MovingSiOrigin(geo)
//	printf "Si is %g µm from the moving origin\r",depth
//	// using these number, the results are:     Si is 243.067 µm from the moving origin
//End



//	Function test()
//		STRUCT microGeometry geo
//		SVAR geoStr=geoStr
//		StructGet/S/B=2 geo, geoStr		// retrieve the information from geoStr (also gets all of geo.wire.xxx)
//	
//		print SelectString(USE_DISTORTION_DEFAULT,"no distortion","using distortion")
//	//	Variable px,py
//		Variable depth, edge = 1
//		Variable i, N=1000
//	
//		Variable ccdy,ccdz
//		Variable ccdy0, dccdy=27					// ccdy ranges from ccdy0±dccdy
//		Variable ccdz0,dccdz						// ccdz ranges from ccdz0±dccdz
//		Variable ccdz1,ccdz2						// ends of z range on ccd
//		Wave ccd=root:Packages:geometry:pixel2depth_ccd
//		pixel2XYZ(geo,1000,1000,ccd)			// returns xyz of pixel on CCD relative to the Si position (Wenge coordinates)
//		ccdy0 = ccd[2]								// change from Wenge to beamline
//		pixel2XYZ(geo,0,0,ccd)
//		ccdz1 = -ccd[1]								// change from wenge coord to beam line
//		pixel2XYZ(geo,2000,2000,ccd)
//		ccdz2 = -ccd[1]
//		ccdz0 = (ccdz1+ccdz2)/2					// center of allowed ccdz values
//		dccdz = abs(ccdz1-ccdz2)/2				// range of allowed ccdz values
//	
//		Wave Xw = $MakeUnique3Vector($"")
//		Xw = {.5, -6119, -800}
//	
//	 	Variable timer = startMSTimer
//		for (i=0;i<N;i+=1)
//	//		px = round(enoise(1000))+1000
//	//		py = round(enoise(1000))+1000
//	//		depth = pixel2depth(geo,px,py,Xw,edge)
//			ccdy = ccdy0+enoise(dccdy)
//			ccdz = ccdz0+enoise(dccdz)
//			depth = pixelXYZ2depth(geo,ccdy,ccdz,Xw,edge)
//		endfor
//		Variable seconds = stopMSTimer(timer)*1e-6
//		Variable total = (2084^2)/N * seconds
//		printf "execution took %g sec,     a whole image will take %s\r",seconds,Secs2Time(total,5,1)
//		KillWaves Xw
//	End


//Function xxxTemp()
//	STRUCT microGeometry geo
//	SVAR geoStr=geoStr
//	StructGet/S/B=2 geo, geoStr		// retrieve the information from geoStr (also gets all of geo.wire.xxx)
//
//	Wave ki = $MakeUnique3Vector($"")
//	ki = geo.ki[p]
//	printf "ki = {%g, %g, %g},   |ki| = %g\r",ki[0],ki[1],ki[2],norm(ki)
//	KillWaves ki
//End



// Given the three points, P=pixel, S=source, W=wire:
// The next three routines do calculations for the unknown point when the other two are known
//
Function calcWireH(geo,depth,py,pz,edge)	// calc H of wire that is tangent to line from depth on beam to pixel on ccd
	STRUCT microGeometry &geo
	Variable depth						// current Z coordinate in sample relative to the Si position, not distance from Si (µm)
	Variable py,pz						// y and z component of pixel (xyz) position in beam line coords (origin at Si)
	Variable edge						// 1=leading edge (usual),  0=trailing edge

	Variable R=(geo.wire.dia)/2, Fo=geo.wire.F		// wire radius, Fo of the wire
	Variable kiy=geo.ki[2], kiz=-geo.ki[1]			// convert from Wenge coords to beam line
	Variable sy,sz									// y,z coordinates of sample point beamline coords (origin at Si)
	sy = kiy/kiz*depth								// project along ki a distance depth,  sz = kiz*depth/kiz,  sx = kix*depth/kiz
	sz = depth

	Variable Xoy, Xoz								// point where wire intersects incident beam relative to Si position (µm)
	Xoy = geo.wire.xyz0[1] - geo.wire.xyzSi[1]		// origin at Si position
	Xoz = geo.wire.xyz0[2] - geo.wire.xyzSi[2]

	// eqn of line that wire travels along is:   y = z*tan(theta) + b
	Variable tanTheta,b		
	Variable cosTheta = NumVarOrDefault("root:Packages:geometry:cosThetaWire",cos(PI/4))
	Variable sinTheta = NumVarOrDefault("root:Packages:geometry:sinThetaWire",sin(PI/4))
	tanTheta = sinTheta/cosTheta					// = tan(theta), the slope
	b = Xoy - Xoz*tanTheta							// b comes from point that wire passes through on the incident beam

	// need eqn of line parallel to line from S to P, but moved over by distance r (=wire radius):     z = m*y + beta1
	Variable sinphi = (py-sy) / sqrt((py-sy)^2 + (pz-sz)^2)	// sin(phi), phi is angle of ray to pixel measured from z toward y
	Variable m = (pz-sz)/(py-sy)					// slope of line (it is translated, but parallel to P,S line
	Variable beta0 = sz-m*sy						// intercept for line through center of parallel to line from S to P
	Variable beta1 = (edge ? beta0-r/sinphi : beta0+r/sinphi)	// intercept of line from S to P tangent to leading (or trailing) edge

	Variable z = (m*b+beta1) / (1-m*tanTheta)	// z of wire at intersection of two lines
	z += geo.wire.xyzSi[2]							// convert from Si at origin to PM500 origin

	Variable H = (z-Fo*sinTheta)/cosTheta			// from:   Z =  H*cos(angle) + F*sin(angle)
	return H										// H of the wire, the returned value (wire origin)
End
//
//
// for a ray coming from depth, and tangent to wire, where will it hit detector, returns CCD(y,z) position as complex value
Function/C depthWire2pixelXYZ(geo,depth,Xw,xyzMin,xyzMax,edge)
	STRUCT microGeometry &geo
	Variable depth									// depth along incident beam, relative to the Si position (µm)
	Wave Xw										// wire (x,y,z) beam line coordinates in wire units (µm)
	Wave xyzMin,xyzMax							// two points on the detector (assumed to be ends) relative to Si
	Variable edge									// 1=leading edge (usual),  0=trailing edge
	Variable testing = (WhichListItem("TestdepthWire2pixelXYZ", GetRTStackInfo(0))>=0)

	Variable R = (geo.wire.dia)/2					// radius of wire (µm)
	Variable kiy=geo.ki[2], kiz=-geo.ki[1]			// components of ki[] in beam line coords (not Wenge coords)
	Variable sz = depth, sy = kiy/kiz * depth		// source position, relative to the Si position (beam line coords)

	Variable vy,vz,v								// (vy,vz) vector from source (at depth) to center of wire
	vy = (Xw[1] - geo.wire.xyzSi[1]) - sy
	vz = (Xw[2] - geo.wire.xyzSi[2]) - sz
	v = sqrt(vy*vy + vz*vz)						// length of v

	Variable phi0 = atan2(vy,vz)					// angle from positive z to tangent ray measured from z toward y
	Variable dphi =asin(R/v)						// angle between v and ray at tangent point (always positive)
	Variable phi = phi0 + (edge ? -dphi : dphi)
	Variable cotPhi = 1/tan(phi)

	// line from P to S is z = y*cot(phi) + b
	Variable b = sz - sy*cotPhi

	// line along surface of detector, line between xyzMin[] & xyzMax[],   y = m*z + beta1
	Variable m = (xyzMax[1]-xyzMin[1])/(xyzMax[2]-xyzMin[2])
	Variable beta1 = xyzMin[1] - m*xyzMin[2]

	Variable ccdy, ccdz								// y, z parts of the pixel location relative to Si (beam line coords)
	ccdy = (m*b + beta1) / (1-m*cotPhi)			// y value of intersection
	ccdz = ccdy*cotPhi + b							// from line P to S,   z = y*cot(phi) + b

	if (testing)
		Variable cosTheta = NumVarOrDefault("root:Packages:geometry:cosThetaWire",cos(PI/4))
		Variable sinTheta = NumVarOrDefault("root:Packages:geometry:sinThetaWire",sin(PI/4))
		printf "Hwire=%g, Fwire=%g,   Xw=(%g, %g, %g)\r",Xw[1]*sinTheta + Xw[2]*cosTheta, -Xw[1]*cosTheta + Xw[2]*sinTheta,Xw[0], Xw[1], Xw[2]
		printf "s = (%g, %g),   v = (%g, %g),   R=%gµm,  |wire-sample| = %gµm\r",sy,sz,vy,vz,R,v
		printf "line from source to tangent is:   z = %g*y + %g\r",cotPhi,b
		printf "line between xyzMin & xyzMax is:   y = %g*z + %g\r",m,beta1
		// check that ccdy,ccdz satisfy both equations
		printf " is z-y*cot(phi) = %g == %g=b?, Æ=%g\r",ccdz-ccdy*cotPhi,b,ccdz-ccdy*cotPhi-b
		printf " is y-m*z = %g == %g=beta?,  Æ=%g\r",ccdy-m*ccdz,beta1,ccdy-m*ccdz-beta1
	endif
	return cmplx(ccdy,ccdz)
End
//
//
// returns depth (starting point) of ray ending at pixel (px,py) that is tangent to leading edge of the wire
// the returned depth is relative to the Si position (µm)
Function pixel2depth(geo,px,py,Xw,edge)
	STRUCT microGeometry &geo
	Variable px,py						// uncorrected unbinned pixel position (pixels) (zero based)
	Wave Xw							// wire (x,y,z) beam line coordinates in wire units (µm)
	Variable edge						// 1=leading edge (usual),  0=trailing edge
	Variable depth=NaN					// computed depth relative to the Si position (µm)

	Wave ccd=root:Packages:geometry:pixel2depth_ccd

	pixel2XYZ(geo,px,py,ccd)			// returns xyz of pixel on CCD (µm) relative to the Si position (Wenge coordinates)
	Variable swap
	swap = ccd[1]							// convert xyz from Wenge to beam line coordinates
	ccd[1] = ccd[2]
	ccd[2] = -swap

	depth = pixelXYZ2depth(geo,ccd[1],ccd[2],Xw,edge)
	return depth						// relative to the Si position (µm)
End
//
// Returns depth (starting point) of ray ending at point xyz (probably a pixel location) that is tangent
// to leading edge of the wire.  The returned depth is relative to the Si position (µm)
// depth is the change in z value from the Si position, not the distance
Function pixelXYZ2depth(geo,ccdy,ccdz,Xw,edge)
	STRUCT microGeometry &geo
	Variable ccdy, ccdz					// y, z parts of the xyz location of pixel (beam line coords)(µm), origin at Si
	Wave Xw								// wire (x,y,z) beam line coordinates in wire PM500 units (µm)
	Variable edge							// 1=leading edge (usual),  0=trailing edge
	Variable depth							// computed depth relative to the Si position (µm)

	Variable kiy,kiz						// components of ki[] in beam line coords (not Wenge coords)
	kiy = geo.ki[2]
	kiz = -geo.ki[1]

	Variable Xwy,Xwz						// wire center relative to the Si position (using beam line coords)
	Xwy = Xw[1] - geo.wire.xyzSi[1]
	Xwz = Xw[2] - geo.wire.xyzSi[2]

	Variable Xoy, Xoz						// point where wire intersects incident beam relative to Si position (µm)
	Xoy = geo.wire.xyz0[1] - geo.wire.xyzSi[1]
	Xoz = geo.wire.xyz0[2] - geo.wire.xyzSi[2]

	Variable vy,vz, v
	vy = Xwy-ccdy						// (vy,vz) vector from pixel on ccd toward the wire center
	vz = Xwz-ccdz
	v = sqrt(vy*vy + vz*vz)				// length of vector (0,vy,vz) = Xw[]-P[]

	Variable R=(geo.wire.dia)/2			// wire radius
	Variable dphi = asin(R/v)			// angle between wire center and tangent point of wire (from ccd pixel)
	Variable phi0 = atan2(vz,vy)		// angle from yhat to V (V is vector from pixel to wire center)
	Variable tanphi = tan((edge ? phi0-dphi : phi0+dphi))

	// line from pixel to tangent point is:   z = y*tan(phio±dphi) + b
	Variable b = ccdz - ccdy*tanphi

	// line of the beam goes through Xo and with direction given by ki[]
	// line is:  y = z * (kiy/kiz) + beta
	Variable beta1 = Xoy - Xoz*(kiy/kiz)

	depth = (beta1*tanphi+b) / (1-tanphi*kiy/kiz)	// this is the z value of the intercept
	return depth						// this is depth (along z) relative to the Si
End
//
//
// Returns depth (starting point) of ray ending at point xyz (probably a pixel location) that is tangent
// to leading edge of the wire.  The returned depth is relative to the Si position (µm)
//	Function pixelXYZ2depth(geo,ccdy,ccdz,Xw,edge)
//		STRUCT microGeometry &geo
//		Variable ccdy, ccdz					// y, z parts of the xyz location of pixel (beam line coords)
//		Wave Xw							// wire (x,y,z) beam line coordinates in wire units (µm)
//		Variable edge						// 1=leading edge (usual),  0=trailing edge
//		Variable depth=NaN					// computed depth relative to the Si position (µm)
//	
//		Variable Xwy,Xwz					// wire center relative to the Si position (using beam line coords)
//		Xwy = Xw[1] - geo.wire.xyzSi[1]
//		Xwz = Xw[2] - geo.wire.xyzSi[2]
//	
//		ccdy -= geo.wire.xyzSi[1]			// change point on ccd to be relative to the Si position
//		ccdz -= geo.wire.xyzSi[2]
//	
//		Variable kiy,kiz						// components of ki[] in beam line coords (not Wenge coords)
//		kiy = geo.ki[2]
//		kiz = -geo.ki[1]
//	
//		Variable Xoy, Xoz					// point where wire intersects incident beam relative to Si position (µm)
//		Xoy = geo.wire.xyz0[1] - geo.wire.xyzSi[1]
//		Xoz = geo.wire.xyz0[2] - geo.wire.xyzSi[2]
//	
//		Variable r = (geo.wire.dia)/2		// radius of wire (µm)
//	
//		// want to find intersction of two lines
//		//	1) line going through Xo[3] along the ki[3] direction
//		//	2) line going through ccd[3] and tangent to wire at Xw[3] postion
//	
//		Variable d							// distance from ccd to wire center in yz plane
//		d = sqrt((Xwy-ccdy)^2 + (Xwz-ccdz)^2)	// distance from ccd to wire center in yz plane
//		Variable phi = asin(r/d)			// angle between wire center and tangent point (from ccd pixel)
//	
//		Variable vy,vz
//		vy = Xwy-ccdy						// (vy,vz) vector from pixel on ccd toward the wire center
//		vz = Xwz-ccdz
//	
//		// calculate theta, the angle of the correction vector, v[] (whose length is r)
//		// angle measured from  y^ (toward z^) and then points from ccd[] toward the tangent point
//		// atan2(vz,vy) is angle from ccd to the wire center, then correct to reach tangent point
//		// leading edge atan2()-(phi+90¡),  trailing atan2+(phi+90¡)
//		Variable theta =  atan2(vz,vy) + (edge ? -phi-PI/2 : phi+PI/2)
//		vy += r*cos(theta)					// and add the correction
//		vz += r*sin(theta)
//		// v now points from ccd[] toward the leading tangent point
//	
//		// so solve the following vector equation, line from ccd intersects incident beam.
//		//  the intersection is ccd+a*v == Xo+b*ki,   where  c[] = Xo[]-ccd[]
//		//
//		// to solve for a & b we need two equations, so use the y-z plane, giving two equations in a and b
//		//			a*vy -b*kiy = cy		// the y and z components of above vector equation
//		//			a*vz - b*kiz = cz
//		//
//		// multiply second eqn by vy/vz and subtract to eliminate the a term, and solve for b
//		//			b * [ kiz*vy/vz - kiy] = cz*vy/vz - py
//		//
//		Variable cy=Xoy-ccdy, cz=Xoz-ccdz, b
//		b = (cz*vy/vz-cy) / (kiz*vy/vz-kiy)
//		depth = Xoz+b*kiz					// from equation for the incident beam:  xyz = Xo[]+b*ki[]
//		return depth						// relative to the Si position (µm)
//	End









//	Function test_convert()
//		STRUCT microGeometry geo
//		// first set the geometry structure
//		geo.NxCCD=2084	; geo.NyCCD=2084		// # of un-binned pixels in CCD, also called xdim, ydim
//		geo.dpsx=50010	; geo.dpsy=50168		// size of CCD (µm)
//		geo.xcent=1091.53	; geo.ycent=939.941	// center of CCD in un-binned pixels
//		geo.dd=30017.3							// height of CCD (µm)
//		geo.xbet=0.31037	; geo.xgam=0.75037	// incident beam tilt (degrees)
//		geo.xalfd=0.10525	; geo.xbetd=0.171		// CCD tilt (degrees)
//		geo.wire.F = 3830
//		geo.wire.H0 = -4823
//		geo.wire.Hyc = -5000
//		geo.wire.X = 0
//		geo.wire.dia = 51
//		GeometryUpdateCalc(geo)					// this sets rest of geometry, formerly called detectortilt()
//		String/G geoStr
//		StructPut/S/B=2 geo, geoStr				// store the information into geoStr (this stores all of geo.wire too)
//		//
//		//	you can read the struct with:
//		//		SVAR geoStr=geoStr
//		//		StructGet/S/B=2 geo, geoStr		// retrieve the information from geoStr (also gets all of geo.wire.xxx)
//		//
//		printf "using parameters:\r"
//		printf "	NxCCD=%d, NyCCD=%d		// # of un-binned pixels in CCD, also called xdim, ydim\r",geo.NxCCD,geo.NyCCD
//		printf "	dpsx=%g, dpsy=%g			// size of CCD (µm)\r",geo.dpsx,geo.dpsy
//		printf "	xcent=%.2f, ycent=%g	// center of CCD in un-binned pixels\r",geo.xcent,geo.ycent
//		printf "	dd=%g						// height of CCD (µm)\r",geo.dd
//		printf "	xbet=%g, xgam=%g		// incident beam tilt (degrees)\r",geo.xbet,geo.xgam
//		printf "	xalfd=%g, xbetd=%g		// CCD tilt (degrees)\r",geo.xalfd,geo.xbetd
//		printf "		for the wire:\r"
//		printf "	F=%g							// F of the wire for the Si calibration wire scan\r",geo.wire.F
//		printf "	H0=%g							// H of wire centered on the direct beam\r",geo.wire.H0
//		printf "	Hyc=%g						// H of wire centered on the vertical over the Si\r",geo.wire.Hyc
//		printf "	X=%g								// X of wire for Si\r",geo.wire.X
//		printf "	diameter=%g						// diameter of wire \r",geo.wire.dia
//		printf "	Si = {%.2f, %.2f, %.2f}	// wire coordinates of the Si position (µm)\r",geo.wire.xyzSi[0],geo.wire.xyzSi[1],geo.wire.xyzSi[2]
//		printf "	@beam = {%.2f, %.2f, %.2f}	// wire coordinates on the incident beam (µm)\r\r",geo.wire.xyz0[0],geo.wire.xyz0[1],geo.wire.xyz0[2]
//	
//		Make/O/N=(2,2) testXY					// test values
//		testXY= {{1024,1098.4,100,200,300,1042,0},{1024,917.88,1500,1600,350,2040,2040}}
//		Make/O/D testTheta = {42.8988911219471,45.3535835733482,34.9344263723845,32.976490326739,55.5618162886343,24.0805,27.8326}
//	
//		Make/N=3/O/D qhat
//		String errStr
//		Variable/C zpix
//		Variable i,theta, px,py
//		for (i=0;i<DimSize(testXY,0);i+=1)
//			px = testXY[i][0]
//			py = testXY[i][1]
//			theta = pixel2q(geo,px,py,qhat)*180/PI
//			sprintf errStr, "\t\ttheta ERROR = %g",theta-testTheta[i]
//			errStr = SelectString(abs(theta-testTheta[i])>1e-13,"",errStr)
//			printf "theta(%g, %g) = %g¡,		q = {%.5f, %.5f, %.5f}%s\r",px,py,theta,qhat[0],qhat[1],qhat[2],errStr
//			zpix = q2pixel(geo,qhat)				// returns pixel position as a complex number cmplx(px,py)
//			if (abs(	px-real(zpix))+abs(py-imag(zpix)) > 1e-6)		// cannot use 1e-12 with the distortion inversion
//				printf "ERROR, recomputed pixel = (%g, %g)\r",real(zpix),imag(zpix) ; beep
//			endif
//		endfor
//		KillWaves/Z testXY, qhat, testTheta
//	End



//	Function tt1()		// a test for peakcorrection2() and its inverse peakUncorrect()
//	//	Wave xymap = $loadCCDCorrectionTable("")
//	//	printf "loaded distortion table into '%s'\r", GetWavesDataFolder(xymap,2)
//	
//		Wave xymap=root:Packages:geometry:xymap
//	//	Wave xymap = xymap
//	
//		if (!USE_DISTORTION_DEFAULT)
//			Abort "not using distortion, see 'Constant distort =' at top"
//		endif
//	
//	
//		Variable px0,py0,px1,py1
//		px0=200  ;  py0=200
//		px1=px0  ;  py1=py0
//		peakcorrection2(xymap,px1,py1)
//		printf "x = %g  -->  %g,    y = %g  -->  %g,     Æx = %g,  Æy = %g",px0,px1,py0,py1,px1-px0,py1-py0
//		peakUncorrect(xymap,px1,py1)
//		printf "     and invert to get (%g, %g),    Æ=(%g, %g)\r",px1,py1,px0-px1,py0-py1
//	
//		px0=942.054  ;  py0=1673.23
//		px1=px0  ;  py1=py0
//		peakcorrection2(xymap,px1,py1)
//		printf "x = %g  -->  %g,    y = %g  -->  %g,     Æx = %g,  Æy = %g",px0,px1,py0,py1,px1-px0,py1-py0
//		peakUncorrect(xymap,px1,py1)
//		printf "     and invert to get (%g, %g),    Æ=(%g, %g)\r",px1,py1,px0-px1,py0-py1
//	
//		px0=500  ;  py0=500
//		px1=px0  ;  py1=py0
//		peakcorrection2(xymap,px1,py1)
//		printf "x = %g  -->  %g,    y = %g  -->  %g,     Æx = %g,  Æy = %g",px0,px1,py0,py1,px1-px0,py1-py0
//		peakUncorrect(xymap,px1,py1)
//		printf "     and invert to get (%g, %g),    Æ=(%g, %g)\r",px1,py1,px0-px1,py0-py1
//	
//		px0=1042  ;  py0=1042
//		px1=px0  ;  py1=py0
//		peakcorrection2(xymap,px1,py1)
//		printf "x = %g  -->  %g,    y = %g  -->  %g,     Æx = %g,  Æy = %g",px0,px1,py0,py1,px1-px0,py1-py0
//		peakUncorrect(xymap,px1,py1)
//		printf "     and invert to get (%g, %g),    Æ=(%g, %g)\r",px1,py1,px0-px1,py0-py1
//	
//	
//		px0=93.09*2  ;  py0=143.96*2
//		px1=px0  ;  py1=py0
//		peakcorrection2(xymap,px1,py1)
//		printf "x = %g  -->  %g,    y = %g  -->  %g,     Æx = %g,  Æy = %g",px0,px1,py0,py1,px1-px0,py1-py0
//		peakUncorrect(xymap,px1,py1)
//		printf "     and invert to get (%g, %g),    Æ=(%g, %g)\r",px1,py1,px0-px1,py0-py1
//	End


// convert qvector to pixel position
Function/C q2pixel(geo,qvec)						// returns pixel position as a complex number cmplx(px,py)
	STRUCT microGeometry &geo
	Wave qvec										// qvec need not be normalized
	Variable px,py									// final pixel position, full chip unbinned 0 based pixels

	Wave qhat=root:Packages:geometry:q2pixel_qhat, ki=root:Packages:geometry:q2pixel_ki
	Wave kout=root:Packages:geometry:q2pixel_kout, vecB=root:Packages:geometry:q2pixel_vecB
	Wave mat=root:Packages:geometry:q2pixel_mat
	qhat = qvec
	normalize(qhat)

	ki = geo.ki[p]
	normalize(ki)
	Variable qLen = -2*MatrixDot(qhat,ki)			// length of qhat, note (q^ dot -ki) always positive
	kout = qhat*qLen + ki							// ko - ki = q
	normalize(kout)

	//	k k^ = a px^ + b py^ + dd {001}				// need to solve for a & b, (know all unit vectors and dd)
	//  -dd {001} = a px^ + b py^ - k k^			// rewrite this way to solve for a, b, k, (I know dd)
	// vecB = mat x {a,b,k}
	vecB = {0,0,-geo.dd}
	mat[0][0] = geo.rded00
	mat[1][0] = geo.rded01
	mat[2][0] = geo.rded02
	mat[0][1] = geo.rded10
	mat[1][1] = geo.rded11
	mat[2][1] = geo.rded12
	mat[][2] = -kout[p]
//	MatrixOp/O M_x = Inv(mat) x vecB
	MatrixLUD mat
	MatrixLUBkSub M_Lower, M_Upper, W_LUPermutation, vecB
	Wave M_x = M_x
	M_x[0] = CCD_X_REVERSE ? -M_x[0] : M_x[0]
	px = geo.xcent + M_x[0]*geo.NxCCD/geo.dpsx	// this formula is for 1 based pixels
	py = geo.ycent - M_x[1]*geo.NyCCD/geo.dpsy
	px -= 1											// convert from 1 based to 0 based pixels
	py -= 1
	Wave xymap=root:Packages:geometry:xymap
	peakUncorrect(xymap,px,py)	// go from true peak position to the measured peak position (put distortion back in), just invert peakcorrection2()

//	KillWaves/Z M_Lower,M_Upper,W_LUPermutation,M_x, kout,ki,qhat,vecB,mat
	return cmplx(px,py)
End



// convert xp,yp positions on screen into absolute coordinates x,y,z (Wenge coords)
Function pixel2q(geo,px,py,qhat)
	STRUCT microGeometry &geo
	Variable px,py						// pixel position, 0 based, first pixel is (0,0), NOT (1,1)
	Wave qhat							// q-vector in Wenge coordinates

	Wave ki=root:Packages:geometry:pixel2q_ki, kout=root:Packages:geometry:pixel2q_kout
	ki = geo.ki[p]									// incident beam direction
	pixel2XYZ(geo,px,py,kout)						// kout is in direction of pixel in xyz coords
	normalize(kout)

	Variable theta = acos(MatrixDot(kout,ki))/2	// Bragg angle (radians)
	if (WaveExists(qhat))
		qhat = kout-ki								// qhat bisects kout and -ki
		normalize(qhat)
	endif
	return theta
End
//
// convert px,py positions on detector into x,y,z (Wenge coords) with origin at Si position (µm)
Function pixel2XYZ(geo,px,py,xyz)
	STRUCT microGeometry &geo
	Variable px,py								// pixel is zero based on input, first pixel is (0,0), NOT (1,1)
	Wave xyz										// absolute position of pixel in Wenge coords (a 3-vector) (µm)

	Wave xymap=root:Packages:geometry:xymap
	peakcorrection2(xymap,px,py)				// convert pixel on CCD to a true distance in pixels (takes zero based pixels)

	Variable xxst,yyst
	xxst = (1+px-geo.xcent)*geo.dpsx/geo.NxCCD	// the "1+" is to make 1 based pixels
	xxst = CCD_X_REVERSE ? -xxst : xxst
	yyst = -(1+py-geo.ycent)*geo.dpsy/geo.NyCCD
	Variable xcc,ycc,ddc
	xcc = geo.rded00*xxst + geo.rded10*yyst
	ycc = geo.rded01*xxst + geo.rded11*yyst
	ddc = geo.rded02*xxst + geo.rded12*yyst	// dd is up
	xyz = {xcc,ycc,ddc+geo.dd}					// position in mm
	xyz *= 1000									// convert to µm
	return 0
End


// DEPRECATED		This routine should probably not be used.   And it is probably not used by anyone else either.
// convert pixel location xyz[3](µm) to pixel values
Function/C xyz2pixel(geo,xyz)						// returns pixel values as a complex number cmplx(px,py)
	STRUCT microGeometry &geo
	Wave xyz										// absolute position of pixel in Wenge coords (a 3-vector) (µm), origin at Si
	Variable px,py									// final pixel position, returned (zero based pixels)

	Wave vecB=root:Packages:geometry:q2pixel_vecB
	Wave mat=root:Packages:geometry:q2pixel_mat

	Variable xcc,ycc,ddc								//	xyz = {xcc,ycc,ddc+geo.dd}
	xcc = xyz[0]									// coordinates with origin at center of ccd, (xc,yc)
	ycc = xyz[1]
	ddc = xyz[2] - geo.dd*1000
	vecB = {xcc,ycc,ddc}
	vecB /= 1000									// convert from µm to mm

	mat[0][0] = geo.rded00
	mat[1][0] = geo.rded01
	mat[2][0] = geo.rded02
	mat[0][1] = geo.rded10
	mat[1][1] = geo.rded11
	mat[2][1] = geo.rded12
	mat[0][2] = geo.rded20
	mat[1][2] = geo.rded21
	mat[2][2] = geo.rded22

	MatrixOp/O vecInCCD = Inv(mat) x vecB
	vecInCCD[0] = CCD_X_REVERSE ? -vecInCCD[0] : vecInCCD[0]
	px = geo.xcent + vecInCCD[0]*geo.NxCCD/geo.dpsx	// here pixels are 1 based
	py = geo.ycent - vecInCCD[1]*geo.NyCCD/geo.dpsy
	KillWaves/Z vecInCCD
	px -= 1											// convert from 1 to 0 based pixels
	py -= 1

	Wave xymap=root:Packages:geometry:xymap
	peakUncorrect(xymap,px,py)	// go from true peak position to the measured peak position (put distortion back in), just invert peakcorrection2()
	return cmplx(px,py)							// return 0 based pixels
End



//correct the peak positon from measured to actual, Wenge Yang 5/19/2003
// This assumes a zero based coordinate on input, origin is (0,0), not (1,1)
Function/C peakcorrection2(xymap,px,py)		// returns cmplx(dx,dy)
	Wave xymap								// distortion map
	Variable &px,&py							// un-binned full frame pixel location on input, changed to distorted value at end

	Variable useDistortion = NumVarOrDefault("root:Packages:geometry:useDistortion",USE_DISTORTION_DEFAULT)
	if (!useDistortion || !WaveExists(xymap))	// do not distort
		return cmplx(0,0)
	endif
	Variable dx, dy								// the calculated corrections
	// here we start the calculation of distortion correction just like Wenge
	px += 1									// convert origin from (0,0) to (1,1)
	py += 1


	// This section added to speed up the distortion correction when correcting every pixel, not just a few fitted peaks
	if (exists("root:Packages:geometry:tempCachedDistortionMap")==1)	// this section is when using precomputed distortion calculations
		Wave distortionMap = root:Packages:geometry:tempCachedDistortionMap
		if (NumberByKey("use",note(distortionMap),"="))
			if (WaveDims(distortionMap)!=3 || DimSize(distortionMap,2)!=2)
				Abort "peakcorrection2(), The wave root:Packages:geometry:tempCachedDistortionMap, has an illegal size"
			endif
			Variable i,j, Ni = DimSize(distortionMap,0), Nj = DimSize(distortionMap,1)
			i = (px - DimOffset(distortionMap,0)) / DimDelta(distortionMap,0)
			j = (py - DimOffset(distortionMap,1)) / DimDelta(distortionMap,1)
			if (i<0 || j<0 || i>=Ni || j>=Nj)
				Abort "peakcorrection2(), 'root:Packages:geometry:tempCachedDistortionMap' has a scaling that does not fit input pixels"
			endif
			dx = distortionMap[i][j][0]
			dy = distortionMap[i][j][1]
			px += dx - 1							// add correction and convert origin from (1,1) to (0,0)
			py += dy - 1
			return cmplx(dx,dy)
		endif
	endif


	Variable nx=DimSize(xymap,0)-1, ny=DimSize(xymap,1)-1
	String noteStr = note(xymap)
	if (!stringmatch(StringByKey("CCDname",noteStr,"="),"White2084x2084"))
		px -= 1									// convert pixels back to 0 based
		py -= 1
		return cmplx(0,0)						// this routine is only for (2084 x 2084) CCD
	endif
	Variable cornerX0,cornerY0,cornerX1,cornerY1
	cornerX0 = NumberByKey("cornerX0",noteStr,"=")-1
	cornerY0 = NumberByKey("cornerY0",noteStr,"=")-1
	cornerX1 = NumberByKey("cornerX1",noteStr,"=")-1
	cornerY1 = NumberByKey("cornerY1",noteStr,"=")-1
//	printf "cornerX0=%d,  cornerY0=%d,  cornerX1=%d,  cornerY1=%d\r",cornerX0,cornerY0,cornerX1,cornerY1

	Make/O/D/N=4 cornerxy
	Make/N=4/O tempx,tempy
	tempx = {0,nx,0,nx}
	tempy = {0,0,ny,ny}
	cornerxy[0] = xymap[tempx[cornerX0]][tempy[cornerX0]][0]
	cornerxy[1] = xymap[tempx[cornerY0]][tempy[cornerY0]][1]
	cornerxy[2] = xymap[tempx[cornerX1]][tempy[cornerX1]][0]
	cornerxy[3] = xymap[tempx[cornerY1]][tempy[cornerY1]][1]
	KillWaves/Z tempx,tempy
//	printWave(cornerxy)

	Variable x0=xymap[0][0][0], y0=xymap[0][0][1]
	Variable xxstep, xystep, yxstep, yystep
	xxstep = (xymap[nx][0][0]-xymap[0][0][0]) / nx
	xystep = (xymap[nx][0][1]-xymap[0][0][1]) / nx
	yxstep = (xymap[0][ny][0]-xymap[0][0][0]) / ny
	yystep = (xymap[0][ny][1]-xymap[0][0][1]) / ny

	Variable xcent, ycent
//	xcent = (round(((px-x0)*yystep-(py-y0)*yxstep)/(xxstep*yystep-xystep*yxstep))<nx)>0
//	ycent = (round(((px-x0)*xystep-(py-y0)*xxstep)/(xystep*yxstep-yystep*xxstep))<ny)>0
	xcent = round(((px-x0)*yystep-(py-y0)*yxstep)/(xxstep*yystep-xystep*yxstep))
	xcent = limit(xcent,0,nx)
	ycent = round(((px-x0)*xystep-(py-y0)*xxstep)/(xystep*yxstep-yystep*xxstep))
	ycent = limit(ycent,0,ny)

//	printf "px=%g,  py=%g,  xymap[xcent=%d][ycent=%d][0]=%g,  xymap[%d][%d][2]=%g\r",px, py,xcent,ycent,xymap[xcent][ycent][0],xcent,ycent,xymap[xcent][ycent][1]

	dx=0  ;  dy=0								// initialize the correction
	if (px <= cornerxy[0] || px >= cornerxy[2] || py < cornerxy[1] || py >= cornerxy[3])
		dx = xymap[xcent][ycent][2]			// outside measured range, use default
		dy = xymap[xcent][ycent][3]
	else
		Variable xcent2, ycent2
		Variable r1,r2,r3,r4
		if (px >= xymap[xcent][ycent][0] && py >= xymap[xcent][ycent][1])
//			xcent2=(xcent+1)<nx>0
//			ycent2=(ycent+1)<ny>0
			xcent2 = limit(xcent+1,0,nx)
			ycent2 = limit(ycent+1,0,ny)
			r1=sqrt((px-xymap[xcent][ycent][0])^2+(py-xymap[xcent][ycent][1])^2)
			r2=sqrt((px-xymap[xcent][ycent2][0])^2+(py-xymap[xcent][ycent2][1])^2)
			r3=sqrt((px-xymap[xcent2][ycent][0])^2+(py-xymap[xcent2][ycent][1])^2)
			r4=sqrt((px-xymap[xcent2][ycent2][0])^2+(py-xymap[xcent2][ycent2][1])^2)
			dx=(xymap[xcent][ycent][2]/r1+xymap[xcent][ycent2][2]/r2+xymap[xcent2][ycent][2]/r3+xymap[xcent2][ycent2][2]/r4)
			dx=dx/(1/r1+1/r2+1/r3+1/r4)
			dy=(xymap[xcent][ycent][3]/r1+xymap[xcent][ycent2][3]/r2+xymap[xcent2][ycent][3]/r3+xymap[xcent2][ycent2][3]/r4)
			dy=dy/(1/r1+1/r2+1/r3+1/r4)
		endif
		if (px >= xymap[xcent][ycent][0] && py < xymap[xcent][ycent][1])
//			xcent2=(xcent+1)<nx>0
//			ycent2=(ycent-1)<ny>0
			xcent2 = limit(xcent+1,0,nx)
			ycent2 = limit(ycent-1,0,ny)
			r1=sqrt((px-xymap[xcent][ycent][0])^2+(py-xymap[xcent][ycent][1])^2)
			r2=sqrt((px-xymap[xcent][ycent2][0])^2+(py-xymap[xcent][ycent2][1])^2)
			r3=sqrt((px-xymap[xcent2][ycent][0])^2+(py-xymap[xcent2][ycent][1])^2)
			r4=sqrt((px-xymap[xcent2][ycent2][0])^2+(py-xymap[xcent2][ycent2][1])^2)
			dx=(xymap[xcent][ycent][2]/r1+xymap[xcent][ycent2][2]/r2+xymap[xcent2][ycent][2]/r3+xymap[xcent2][ycent2][2]/r4)
			dx=dx/(1/r1+1/r2+1/r3+1/r4)
			dy=(xymap[xcent][ycent][3]/r1+xymap[xcent][ycent2][3]/r2+xymap[xcent2][ycent][3]/r3+xymap[xcent2][ycent2][3]/r4)
			dy=dy/(1/r1+1/r2+1/r3+1/r4)
		endif
		if (px < xymap[xcent][ycent][0] && py >= xymap[xcent][ycent][1])
//			xcent2=(xcent-1)<nx>0
//			ycent2=(ycent+1)<ny>0
			xcent2 = limit(xcent-1,0,nx)
			ycent2 = limit(ycent+1,0,ny)
			r1=sqrt((px-xymap[xcent][ycent][0])^2+(py-xymap[xcent][ycent][1])^2)
			r2=sqrt((px-xymap[xcent][ycent2][0])^2+(py-xymap[xcent][ycent2][1])^2)
			r3=sqrt((px-xymap[xcent2][ycent][0])^2+(py-xymap[xcent2][ycent][1])^2)
			r4=sqrt((px-xymap[xcent2][ycent2][0])^2+(py-xymap[xcent2][ycent2][1])^2)
			dx=(xymap[xcent][ycent][2]/r1+xymap[xcent][ycent2][2]/r2+xymap[xcent2][ycent][2]/r3+xymap[xcent2][ycent2][2]/r4)
			dx=dx/(1/r1+1/r2+1/r3+1/r4)
			dy=(xymap[xcent][ycent][3]/r1+xymap[xcent][ycent2][3]/r2+xymap[xcent2][ycent][3]/r3+xymap[xcent2][ycent2][3]/r4)
			dy=dy/(1/r1+1/r2+1/r3+1/r4)
		endif
		if (px < xymap[xcent][ycent][0] && py < xymap[xcent][ycent][1])
//			xcent2=(xcent-1)<nx>0
//			ycent2=(ycent-1)<ny>0
			xcent2 = limit(xcent-1,0,nx)
			ycent2 = limit(ycent-1,00,ny)
			r1=sqrt((px-xymap[xcent][ycent][0])^2+(py-xymap[xcent][ycent][1])^2)
			r2=sqrt((px-xymap[xcent][ycent2][0])^2+(py-xymap[xcent][ycent2][1])^2)
			r3=sqrt((px-xymap[xcent2][ycent][0])^2+(py-xymap[xcent2][ycent][1])^2)
			r4=sqrt((px-xymap[xcent2][ycent2][0])^2+(py-xymap[xcent2][ycent2][1])^2)
			dx=(xymap[xcent][ycent][2]/r1+xymap[xcent][ycent2][2]/r2+xymap[xcent2][ycent][2]/r3+xymap[xcent2][ycent2][2]/r4)
			dx=dx/(1/r1+1/r2+1/r3+1/r4)
			dy=(xymap[xcent][ycent][3]/r1+xymap[xcent][ycent2][3]/r2+xymap[xcent2][ycent][3]/r3+xymap[xcent2][ycent2][3]/r4)
			dy=dy/(1/r1+1/r2+1/r3+1/r4)
		endif
	endif
	px += dx - 1								// add correction and convert origin from (1,1) to (0,0)
	py += dy - 1
	KillWaves/Z cornerxy
	return cmplx(dx,dy)
End
//
//	//correct the peak positon from measured to actual, Wenge Yang 5/19/2003
//	// This assumes a zero based coordinate on input, origin is (0,0), not (1,1)
//	Function/C peakcorrection2(xymap,px,py)		// returns cmplx(dx,dy)
//		Wave xymap								// distortion map
//		Variable &px,&py							// un-binned full frame pixel location on input, changed to distorted value at end
//	
//		Variable useDistortion = NumVarOrDefault("root:Packages:geometry:useDistortion",USE_DISTORTION_DEFAULT)
//		if (!useDistortion || !WaveExists(xymap))	// do not distort
//			return cmplx(0,0)
//		endif
//		px += 1									// convert origin from (0,0) to (1,1)
//		py += 1
//	
//		Variable nx=DimSize(xymap,0)-1, ny=DimSize(xymap,1)-1
//		String noteStr = note(xymap)
//		if (!stringmatch(StringByKey("CCDname",noteStr,"="),"White2084x2084"))
//			return cmplx(0,0)						// this routine is only for (2084 x 2084) CCD
//		endif
//		Variable cornerX0,cornerY0,cornerX1,cornerY1
//		cornerX0 = NumberByKey("cornerX0",noteStr,"=")-1
//		cornerY0 = NumberByKey("cornerY0",noteStr,"=")-1
//		cornerX1 = NumberByKey("cornerX1",noteStr,"=")-1
//		cornerY1 = NumberByKey("cornerY1",noteStr,"=")-1
//	//	printf "cornerX0=%d,  cornerY0=%d,  cornerX1=%d,  cornerY1=%d\r",cornerX0,cornerY0,cornerX1,cornerY1
//	
//		Make/O/D/N=4 cornerxy
//		Make/N=4/O tempx,tempy
//		tempx = {0,nx,0,nx}
//		tempy = {0,0,ny,ny}
//		cornerxy[0] = xymap[tempx[cornerX0]][tempy[cornerX0]][0]
//		cornerxy[1] = xymap[tempx[cornerY0]][tempy[cornerY0]][1]
//		cornerxy[2] = xymap[tempx[cornerX1]][tempy[cornerX1]][0]
//		cornerxy[3] = xymap[tempx[cornerY1]][tempy[cornerY1]][1]
//		KillWaves/Z tempx,tempy
//	//	printWave(cornerxy)
//	
//		Variable x0=xymap[0][0][0], y0=xymap[0][0][1]
//		Variable xxstep, xystep, yxstep, yystep
//		xxstep = (xymap[nx][0][0]-xymap[0][0][0]) / nx
//		xystep = (xymap[nx][0][1]-xymap[0][0][1]) / nx
//		yxstep = (xymap[0][ny][0]-xymap[0][0][0]) / ny
//		yystep = (xymap[0][ny][1]-xymap[0][0][1]) / ny
//	
//		Variable xcent, ycent
//	//	xcent = (round(((px-x0)*yystep-(py-y0)*yxstep)/(xxstep*yystep-xystep*yxstep))<nx)>0
//	//	ycent = (round(((px-x0)*xystep-(py-y0)*xxstep)/(xystep*yxstep-yystep*xxstep))<ny)>0
//		xcent = round(((px-x0)*yystep-(py-y0)*yxstep)/(xxstep*yystep-xystep*yxstep))
//		xcent = limit(xcent,0,nx)
//		ycent = round(((px-x0)*xystep-(py-y0)*xxstep)/(xystep*yxstep-yystep*xxstep))
//		ycent = limit(ycent,0,ny)
//	
//	//	printf "px=%g,  py=%g,  xymap[xcent=%d][ycent=%d][0]=%g,  xymap[%d][%d][2]=%g\r",px, py,xcent,ycent,xymap[xcent][ycent][0],xcent,ycent,xymap[xcent][ycent][1]
//	
//		Variable dx, dy								// the corrections
//		if (px <= cornerxy[0] || px >= cornerxy[2] || py < cornerxy[1] || py >= cornerxy[3])
//			dx = xymap[xcent][ycent][2]			// outside measured range, use default
//			dy = xymap[xcent][ycent][3]
//		else
//			Variable xcent2, ycent2
//			Variable r1,r2,r3,r4
//			if (px >= xymap[xcent][ycent][0] && py >= xymap[xcent][ycent][1])
//	//			xcent2=(xcent+1)<nx>0
//	//			ycent2=(ycent+1)<ny>0
//				xcent2 = limit(xcent+1,0,nx)
//				ycent2 = limit(ycent+1,0,ny)
//				r1=sqrt((px-xymap[xcent][ycent][0])^2+(py-xymap[xcent][ycent][1])^2)
//				r2=sqrt((px-xymap[xcent][ycent2][0])^2+(py-xymap[xcent][ycent2][1])^2)
//				r3=sqrt((px-xymap[xcent2][ycent][0])^2+(py-xymap[xcent2][ycent][1])^2)
//				r4=sqrt((px-xymap[xcent2][ycent2][0])^2+(py-xymap[xcent2][ycent2][1])^2)
//				dx=(xymap[xcent][ycent][2]/r1+xymap[xcent][ycent2][2]/r2+xymap[xcent2][ycent][2]/r3+xymap[xcent2][ycent2][2]/r4)
//				dx=dx/(1/r1+1/r2+1/r3+1/r4)
//				dy=(xymap[xcent][ycent][3]/r1+xymap[xcent][ycent2][3]/r2+xymap[xcent2][ycent][3]/r3+xymap[xcent2][ycent2][3]/r4)
//				dy=dy/(1/r1+1/r2+1/r3+1/r4)
//			endif
//			if (px >= xymap[xcent][ycent][0] && py < xymap[xcent][ycent][1])
//	//			xcent2=(xcent+1)<nx>0
//	//			ycent2=(ycent-1)<ny>0
//				xcent2 = limit(xcent+1,0,nx)
//				ycent2 = limit(ycent-1,0,ny)
//				r1=sqrt((px-xymap[xcent][ycent][0])^2+(py-xymap[xcent][ycent][1])^2)
//				r2=sqrt((px-xymap[xcent][ycent2][0])^2+(py-xymap[xcent][ycent2][1])^2)
//				r3=sqrt((px-xymap[xcent2][ycent][0])^2+(py-xymap[xcent2][ycent][1])^2)
//				r4=sqrt((px-xymap[xcent2][ycent2][0])^2+(py-xymap[xcent2][ycent2][1])^2)
//				dx=(xymap[xcent][ycent][2]/r1+xymap[xcent][ycent2][2]/r2+xymap[xcent2][ycent][2]/r3+xymap[xcent2][ycent2][2]/r4)
//				dx=dx/(1/r1+1/r2+1/r3+1/r4)
//				dy=(xymap[xcent][ycent][3]/r1+xymap[xcent][ycent2][3]/r2+xymap[xcent2][ycent][3]/r3+xymap[xcent2][ycent2][3]/r4)
//				dy=dy/(1/r1+1/r2+1/r3+1/r4)
//			endif
//			if (px < xymap[xcent][ycent][0] && py >= xymap[xcent][ycent][1])
//	//			xcent2=(xcent-1)<nx>0
//	//			ycent2=(ycent+1)<ny>0
//				xcent2 = limit(xcent-1,0,nx)
//				ycent2 = limit(ycent+1,0,ny)
//				r1=sqrt((px-xymap[xcent][ycent][0])^2+(py-xymap[xcent][ycent][1])^2)
//				r2=sqrt((px-xymap[xcent][ycent2][0])^2+(py-xymap[xcent][ycent2][1])^2)
//				r3=sqrt((px-xymap[xcent2][ycent][0])^2+(py-xymap[xcent2][ycent][1])^2)
//				r4=sqrt((px-xymap[xcent2][ycent2][0])^2+(py-xymap[xcent2][ycent2][1])^2)
//				dx=(xymap[xcent][ycent][2]/r1+xymap[xcent][ycent2][2]/r2+xymap[xcent2][ycent][2]/r3+xymap[xcent2][ycent2][2]/r4)
//				dx=dx/(1/r1+1/r2+1/r3+1/r4)
//				dy=(xymap[xcent][ycent][3]/r1+xymap[xcent][ycent2][3]/r2+xymap[xcent2][ycent][3]/r3+xymap[xcent2][ycent2][3]/r4)
//				dy=dy/(1/r1+1/r2+1/r3+1/r4)
//			endif
//			if (px < xymap[xcent][ycent][0] && py < xymap[xcent][ycent][1])
//	//			xcent2=(xcent-1)<nx>0
//	//			ycent2=(ycent-1)<ny>0
//				xcent2 = limit(xcent-1,0,nx)
//				ycent2 = limit(ycent-1,00,ny)
//				r1=sqrt((px-xymap[xcent][ycent][0])^2+(py-xymap[xcent][ycent][1])^2)
//				r2=sqrt((px-xymap[xcent][ycent2][0])^2+(py-xymap[xcent][ycent2][1])^2)
//				r3=sqrt((px-xymap[xcent2][ycent][0])^2+(py-xymap[xcent2][ycent][1])^2)
//				r4=sqrt((px-xymap[xcent2][ycent2][0])^2+(py-xymap[xcent2][ycent2][1])^2)
//				dx=(xymap[xcent][ycent][2]/r1+xymap[xcent][ycent2][2]/r2+xymap[xcent2][ycent][2]/r3+xymap[xcent2][ycent2][2]/r4)
//				dx=dx/(1/r1+1/r2+1/r3+1/r4)
//				dy=(xymap[xcent][ycent][3]/r1+xymap[xcent][ycent2][3]/r2+xymap[xcent2][ycent][3]/r3+xymap[xcent2][ycent2][3]/r4)
//				dy=dy/(1/r1+1/r2+1/r3+1/r4)
//			endif
//		endif
//		px += dx - 1								// add correction and convert origin from (1,1) to (0,0)
//		py += dy - 1
//		KillWaves/Z cornerxy
//		return cmplx(dx,dy)
//	End
//
Function/C peakUncorrect(xymap,px,py)		// go from true peak to the measured peak position (put distortion back in), just invert peakcorrection2()
	Wave xymap								// distortion map
	Variable &px,&py							// un-binned full frame 0 based pixel with distortion correction, changed to un-distorted value at end

	Variable useDistortion = NumVarOrDefault("root:Packages:geometry:useDistortion",USE_DISTORTION_DEFAULT)
	if (!useDistortion || !WaveExists(xymap))	// do not distort
		return cmplx(px,py)
	endif

	Variable maxIter=6, tol=1e-4				// at most 6 iterations, and a tolerance of 1e-4

	Variable pxd=px, pyd=py					// save the starting point (the required distorted position)
	Variable/C dz
	Variable i=0,err = Inf
	for (i=0; i<maxIter && err>tol; i+=1)		// end when it changes by less than tol or after maxIter iterations
		dz = peakcorrection2(xymap,px,py)		// returns changed px,py
		err = abs(px-pxd)+abs(py-pyd)
		px = pxd - real(dz)
		py = pyd - imag(dz)
	endfor
	return cmplx(px,py)
End



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



//Constant NxCCD=2084, NyCCD=2084		// # of un-binned pixels in CCD, also called xdim, ydim
//Constant dpsx=50010, dpsy=50168			// size of CCD (µm)
//Constant xcent=1091.53, ycent=939.941	// center of CCD in un-binned pixels
//Constant dd=30017.3						// height of CCD (µm)
//Constant xbet=0.31037, xgam=0.75037	// incident beam tilt (degrees)
//Constant xalfd=0.10525, xbetd=0.171		// CCD tilt (degrees)

// output from test(), use these to check the calculation
// test()
//	using parameters:
//		NxCCD=2084, NyCCD=2084		// # of un-binned pixels in CCD, also called xdim, ydim
//		dpsx=50010, dpsy=50168			// size of CCD (µm)
//		xcent=1091.53, ycent=939.941	// center of CCD in un-binned pixels
//		dd=30017.3						// height of CCD (µm)
//		xbet=0.31037, xgam=0.75037		// incident beam tilt (degrees)
//		xalfd=0.10525, xbetd=0.171		// CCD tilt (degrees)
//
//	theta(1024, 1024) = 42.8989¡
//	theta(1098.4, 917.88) = 45.3536¡
//	theta(100, 1500) = 34.9344¡
//	theta(200, 1600) = 32.9765¡
//	theta(300, 350) = 55.5618¡








// calculate rotation matrix from detector tilts angles
// note: in IDL matrices are defined as transpose of the usual matrices
// 
// usese geo to set ki[3], rded[3][3], and rde[3][3]
Function GeometryUpdateCalc(geo)			// formerly callled detectortilt()
	STRUCT microGeometry &geo

	Init_microGeo()
	// first convert all angles from degrees --> radians
	Variable xbr = -geo.xbet*PI/180, xcr=-geo.xgam*PI/180
	Variable xard=geo.xalfd*PI/180, xbrd=geo.xbetd*PI/180

	// rotation matrix for incident beam, use it to set ki[]
	//
	// rot(xbr) =	1		0			0
	//				0	  cos(xbr)	sin(xbr)
	//				0	-sin(xbr)	cos(xbr)
	// rot(xbr) rotates a vector about x negative. So +xbetd move the beam upward as it should.
	//
	// rot(xcr) =	  cos(xcr)	sin(xcr)	0
	//				-sin(xcr)	cos(xcr)	0
	//				  	0			0		1
	// rot(xcr) rotates a vector about z (which points upward) in a negative direction.  So +xgam is backward.
	//
	// rde = rot(xbr) x root(xcr)
	//
	Wave rde = root:Packages:geometry:GeoUpdate_rde
	rde[0][0] = cos(xcr)							// rotation matrix for incident beam
	rde[1][0] = -sin(xcr)*cos(xbr)
	rde[2][0] = sin(xbr)*sin(xcr)
	rde[0][1] = sin(xcr)
	rde[1][1] = cos(xbr)*cos(xcr)
	rde[2][1] = -cos(xcr)*sin(xbr)
	rde[0][2] = 0
	rde[1][2] = sin(xbr)
	rde[2][2] = cos(xbr)
	Wave ki = root:Packages:geometry:GeoUpdate_ki
//	ki = -rde[1][p]									// ki = (0,-1,0) x rde
	ki = {0,-1,0}
	MatrixMultiply ki/T, rde
	Wave M_product=M_product
	ki = M_product
	KillWaves M_product
	normalize(ki)
	geo.ki[0] = ki[0]								// in Wenge coordinates (X=out door, Y=upstream, Z=up)
	geo.ki[1] = ki[1]								// note, cannot automatically iterate a vector assignment when it is in a structure
	geo.ki[2] = ki[2]

	geo.rde00 =   rde[0][0]							// rotation matrix that goes from Wenge system to IDEAL Beam-Line
	geo.rde10 =   rde[2][0]							// this is like rde[][], but includes the swap of Y=-z and Z=y
	geo.rde20 = -rde[1][0]							// NOTE, this is a rotation matrix, so its Transpose is its inverse
	geo.rde01 =   rde[0][1]
	geo.rde11 =   rde[2][1]
	geo.rde21 = -rde[1][1]
	geo.rde02 =   rde[0][2]
	geo.rde12 =   rde[2][2]
	geo.rde22 = -rde[1][2]
	//				// an example of how to use geo.rdexx
	//			qBL[0] = qW[0]*geo.rde00 + qW[1]*geo.rde01 + qW[2]*geo.rde02	// transform Wenge to Ideal Beam-line
	//			qBL[1] = qW[0]*geo.rde10 + qW[1]*geo.rde11 + qW[2]*geo.rde12	//   with incident beam along {0,0,1}
	//			qBL[2] = qW[0]*geo.rde20 + qW[1]*geo.rde21 + qW[2]*geo.rde22

	// rotation matrix for detector
	//
	// rot(a) =  cos(a)		0	  sin(a)			// rotation about y axis
	//			  	0			1		1
	//			-sin(a)		0	  cos(a)
	//
	// rot(b) =	1		0			0				// rotation about x axis
	//				0	  cos(xbr)	sin(xbr)
	//				0	-sin(xbr)	cos(xbr)
	//
	// rded = rot(xard) x root(xbrd)
	//
	geo.rded00 = cos(xard)
	geo.rded10 = 0
	geo.rded20 = -sin(xard)
	geo.rded01 = -sin(xard)*sin(xbrd)
	geo.rded11 = cos(xbrd)
	geo.rded21 = -sin(xbrd)*cos(xard)
	geo.rded02 = sin(xard)*cos(xbrd)
	geo.rded12 = sin(xbrd)
	geo.rded22 = cos(xard)*cos(xbrd)

	Variable len = sqrt(geo.wire.axis[0]*geo.wire.axis[0] + geo.wire.axis[1]*geo.wire.axis[1] + geo.wire.axis[2]*geo.wire.axis[2])
	if (numtype(len)==0)
		geo.wire.axis[0] /= len  ;  geo.wire.axis[1] /= len  ;  geo.wire.axis[2] /= len
	else
		geo.wire.axis[0] = 1  ;  geo.wire.axis[1] = 0  ;  geo.wire.axis[2] = 0
	endif

	// compute wire coordinates that would put wire center at the Si position (µm)
	Variable ir2=1/sqrt(2), dZ
	Variable H0=geo.wire.H0, Hyc=geo.wire.Hyc, F=geo.wire.F

	geo.wire.xyz0[0] = geo.wire.X						// wire coordinates when wire is on the incident beam, computed (µm)
	geo.wire.xyz0[1] = (H0-F)*ir2
	geo.wire.xyz0[2] = (H0+F)*ir2

	geo.wire.xyzSi[2] = (Hyc+F)*ir2					// Z of the Si position in wire units & beam line coordinates (µm)
	dZ = (Hyc-H0)*ir2 								// dist (along z) from where wire intersects beam to Si
	geo.wire.xyzSi[1] = (H0-F)*ir2 - dZ*ki[2]/ki[1]	// Y of the Si in wire units & beam line coordinates (µm)
	geo.wire.xyzSi[0] = geo.wire.X - dZ*ki[0]/ki[1]	// X of the Si in wire units & beam line coordinates (µm)
End

//Function MakeGeometryParametersPanel(strStruct)
//	String strStruct									// optional passed value of geo structure, this is used if passed
//	if (strlen(WinList("GeometrySet",";","WIN:64")))	// window alreay exits, bring it to front
//		DoWindow/F GeometrySet
//		return 0
//	endif
//
//	Init_microGeo()
//	NewDataFolder/O root:Packages:geometry:PanelValues// ensure that the needed data folders exist
//	// create globals for the panel
//	Variable/G root:Packages:geometry:PanelValues:xalfd
//	Variable/G root:Packages:geometry:PanelValues:xbetd
//	Variable/G root:Packages:geometry:PanelValues:dd
//	Variable/G root:Packages:geometry:PanelValues:ddOffset
//	Variable/G root:Packages:geometry:PanelValues:NxCCD
//	Variable/G root:Packages:geometry:PanelValues:NyCCD
//	Variable/G root:Packages:geometry:PanelValues:dpsx
//	Variable/G root:Packages:geometry:PanelValues:dpsy
//	Variable/G root:Packages:geometry:PanelValues:xcent
//	Variable/G root:Packages:geometry:PanelValues:ycent
//	Variable/G root:Packages:geometry:PanelValues:xbet
//	Variable/G root:Packages:geometry:PanelValues:xgam
//	Variable/G root:Packages:geometry:PanelValues:dia
//	Variable/G root:Packages:geometry:PanelValues:H0
//	Variable/G root:Packages:geometry:PanelValues:Hyc
//	Variable/G root:Packages:geometry:PanelValues:F
//	NVAR xalfd = root:Packages:geometry:PanelValues:xalfd
//	NVAR xbetd = root:Packages:geometry:PanelValues:xbetd
//	NVAR dd = root:Packages:geometry:PanelValues:dd
//	NVAR ddOffset = root:Packages:geometry:PanelValues:ddOffset
//	NVAR NxCCD = root:Packages:geometry:PanelValues:NxCCD
//	NVAR NyCCD = root:Packages:geometry:PanelValues:NyCCD
//	NVAR dpsx = root:Packages:geometry:PanelValues:dpsx
//	NVAR dpsy = root:Packages:geometry:PanelValues:dpsy
//	NVAR xcent = root:Packages:geometry:PanelValues:xcent
//	NVAR ycent = root:Packages:geometry:PanelValues:ycent
//	NVAR xbet = root:Packages:geometry:PanelValues:xbet
//	NVAR xgam = root:Packages:geometry:PanelValues:xgam
//	NVAR dia = root:Packages:geometry:PanelValues:dia
//	NVAR H0 = root:Packages:geometry:PanelValues:H0
//	NVAR Hyc = root:Packages:geometry:PanelValues:Hyc
//	NVAR F = root:Packages:geometry:PanelValues:F
//
//	STRUCT microGeometry geo
//	geo.NxCCD = 0
//	if (strlen(strStruct)>0)
//		StructGet/S/B=2 geo, strStruct					// found structure information, load into geo
//		Variable/G root:Packages:geometry:PanelValues:dirty=1	// for passed structure
//	else
//		FillGeometryStructDefault(geo)					//fill the geometry structure with current values
//		Variable/G root:Packages:geometry:PanelValues:dirty=0
//	endif
//	xalfd = geo.xalfd										// set the panel globals from the struct
//	xbetd = geo.xbetd
//	dd = geo.dd
//	ddOffset = geo.ddOffset
//	NxCCD = geo.NxCCD
//	NyCCD = geo.NyCCD
//	dpsx = geo.dpsx
//	dpsy = geo.dpsy
//	xcent = geo.xcent
//	ycent = geo.ycent
//	xbet = geo.xbet
//	xgam = geo.xgam
//	dia = geo.wire.dia
//	H0 = geo.wire.H0
//	Hyc = geo.wire.Hyc
//	F = geo.wire.F
//
//	NewPanel /K=1 /W=(675,60,895,653)
//	DoWindow/C GeometrySet
//	SetDrawLayer UserBack
//	SetDrawEnv fname= "Lucida Grande",fsize= 16
//	DrawText 30,18,"Detector Position / tilt"
//	SetDrawEnv fname= "Lucida Grande",fsize= 16
//	DrawText 30,162,"CCD size / resolution"
//	SetDrawEnv fname= "Lucida Grande",fsize= 16
//	DrawText 30,263,"Incident beam direction"
//	SetDrawEnv fname= "Lucida Grande",fsize= 16
//	DrawText 30,323,"Scanning wire"
//	SetVariable xalfd,pos={78,20},size={92,15},proc=GeoPanelVarChangedProc,title="xalfd¡"
//	SetVariable xalfd,help={"alpha tilt of detector (degree)"}
//	SetVariable xalfd,limits={-1,1,0},value= root:Packages:geometry:PanelValues:xalfd,bodyWidth= 60
//	SetVariable xbetd,pos={71,40},size={99,15},proc=GeoPanelVarChangedProc,title="xbetad¡"
//	SetVariable xbetd,help={"beta tilt of detector (degree)"}
//	SetVariable xbetd,limits={-1,1,0},value= root:Packages:geometry:PanelValues:xbetd,bodyWidth= 60
//	SetVariable dd,pos={69,60},size={101,15},proc=GeoPanelVarChangedProc,title="dd (mm)"
//	SetVariable dd,help={"height of detector (mm)"}
//	SetVariable dd,limits={0,2000,0},value= root:Packages:geometry:PanelValues:dd,bodyWidth= 60
//	SetVariable ddOffset,pos={38,80},size={132,15},proc=GeoPanelVarChangedProc,title="dd Offset (mm)"
//	SetVariable ddOffset,help={"offset between motor and dd,  dd = CCDy+ddOffset"}
//	SetVariable ddOffset,limits={-100,100,0},value= root:Packages:geometry:PanelValues:ddOffset,bodyWidth= 60
//	SetVariable xcent,pos={48,100},size={122,15},proc=GeoPanelVarChangedProc,title="xcent (pixel)",format="%.3f"
//	SetVariable xcent,help={"central pixel in CCD X, 1 based & unbinned"}
//	SetVariable xcent,limits={0,5000,0},value= root:Packages:geometry:PanelValues:xcent,bodyWidth= 60
//	SetVariable ycent,pos={48,120},size={122,15},proc=GeoPanelVarChangedProc,title="ycent (pixel)",format="%.3f"
//	SetVariable ycent,help={"central pixel in CCD y, 1 based & unbinned"}
//	SetVariable ycent,limits={0,5000,0},value= root:Packages:geometry:PanelValues:ycent,bodyWidth= 60
//	SetVariable NxCCD,pos={76,165},size={94,15},proc=GeoPanelVarChangedProc,title="NxCCD"
//	SetVariable NxCCD,help={"x size of CCD (pixels)"},format="%d"
//	SetVariable NxCCD,limits={1,5000,0},value= root:Packages:geometry:PanelValues:NxCCD,bodyWidth= 60
//	SetVariable NyCCD,pos={76,185},size={94,15},proc=GeoPanelVarChangedProc,title="NyCCD"
//	SetVariable NyCCD,help={"y size of CCD (pixels)"},format="%d"
//	SetVariable NyCCD,limits={1,5000,0},value= root:Packages:geometry:PanelValues:NyCCD,bodyWidth= 60
//	SetVariable dpsx,pos={58,204},size={112,15},proc=GeoPanelVarChangedProc,title="dpsx (mm)"
//	SetVariable dpsx,help={"x size of CCD (mm)"}
//	SetVariable dpsx,limits={0,500,0},value= root:Packages:geometry:PanelValues:dpsx,bodyWidth= 60
//	SetVariable dpsy,pos={58,222},size={112,15},proc=GeoPanelVarChangedProc,title="dpsy (mm)"
//	SetVariable dpsy,help={"y size of CCD (mm)"}
//	SetVariable dpsy,limits={0,500,0},value= root:Packages:geometry:PanelValues:dpsy,bodyWidth= 60
//	SetVariable xbet,pos={76,265},size={94,15},proc=GeoPanelVarChangedProc,title="xbeta¡"
//	SetVariable xbet,help={"beta tilt of incident beam (degree)"}
//	SetVariable xbet,limits={-1,1,0},value= root:Packages:geometry:PanelValues:xbet,bodyWidth= 60
//	SetVariable xgam,pos={77,285},size={93,15},proc=GeoPanelVarChangedProc,title="xgam¡"
//	SetVariable xgam,help={"gamma tilt of incident beam (degree)"}
//	SetVariable xgam,limits={-12,12,0},value= root:Packages:geometry:PanelValues:xgam,bodyWidth= 60
//	SetVariable dia,pos={44,325},size={126,15},proc=GeoPanelVarChangedProc,title="wire dia (µm)"
//	SetVariable dia,help={"diameter of wire (micron)"}
//	SetVariable dia,limits={1,500,0},value= root:Packages:geometry:PanelValues:dia,bodyWidth= 60
//	SetVariable H0,pos={45,345},size={125,15},proc=GeoPanelVarChangedProc,title="wire H0 (µm)"
//	SetVariable H0,help={"H of wire where it cuts incident beam"}
//	SetVariable H0,limits={-50000,50000,0},value= root:Packages:geometry:PanelValues:H0,bodyWidth= 60
//	SetVariable Hyc,pos={40,365},size={130,15},proc=GeoPanelVarChangedProc,title="wire Hyc (µm)"
//	SetVariable Hyc,help={"H of wire where it cuts vetical ray to (xc,yc)"}
//	SetVariable Hyc,limits={-50000,50000,0},value= root:Packages:geometry:PanelValues:Hyc,bodyWidth= 60
//	SetVariable wireF,pos={52,385},size={118,15},proc=GeoPanelVarChangedProc,title="wire F (µm)"
//	SetVariable wireF,help={"F of wire at H0 & Hyc"}
//	SetVariable wireF,limits={-50000,50000,0},value= root:Packages:geometry:PanelValues:F,bodyWidth= 60
//	Button buttonSaveDefault,pos={35,415},size={150,20},proc=GeometryPanelButtonProc,title="Save as Default"
//	Button buttonSaveDefault,help={"saves these values as default in the current data folder"}
//	Button buttonSaveLocal,pos={35,440},size={150,20},proc=GeometryPanelButtonProc,title="Save in Current Fldr"
//	Button buttonSaveLocal,help={"saves these values as default in the current data folder"}
//	Button buttonEPICS,pos={35,465},size={150,20},proc=GeometryPanelButtonProc,title="Load from EPICS"
//	Button buttonEPICS,help={"Fill the values in this panel from EPICS"}
//	Button buttonFillFromFile,pos={35,490},size={150,20},proc=GeometryPanelButtonProc,title="Load from a file"
//	Button buttonFillFromFile,help={"Fill the values in this panel from a file"}
//	Button buttonWriteToFile,pos={35,515},size={150,20},proc=GeometryPanelButtonProc,title="Write to a file"
//	Button buttonWriteToFile,help={"Write values in this panel to a file"}
//	Button buttonPrintGeo,pos={35,540},size={150,20},proc=GeometryPanelButtonProc,title="Print Geo to History"
//	Button buttonPrintGeo,help={"print geometry values to the history"}
//	Button buttonCancel,pos={35,565},size={150,20},proc=GeometryPanelButtonProc,title="Cancel"
//	Button buttonCancel,help={"Close this panel, do not save any changes"}
//	Button buttonSetdd,pos={171,57},size={45,18},proc=GeometryPanelButtonProc,title="image"
//	Button buttonSetdd,help={"set dd using CCDy from top image"}
//	SetWindow kwTopWin,hook(close)=SetGeoPanelHook
//End
Function MakeGeometryParametersPanel(strStruct)
	String strStruct									// optional passed value of geo structure, this is used if passed
	if (strlen(WinList("GeometrySet",";","WIN:64")))	// window alreay exits, bring it to front
		DoWindow/F GeometrySet
		return 0
	endif
	NewPanel /K=1 /W=(675,60,895,653)
	DoWindow/C GeometrySet
	FillGeometryParametersPanel(strStruct,"GeometrySet",0,0)
	SetWindow kwTopWin,hook(close)=SetGeoPanelHook
End
//
Function/T FillGeometryParametersPanel(strStruct,hostWin,left,top)
	String strStruct									// optional passed value of geo structure, this is used if passed
	String hostWin										// name of home window
	Variable left, top									// offsets from the left and top

	Init_microGeo()
	NewDataFolder/O root:Packages:geometry:PanelValues// ensure that the needed data folders exist
	// create globals for the panel
	Variable/G root:Packages:geometry:PanelValues:xalfd
	Variable/G root:Packages:geometry:PanelValues:xbetd
	Variable/G root:Packages:geometry:PanelValues:dd
	Variable/G root:Packages:geometry:PanelValues:ddOffset
	Variable/G root:Packages:geometry:PanelValues:NxCCD
	Variable/G root:Packages:geometry:PanelValues:NyCCD
	Variable/G root:Packages:geometry:PanelValues:dpsx
	Variable/G root:Packages:geometry:PanelValues:dpsy
	Variable/G root:Packages:geometry:PanelValues:xcent
	Variable/G root:Packages:geometry:PanelValues:ycent
	Variable/G root:Packages:geometry:PanelValues:xbet
	Variable/G root:Packages:geometry:PanelValues:xgam
	Variable/G root:Packages:geometry:PanelValues:dia
	Variable/G root:Packages:geometry:PanelValues:H0
	Variable/G root:Packages:geometry:PanelValues:Hyc
	Variable/G root:Packages:geometry:PanelValues:F
	Variable/G root:Packages:geometry:PanelValues:axisX
	Variable/G root:Packages:geometry:PanelValues:axisY
	Variable/G root:Packages:geometry:PanelValues:axisZ

	NVAR xalfd = root:Packages:geometry:PanelValues:xalfd
	NVAR xbetd = root:Packages:geometry:PanelValues:xbetd
	NVAR dd = root:Packages:geometry:PanelValues:dd
	NVAR ddOffset = root:Packages:geometry:PanelValues:ddOffset
	NVAR NxCCD = root:Packages:geometry:PanelValues:NxCCD
	NVAR NyCCD = root:Packages:geometry:PanelValues:NyCCD
	NVAR dpsx = root:Packages:geometry:PanelValues:dpsx
	NVAR dpsy = root:Packages:geometry:PanelValues:dpsy
	NVAR xcent = root:Packages:geometry:PanelValues:xcent
	NVAR ycent = root:Packages:geometry:PanelValues:ycent
	NVAR xbet = root:Packages:geometry:PanelValues:xbet
	NVAR xgam = root:Packages:geometry:PanelValues:xgam
	NVAR dia = root:Packages:geometry:PanelValues:dia
	NVAR H0 = root:Packages:geometry:PanelValues:H0
	NVAR Hyc = root:Packages:geometry:PanelValues:Hyc
	NVAR F = root:Packages:geometry:PanelValues:F
	NVAR axisX = 	root:Packages:geometry:PanelValues:axisX
	NVAR axisY = 	root:Packages:geometry:PanelValues:axisY
	NVAR axisZ = 	root:Packages:geometry:PanelValues:axisZ

	STRUCT microGeometry geo
	geo.NxCCD = 0
	if (strlen(strStruct)>0)
		StructGet/S/B=2 geo, strStruct					// found structure information, load into geo
		Variable/G root:Packages:geometry:PanelValues:dirty=1	// for passed structure
	else
		FillGeometryStructDefault(geo)					//fill the geometry structure with current values
		Variable/G root:Packages:geometry:PanelValues:dirty=0
	endif
	xalfd = geo.xalfd										// set the panel globals from the struct
	xbetd = geo.xbetd
	dd = geo.dd
	ddOffset = geo.ddOffset
	NxCCD = geo.NxCCD
	NyCCD = geo.NyCCD
	dpsx = geo.dpsx
	dpsy = geo.dpsy
	xcent = geo.xcent
	ycent = geo.ycent
	xbet = geo.xbet
	xgam = geo.xgam
	dia = geo.wire.dia
	H0 = geo.wire.H0
	Hyc = geo.wire.Hyc
	F = geo.wire.F
	axisX = geo.wire.axis[0]
	axisY = geo.wire.axis[1]
	axisZ = geo.wire.axis[2]

	SetWindow kwTopWin,userdata(GeoPanelName)=hostWin+"#geoPanel"
	NewPanel/K=1/W=(left,top,left+221,top+603)/HOST=$hostWin
	ModifyPanel frameStyle=0, frameInset=0
	RenameWindow #,geoPanel
	SetDrawLayer UserBack
	SetDrawEnv fname= "Lucida Grande",fsize= 16
	DrawText 30,25,"Detector Position / tilt"
	SetDrawEnv fname= "Lucida Grande",fsize= 16
	DrawText 30,169,"CCD size / resolution"
	SetDrawEnv fname= "Lucida Grande",fsize= 16
	DrawText 30,270,"Incident beam direction"
	SetDrawEnv fname= "Lucida Grande",fsize= 16
	DrawText 30,330,"Scanning wire"
	SetVariable xalfd,pos={78,27},size={92,15},proc=GeoPanelVarChangedProc,title="xalfd¡"
	SetVariable xalfd,help={"alpha tilt of detector (degree)"}
	SetVariable xalfd,limits={-1,1,0},value= root:Packages:geometry:PanelValues:xalfd,bodyWidth= 60
	SetVariable xbetd,pos={71,47},size={99,15},proc=GeoPanelVarChangedProc,title="xbetad¡"
	SetVariable xbetd,help={"beta tilt of detector (degree)"}
	SetVariable xbetd,limits={-1,1,0},value= root:Packages:geometry:PanelValues:xbetd,bodyWidth= 60
	SetVariable dd,pos={69,67},size={101,15},proc=GeoPanelVarChangedProc,title="dd (mm)"
	SetVariable dd,help={"height of detector (mm)"}
	SetVariable dd,limits={0,2000,0},value= root:Packages:geometry:PanelValues:dd,bodyWidth= 60
	SetVariable ddOffset,pos={38,87},size={132,15},proc=GeoPanelVarChangedProc,title="dd Offset (mm)"
	SetVariable ddOffset,help={"offset between motor and dd,  dd = CCDy+ddOffset"}
	SetVariable ddOffset,limits={-100,100,0},value= root:Packages:geometry:PanelValues:ddOffset,bodyWidth= 60
	SetVariable xcent,pos={48,107},size={122,15},proc=GeoPanelVarChangedProc,title="xcent (pixel)",format="%.3f"
	SetVariable xcent,help={"central pixel in CCD X, 1 based & unbinned"}
	SetVariable xcent,limits={0,5000,0},value= root:Packages:geometry:PanelValues:xcent,bodyWidth= 60
	SetVariable ycent,pos={48,127},size={122,15},proc=GeoPanelVarChangedProc,title="ycent (pixel)",format="%.3f"
	SetVariable ycent,help={"central pixel in CCD y, 1 based & unbinned"}
	SetVariable ycent,limits={0,5000,0},value= root:Packages:geometry:PanelValues:ycent,bodyWidth= 60
	SetVariable NxCCD,pos={76,172},size={94,15},proc=GeoPanelVarChangedProc,title="NxCCD"
	SetVariable NxCCD,help={"x size of CCD (pixels)"},format="%d"
	SetVariable NxCCD,limits={1,5000,0},value= root:Packages:geometry:PanelValues:NxCCD,bodyWidth= 60
	SetVariable NyCCD,pos={76,192},size={94,15},proc=GeoPanelVarChangedProc,title="NyCCD"
	SetVariable NyCCD,help={"y size of CCD (pixels)"},format="%d"
	SetVariable NyCCD,limits={1,5000,0},value= root:Packages:geometry:PanelValues:NyCCD,bodyWidth= 60
	SetVariable dpsx,pos={58,211},size={112,15},proc=GeoPanelVarChangedProc,title="dpsx (mm)"
	SetVariable dpsx,help={"x size of CCD (mm)"}
	SetVariable dpsx,limits={0,500,0},value= root:Packages:geometry:PanelValues:dpsx,bodyWidth= 60
	SetVariable dpsy,pos={58,229},size={112,15},proc=GeoPanelVarChangedProc,title="dpsy (mm)"
	SetVariable dpsy,help={"y size of CCD (mm)"}
	SetVariable dpsy,limits={0,500,0},value= root:Packages:geometry:PanelValues:dpsy,bodyWidth= 60
	SetVariable xbet,pos={76,272},size={94,15},proc=GeoPanelVarChangedProc,title="xbeta¡"
	SetVariable xbet,help={"beta tilt of incident beam (degree)"}
	SetVariable xbet,limits={-1,1,0},value= root:Packages:geometry:PanelValues:xbet,bodyWidth= 60
	SetVariable xgam,pos={77,292},size={93,15},proc=GeoPanelVarChangedProc,title="xgam¡"
	SetVariable xgam,help={"gamma tilt of incident beam (degree)"}
	SetVariable xgam,limits={-12,12,0},value= root:Packages:geometry:PanelValues:xgam,bodyWidth= 60
	SetVariable dia,pos={44,332},size={126,15},proc=GeoPanelVarChangedProc,title="wire dia (µm)"
	SetVariable dia,help={"diameter of wire (micron)"}
	SetVariable dia,limits={1,500,0},value= root:Packages:geometry:PanelValues:dia,bodyWidth= 60
	SetVariable H0,pos={45,352},size={125,15},proc=GeoPanelVarChangedProc,title="wire H0 (µm)"
	SetVariable H0,help={"H of wire where it cuts incident beam"}
	SetVariable H0,limits={-50000,50000,0},value= root:Packages:geometry:PanelValues:H0,bodyWidth= 60
	SetVariable Hyc,pos={40,372},size={130,15},proc=GeoPanelVarChangedProc,title="wire Hyc (µm)"
	SetVariable Hyc,help={"H of wire where it cuts vetical ray to (xc,yc)"}
	SetVariable Hyc,limits={-50000,50000,0},value= root:Packages:geometry:PanelValues:Hyc,bodyWidth= 60
	SetVariable wireF,pos={52,392},size={118,15},proc=GeoPanelVarChangedProc,title="wire F (µm)"
	SetVariable wireF,help={"F of wire at H0 & Hyc"}
	SetVariable wireF,limits={-50000,50000,0},value= root:Packages:geometry:PanelValues:F,bodyWidth= 60

//	SetVariable wireAxisX,pos={52,392},size={30,15},proc=GeoPanelVarChangedProc,title="wire F (µm)"
//	SetVariable wireAxisX,help={"direction of wire axis, beam line coords"}
//	SetVariable wireAxisX,limits={1,1,0},value= root:Packages:geometry:PanelValues:axisX,bodyWidth= 60
//	SetVariable wireAxisY,pos={52,392},size={30,15},proc=GeoPanelVarChangedProc
//	SetVariable wireAxisY,help={"direction of wire axis, beam line coords"}
//	SetVariable wireAxisY,limits={1,1,0},value= root:Packages:geometry:PanelValues:axisY,bodyWidth= 60
//	SetVariable wireAxisZ,pos={52,392},size={30,15},proc=GeoPanelVarChangedProc
//	SetVariable wireAxisZ,help={"direction of wire axis, beam line coords"}
//	SetVariable wireAxisZ,limits={1,1,0},value= root:Packages:geometry:PanelValues:axisZ,bodyWidth= 60

	SetVariable wireAxisX,pos={9,407},size={68,15},proc=GeoPanelVarChangedProc,title="axis"
	SetVariable wireAxisX,help={"direction of wire axis, beam line coords"},format="%.6f"
	SetVariable wireAxisX,limits={1,1,0},value= root:Packages:geometry:PanelValues:axisX,bodyWidth= 45
	SetVariable wireAxisY,pos={84,407},size={45,15},proc=GeoPanelVarChangedProc,title=" "
	SetVariable wireAxisY,help={"direction of wire axis, beam line coords"},format="%.6f"
	SetVariable wireAxisY,limits={1,1,0},value= root:Packages:geometry:PanelValues:axisY,bodyWidth= 45
	SetVariable wireAxisZ,pos={138,407},size={45,15},proc=GeoPanelVarChangedProc,title=" "
	SetVariable wireAxisZ,help={"direction of wire axis, beam line coords"},format="%.6f"
	SetVariable wireAxisZ,limits={1,1,0},value= root:Packages:geometry:PanelValues:axisZ,bodyWidth= 45

	Button buttonSaveDefault,pos={35,422},size={150,20},proc=GeometryPanelButtonProc,title="Save as Default"
	Button buttonSaveDefault,help={"saves these values as default in the current data folder"}
	Button buttonSaveLocal,pos={35,447},size={150,20},proc=GeometryPanelButtonProc,title="Save in Current Fldr"
	Button buttonSaveLocal,help={"saves these values as default in the current data folder"}
	Button buttonEPICS,pos={35,472},size={150,20},proc=GeometryPanelButtonProc,title="Load from EPICS"
	Button buttonEPICS,help={"Fill the values in this panel from EPICS"}
//	Button buttonFillFromFile,pos={35,497},size={150,20},proc=GeometryPanelButtonProc,title="Load from a file"
//	Button buttonFillFromFile,help={"Fill the values in this panel from a file"}
	Button buttonFillFromFile,pos={35,497},size={105,20},proc=GeometryPanelButtonProc,title="Load from file"
	Button buttonFillFromFile,help={"Fill the values in this panel from a file"}
	Button buttonFillFromWeb,pos={145,497},size={40,20},proc=GeometryPanelButtonProc,title="web"
	Button buttonFillFromWeb,help={"Fill the values in this panel from the web"}
	Button buttonWriteToFile,pos={35,522},size={150,20},proc=GeometryPanelButtonProc,title="Write to a file"
	Button buttonWriteToFile,help={"Write values in this panel to a file"}
	Button buttonPrintGeo,pos={35,547},size={150,20},proc=GeometryPanelButtonProc,title="Print Geo to History"
	Button buttonPrintGeo,help={"print geometry values to the history"}
	Button buttonCancel,pos={35,572},size={150,20},proc=GeometryPanelButtonProc,title="Cancel"
	Button buttonCancel,help={"Close this panel, do not save any changes"}
	Button buttonSetdd,pos={171,64},size={45,18},proc=GeometryPanelButtonProc,title="image"
	Button buttonSetdd,help={"set dd using CCDy from top image"}
	return "#geoPanel"
End
//
// this is called by all of the varables, and it just sets the dirty flag for the panel
Function GeoPanelVarChangedProc(ctrlName,varNum,varStr,varName) : SetVariableControl
	String ctrlName
	Variable varNum
	String varStr
	String varName
	NVAR dirty = root:Packages:geometry:PanelValues:dirty
	dirty = 1
End
//
Function SetGeoPanelHook(s)
	STRUCT WMWinHookStruct &s
	if (s.eventCode !=2)							// only process kill events
		return 0
	endif
	Variable dirty = NumVarOrDefault("root:Packages:geometry:PanelValues:dirty",1)
	if (!dirty)
		return 0								// nothing was changed, do not ask any questions
	endif

	// values were changed in the panel and not saved, ask the user what to do.
	Variable action=1
	Prompt action, "choose an action, do not press cancel",Popup,"Save in default data folder;Save in current data folder;just close window, do not save these values"
	do
		DoPrompt/Help="3D-Xray Diffraction[geometry]" "changes were made, but not saved",action
	while(V_flag)
	// action = V_flag ? 3 : action				// if cancel was pushed, do not save, treat it like a "just quit"
	GeometryPanelButtonProc(StringFromList(action-1,"buttonSaveDefault;buttonSaveLocal;buttonCancel"))
End
//
Function GeometryPanelButtonProc(ctrlName) : ButtonControl
	String ctrlName

	STRUCT microGeometry geo
	NVAR xalfd = root:Packages:geometry:PanelValues:xalfd		// connect to the globals
	NVAR xbetd = root:Packages:geometry:PanelValues:xbetd
	NVAR dd = root:Packages:geometry:PanelValues:dd
	NVAR ddOffset = root:Packages:geometry:PanelValues:ddOffset
	NVAR NxCCD = root:Packages:geometry:PanelValues:NxCCD
	NVAR NyCCD = root:Packages:geometry:PanelValues:NyCCD
	NVAR dpsx = root:Packages:geometry:PanelValues:dpsx
	NVAR dpsy = root:Packages:geometry:PanelValues:dpsy
	NVAR xcent = root:Packages:geometry:PanelValues:xcent
	NVAR ycent = root:Packages:geometry:PanelValues:ycent
	NVAR xbet = root:Packages:geometry:PanelValues:xbet
	NVAR xgam = root:Packages:geometry:PanelValues:xgam
	NVAR dia = root:Packages:geometry:PanelValues:dia
	NVAR H0 = root:Packages:geometry:PanelValues:H0
	NVAR Hyc = root:Packages:geometry:PanelValues:Hyc
	NVAR F = root:Packages:geometry:PanelValues:F
	NVAR axisX = root:Packages:geometry:PanelValues:axisX
	NVAR axisY = root:Packages:geometry:PanelValues:axisY
	NVAR axisZ = root:Packages:geometry:PanelValues:axisZ
	NVAR dirty = root:Packages:geometry:PanelValues:dirty

	if (stringmatch(ctrlName,"buttonCancel"))				// close window, do not save numbers
		String win = StringFromLIst(0,WinList("*",";","WIN:64"))
		SetWindow $win,hook(close)= $""
		DoWindow/K $win
		dirty = 0
		return 0
	elseif (stringmatch(ctrlName,"buttonSaveDefault"))
		String/G root:Packages:geometry:geoStructStr=""	// use default location
		SVAR geoStructStr = root:Packages:geometry:geoStructStr
	elseif (stringmatch(ctrlName,"buttonSaveLocal"))
		String/G geoStructStr=""							// save it in local folder
		SVAR geoStructStr = geoStructStr
	elseif (stringmatch(ctrlName,"buttonEPICS"))
		if (FillGeometryStructFromEpics(geo))
			return 1
		endif
		xalfd = geo.xalfd										// set the panel globals from the struct
		xbetd = geo.xbetd
		dd = geo.dd/1000
		ddOffset = geo.ddOffset/1000
		NxCCD = geo.NxCCD
		NyCCD = geo.NyCCD
		dpsx = geo.dpsx
		dpsy = geo.dpsy
		xcent = geo.xcent
		ycent = geo.ycent
		xbet = geo.xbet
		xgam = geo.xgam
		dia = geo.wire.dia
		H0 = geo.wire.H0
		Hyc = geo.wire.Hyc
		F = geo.wire.F
		axisX = geo.wire.axis[0]
		axisY = geo.wire.axis[1]
		axisZ = geo.wire.axis[2]
		dirty = 1
		return 0
	elseif (stringmatch(ctrlName,"buttonFillFromWeb"))
		Variable fnum
		Open/D/M="Select raw, NOT reconstructed SPE image"/P=home/R/T=".SPE" fnum
		if (strlen(S_fileName)<1)
			return 1
		endif
		GetFileFolderInfo/Q/Z S_fileName
		if (V_Flag)
			return 1
		endif
		if (GeoFromWeb(V_modificationDate,geo))
			return 1
		endif
		xalfd = geo.xalfd										// set the panel globals from the struct
		xbetd = geo.xbetd
		dd = geo.dd
		ddOffset = geo.ddOffset
		NxCCD = geo.NxCCD
		NyCCD = geo.NyCCD
		dpsx = geo.dpsx
		dpsy = geo.dpsy
		xcent = geo.xcent
		ycent = geo.ycent
		xbet = geo.xbet
		xgam = geo.xgam
		dia = geo.wire.dia
		H0 = geo.wire.H0
		Hyc = geo.wire.Hyc
		F = geo.wire.F
		axisX = geo.wire.axis[0]
		axisY = geo.wire.axis[1]
		axisZ = geo.wire.axis[2]
		dirty = 1
		return 0
	elseif (stringmatch(ctrlName,"buttonFillFromFile"))
		if (ReadGeoFromKeyFile("","home",geo))
			return 1
		endif
		xalfd = geo.xalfd										// set the panel globals from the struct
		xbetd = geo.xbetd
		dd = geo.dd
		ddOffset = geo.ddOffset
		NxCCD = geo.NxCCD
		NyCCD = geo.NyCCD
		dpsx = geo.dpsx
		dpsy = geo.dpsy
		xcent = geo.xcent
		ycent = geo.ycent
		xbet = geo.xbet
		xgam = geo.xgam
		dia = geo.wire.dia
		H0 = geo.wire.H0
		Hyc = geo.wire.Hyc
		F = geo.wire.F
		axisX = geo.wire.axis[0]
		axisY = geo.wire.axis[1]
		axisZ = geo.wire.axis[2]
		dirty = 1
		return 0
	elseif (stringmatch(ctrlName,"buttonWriteToFile"))
		geo.xbet = xbet										// set all the geo values from global values
		geo.xgam = xgam
		geo.xalfd = xalfd
		geo.xbetd = xbetd
		geo.dd = dd
		geo.ddOffset = ddOffset
		geo.NxCCD = NxCCD
		geo.NyCCD = NyCCD
		geo.dpsx = dpsx
		geo.dpsy = dpsy
		geo.xcent = xcent
		geo.ycent = ycent
		geo.wire.dia = dia
		geo.wire.F = F
		geo.wire.H0 = H0
		geo.wire.Hyc = Hyc
		geo.wire.axis[0] = axisX
		geo.wire.axis[1] = axisY
		geo.wire.axis[2] = axisZ
		WriteGeoToKeyFile("","home",geo,"")
		return 0
	elseif (stringmatch(ctrlName,"buttonPrintGeo"))		// print shown geometry parameters to the history
		geo.xbet = xbet										// set all the geo values from global values
		geo.xgam = xgam
		geo.xalfd = xalfd
		geo.xbetd = xbetd
		geo.dd = dd
		geo.ddOffset = ddOffset
		geo.NxCCD = NxCCD
		geo.NyCCD = NyCCD
		geo.dpsx = dpsx
		geo.dpsy = dpsy
		geo.xcent = xcent
		geo.ycent = ycent
		geo.wire.dia = dia
		geo.wire.F = F
		geo.wire.H0 = H0
		geo.wire.Hyc = Hyc
		geo.wire.axis[0] = axisX
		geo.wire.axis[1] = axisY
		geo.wire.axis[2] = axisZ
		GeometryUpdateCalc(geo)							// calculate other values
		printGeometry(geo)									// print them all out
		return 0
	elseif (stringmatch(ctrlName,"buttonSetdd"))			// set dd using values from top image
		Wave image=ImageNameToWaveRef("",StringFromList(0,ImageNameList("",";")))
		if (!WaveExists(image))
			DoAlert 0,"image cannot be found for top graph"
			return 1
		endif
		Variable CCDy = NumberByKey("CCDy",note(image),"=")
		if (numtype(CCDy))
			DoAlert 0,"unable to read CCD height from image '"+NameOfWave(image)+"'"
			return 1
		endif
		dd = ddOffset+CCDy
		dirty = 1
		return 0
	else
		Abort "''GeometryPanelButtonProc' called from unsupported control"
	endif

	// store the global values from panel into the 'geoStructStr''
	geo.xalfd = xalfd
	geo.xbetd = xbetd
	geo.dd = dd
	geo.ddOffset = ddOffset
	geo.NxCCD = NxCCD
	geo.NyCCD = NyCCD
	geo.dpsx = dpsx
	geo.dpsy = dpsy
	geo.xcent = xcent
	geo.ycent = ycent
	geo.xbet = xbet
	geo.xgam = xgam
	geo.wire.dia = dia
	geo.wire.H0 = H0
	geo.wire.Hyc = Hyc
	geo.wire.F = F
	geo.wire.axis[0] = axisX
	geo.wire.axis[1] = axisY
	geo.wire.axis[2] = axisZ
	GeometryUpdateCalc(geo)
	StructPut/S/B=2 geo, geoStructStr
	dirty = 0
	SavePackagePreferences/FLSH=1 "microGeo","microGeoPrefs",0,geo
End
//
Function FillGeometryStructFromEpics(geo)		//fill the geometry structure from EPICS (uses caget)
	STRUCT microGeometry &geo

	String geoList = "xalfd;xbetd;dd;ddOffset;NxCCD;NyCCD;dpsx;dpsy;xcent;ycent;xbet;xgam;wireDia;wireH0;wireH1;wireF0"
	String item,pvList=""
	Variable i

	// make list of pv names for each item in geoList (items in geo that have EPICS pv's)
	for (i=0,item=StringFromList(0,geoList); strlen(item); i+=1,item=StringFromList(i,geoList))
		pvList += EPICS_PREFIX+item+";"
	endfor

	// use FUNCREF incase this computer does not have the function get_mult_PV()
	String funcSpec = SelectString(exists("get_mult_PV")==6 , "protoEPICS_geo", "get_mult_PV")
	FUNCREF protoEPICS_geo  func=$funcSpec
	String epicsOut = func(pvList)				// get values from EPICS
	if (strlen(epicsOut)<2)
		return 1
	endif
	// turn result of get_mult_PV() into a keyword=value list
	Variable n=strlen(EPICS_PREFIX)
	pvList = ""
	for (i=0,item=StringFromList(0,epicsOut); strlen(item); i+=1,item=StringFromList(i,epicsOut))
		if (strsearch(item,EPICS_PREFIX,0)==0)
			item = item[n,Inf]
		endif
		pvList += item+";"
	endfor
	geo.xbet = NumberByKey("xbet",pvList,"=")				// set all the geo values from keword list
	geo.xgam = NumberByKey("xgam",pvList,"=")
	geo.xalfd = NumberByKey("xalfd",pvList,"=")
	geo.xbetd = NumberByKey("xbetd",pvList,"=")
	geo.dd = NumberByKey("dd",pvList,"=")*1000				// convert mm to µm
	geo.ddOffset = NumberByKey("ddOffset",pvList,"=")*1000	// convert mm to µm
	geo.NxCCD = NumberByKey("NxCCD",pvList,"=")
	geo.NyCCD = NumberByKey("NyCCD",pvList,"=")
	geo.dpsx = NumberByKey("dpsx",pvList,"=")				// convert mm to µm
	geo.dpsy = NumberByKey("dpsy",pvList,"=")				// convert mm to µm
	geo.xcent = NumberByKey("xcent",pvList,"=")
	geo.ycent = NumberByKey("ycent",pvList,"=")

	geo.wire.dia = NumberByKey("wireDia",pvList,"=")
	geo.wire.F = NumberByKey("wireF0",pvList,"=")
	geo.wire.H0 = NumberByKey("wireH0",pvList,"=")
	geo.wire.Hyc = NumberByKey("wireH1",pvList,"=")
	geo.wire.axis[0] = NumberByKey("wireaxisX",pvList,"=")
	geo.wire.axis[1] = NumberByKey("wireaxisY",pvList,"=")
	geo.wire.axis[2] = NumberByKey("wireaxisZ",pvList,"=")
	GeometryUpdateCalc(geo)
	return 0
End
//
Function/S protoEPICS_geo(str)			// proto function for get_mult_PV(), also contains default values for testing
	String str
	DoAlert 0, "EPICS not available, goemetry unchanged"
	return ""
//	if (!stringmatch(StringFromList(0,GetRTStackInfo(0)),"GeometryPanelButtonProc"))
//		print "EPICS not available, Using default geo values"
//	endif
//	DoAlert 1, "EPICS not available, Use default geo values?"
//	if (V_flag!=1)
//		return ""
//	endif
//	str = "yum:geometry:xalfd=0.15605;yum:geometry:xbetd=0.19635;yum:geometry:dd=40.0143;yum:geometry:ddOffset=0.0143;"
//	str += "yum:geometry:NxCCD=2084;yum:geometry:NyCCD=2084;yum:geometry:dpsx=50.01;yum:geometry:dpsy=50.168;"
//	str += "yum:geometry:xcent=1089.7;yum:geometry:ycent=942.901;yum:geometry:xbet=0.30943;yum:geometry:xgam=0.75192;"
//	str += "yum:geometry:wireDia=52;yum:geometry:wireH0=-5305.9;yum:geometry:wireH1=-4889.03;yum:geometry:wireF0=3830; "
//	return str
End
//
//Function testFillGeometryStructFromEpics()
//	STRUCT microGeometry geo
//	FillGeometryStructFromEpics(geo)
//	printGeometry(geo)
//	String/G geoSructStr
//	StructPut/S/B=2 geo,geoSructStr
//End





// write the geometry structure to a keyword file
Function WriteGeoToKeyFile(fileName,path,geo,fileNote)
	String fileName									// full path name to the file
	String path										// name of an Igor path to use
	STRUCT microGeometry &geo					// the structure to fill from the file
	String fileNote									// a note about this file
	GeometryUpdateCalc(geo)						// calculate other values
	Variable f										// file id
	Open/C="R*ch"/M="new geometry parameters file"/P=$path/T="TEXT"/Z=2 f as fileName
	if (V_flag)
		DoAlert 0, "nothing written to file"
		return 1
	endif
	if( (strlen(fileName)+strlen(fileNote))==0 )
		Prompt fileNote, "add a descriptive note to this file?"
		DoPrompt "descriptive note",fileNote
		if (V_flag)
			fileNote=""
		endif
	endif
	fileName = S_fileName
	fprintf f,"$filetype		geometryFile\n"
//	fprintf f,"$geometryFile\n"
	fprintf f, "$dateWritten	%s\n", date()
	fprintf f, "$timeWritten	%s (%g)\n", Secs2Time(DateTime,3,1),date2secs(-1,-1,-1)/3600
	fprintf f, "$EPOCH			%.0f		// seconds from midnight January 1, 1904\n",DateTime
	if (strlen(fileNote))
		fprintf f,"$fileNote		%s\n",ReplaceString("\r",fileNote,"\\r")
	endif
	fprintf f,"$xalfd			%.5f			// CCD tilt (degrees)\n",geo.xalfd
	fprintf f,"$xbetd			%.5f\n",geo.xbetd
	fprintf f,"$dd				%.4f			// height of CCD (µm)\n",geo.dd
	fprintf f,"$ddOffset		%.4f			// dd = CCDy + ddOffset (µm)\n",geo.ddOffset
	fprintf f,"$NxCCD			%d			// number of un-binned pixels in CCD, also called xdim, ydim\n",geo.NxCCD
	fprintf f,"$NyCCD			%d\n",geo.NyCCD
	fprintf f,"$dpsx			%.3f			// size of CCD (mm)\n",geo.dpsx
	fprintf f,"$dpsy			%.3f\n",geo.dpsy
	fprintf f,"$xcent			%.3f		// center of CCD in un-binned 1 based pixels\n",geo.xcent
	fprintf f,"$ycent			%.3f\n",geo.ycent
	fprintf f,"$xbet			%.5f			// incident beam tilt (degrees)\n",geo.xbet
	fprintf f,"$xgam			%.5f\n",geo.xgam
	fprintf f,"$dia			%.2f			// diameter of wire (µm)\n",geo.wire.dia
	fprintf f,"$H0				%.2f		// H of wire centered on the direct beam (µm)\n",geo.wire.H0
	fprintf f,"$Hyc			%.2f		// H of wire centered on the vertical over the Si (µm)\n",geo.wire.Hyc
	fprintf f,"$F				%.2f			// F of the wire for the Si calibration wire scan (µm)\n",geo.wire.F
	fprintf f,"$wireAxis		{%.6f,%.6f,%.6f} // unit vector along wire axis\n",geo.wire.axis[0],geo.wire.axis[1],geo.wire.axis[2]
	Close f
	String str
	sprintf str, "finished wrting geometry to file '%s'",fileName
	DoAlert 0,str
	return 0
End


// read in the geometry structure from a keyword file
Function ReadGeoFromKeyFile(fileName,path,geo)
	String fileName									// full path name to the file
	String path										// name of an Igor path to use
	STRUCT microGeometry &geo					// the structure to fill from the file

	String list = keyStrFromFile(fileName,"geometryFile",path)// read in all of the tagged values into a keyword list
	if (strlen(list)<1)
		return 1
	endif

	Variable xalfd, xbetd, dd, ddOffset, NxCCD, NyCCD, dpsx, dpsy, xcent, ycent, xbet, xgam, dia, H0, Hyc, F,Ax,Ay,Az
	xalfd = NumberByKey("xalfd",list,"=")
	xbetd = NumberByKey("xbetd",list,"=")
	dd = NumberByKey("dd",list,"=")
	ddOffset = NumberByKey("ddOffset",list,"=")
	NxCCD = NumberByKey("NxCCD",list,"=")
	NyCCD = NumberByKey("NyCCD",list,"=")
	dpsx = NumberByKey("dpsx",list,"=")
	dpsy = NumberByKey("dpsy",list,"=")
	xcent = NumberByKey("xcent",list,"=")
	ycent = NumberByKey("ycent",list,"=")
	xbet = NumberByKey("xbet",list,"=")
	xgam = NumberByKey("xgam",list,"=")
	dia = NumberByKey("dia",list,"=")
	H0 = NumberByKey("H0",list,"=")
	Hyc = NumberByKey("Hyc",list,"=")
	F = NumberByKey("F",list,"=")
	sscanf StringByKey("wireAxis",list,"="),"{%g,%g,%g}",Ax,Ay,Az
	if (V_flag!=3)
		Ax=NaN  ;  Ay=NaN  ;  Az=NaN 
	endif

	geo.xbet = numtype(xbet) ? geo.xbet : xbet		// set all the geo if they are present
	geo.xgam = numtype(xgam) ? geo.xgam : xgam
	geo.xalfd = numtype(xalfd) ? geo.xalfd : xalfd
	geo.xbetd = numtype(xbetd) ? geo.xbetd : xbetd
	geo.dd = numtype(dd) ? geo.dd : dd
	geo.ddOffset = numtype(ddOffset) ? geo.ddOffset : ddOffset
	geo.NxCCD = numtype(NxCCD) ? geo.NxCCD : NxCCD
	geo.NyCCD = numtype(NyCCD) ? geo.NyCCD : NyCCD
	geo.dpsx = numtype(dpsx) ? geo.dpsx : dpsx
	geo.dpsy = numtype(dpsy) ? geo.dpsy : dpsy
	geo.xcent = numtype(xcent) ? geo.xcent : xcent
	geo.ycent = numtype(ycent) ? geo.ycent : ycent
	geo.wire.dia = numtype(dia) ? geo.wire.dia : dia
	geo.wire.F = numtype(F) ? geo.wire.F : F
	geo.wire.H0 = numtype(H0) ? geo.wire.H0 : H0
	geo.wire.Hyc = numtype(Hyc) ? geo.wire.Hyc : Hyc
	if (numtype(Ax+Ay+Az)==0)
		geo.wire.axis[0] = Ax
		geo.wire.axis[1] = Ay
		geo.wire.axis[2] = Az
	endif
	GeometryUpdateCalc(geo)						// calculate other values
	printf "Loaded geometry information from '%s'\r",StringByKey("keyStrFileName",list,"=")
	return 0
End
//
//Function testReadGeoFromKeyFile()
//	STRUCT microGeometry geo
//	Variable i = ReadGeoFromKeyFile("geometry.txt",geo)
//	printGeometry(geo)
//End
//
Static Function/S keyStrFromFile(fname,ftype,path)// read in all of the tagged values from a file into a keyword list
	String fname										// full path name to file with tagged geometry values
	String ftype										// the required file identifier, included as a tag (ftype is optional)
	String path										// name of Igor path

	Variable refNum
	Open/M="file containing tagged values"/P=$path/R/Z=2 refNum as fname
	if (strlen(S_fileName)<1 || !refNum)
		return ""
	endif

	Variable i
	String line
	if (strlen(ftype))								// an ftype is present, so check the file
		Variable OK=0
		FReadLine refNum, line
		line[strlen(line)-1,Inf] = " "				// change the trailing term char to a space
		if (strsearch(line,"$filetype",0)==0)		// starts with "$filetype"
			OK = char2num(line[9])<=32			// $filetype ends with whitespace
			i = strsearch(line,"//",0)				// search for comment identifier
			i = (i<0) ? strlen(line) : i				// points to start of comment id (or end of string)
			line = line[0,i-1]						// trim off comment
			line = ReplaceString("\t",line,";")		// change all tabs, commas, and spaces to semi-colons
			line = ReplaceString(",",line,";")
			line = ReplaceString(" ",line,";")
			OK = OK && WhichListItem(ftype,line)>=0
		else
			ftype = "$"+ftype
			OK = (strsearch(line,ftype,0)==0 && char2num(line[strlen(ftype)])<=32)
		endif
		if (!OK)			// if (strsearch(line,ftype,0) || char2num(line[strlen(ftype)])>32)
			printf "The file must start with '%s', but the first line is '%s'\r",ftype,line
			Close refNum
			return ""
		endif
	endif
//	print "loading geometry parameters from ",S_fileName
	String tagName,value,list = ReplaceStringByKey("keyStrFileName","",S_fileName,"=")
	Variable dollar = char2num("$")
	do
		FReadLine refNum, line
		if (char2num(line)==dollar)				// if it starts with a $, probable tag so process it
			i = strsearch(line,"//",0)				// strip off comments
			if (i>=0)
				line = line[0,i-1]
			endif
			for (i=0;char2num(line[i+1])>32;i+=1)// find end of tag, it ends with a space or lower
			endfor
			tagName = line[1,i]
			if (strlen(tagName)<1)
				continue
			endif

			// check if tag is already in list, only take the first instance of a tag, not the last
			if (keyInList(tagName,list,"=",""))			// check if this key is already in list
				continue
			endif

			for (i=i+1;char2num(line[i])<=32;i+=1)	// find first non-white space
			endfor
			value = line[i,Inf]							// value associated with tagName

			for (i=strlen(value)-1;char2num(value[i])<=32 && i>0;i-=1)	// strip off trailing whitespace
			endfor
			value = ReplaceString(";",value[0,i],":")
			list = ReplaceStringByKey(tagName,list,value,"=")
		endif
	while (strlen(line))							// go until EOF
	Close refNum
	return list
End
//
//	$geometryFile
//	$NxCCD		2084					// number of un-binned pixels in CCD, also called xdim, ydim
//	$NyCCD		2084
//	$dpsx		50.01						// size of CCD (mm)
//	$dpsy		50.168
//	$xcent		1089.70					// center of CCD in un-binned pixels
//	$ycent		942.901
//	$dd			40.0143					// height of CCD (µm)
//	$ddOffset	0.0143					// dd = CCDy + ddOffset (µm)
//	$xbet		0.30943					// incident beam tilt (degrees)
//	$xgam		0.75192
//	$xalfd		0.15605					// CCD tilt (degrees)
//	$xbetd		0.19635
//	//	for the wire:
//	$F			3830						// F of the wire for the Si calibration wire scan (µm)
//	$H0		-5305.9					// H of wire centered on the direct beam (µm)
//	$Hyc		-4889.03					// H of wire centered on the vertical over the Si (µm)
//	$X			0							// X of wire for Si (µm)
//	$dia		52							// diameter of wire (µm)
//	$wireAxis	{1, 0, 0.008727}		// direction of wire axis, should be close to {1,0,0}


// with an Igor epoch, retrieve the geo data form the web server
Function GeoFromWeb(epoch,geo)
	Variable epoch					// Igor epoch, also value of V_modificationDate from GetFileFolderInfo
	STRUCT microGeometry &geo

	String ddate, ttime				// ddate="2007-03-16",ttime="14:32:55"
	Variable im,id,iy
	sscanf Secs2Date(epoch,-1),"%d/%d/%d",id,im,iy
	sprintf ddate "%04d-%02d-%02d",iy,im,id
	ttime = Secs2Time(epoch,3)

	String cmd
	sprintf cmd "do shell script \"curl -s -f -m 20 http://www.uni.aps.anl.gov/34ide/index.php\\\?tag=geo\\\&date=%sT%s\"",ddate,ttime
	//	printf "cmd = %s\r\r",cmd

	ExecuteScriptText cmd
	// printf "%s", S_value  ;  print " \r\r"
	if (strsearch(S_value,"ERROR",0)==1)
		return 1
	endif

	Variable N = strlen(S_value)
	if (strsearch(S_value, "\"",0) == 0)
		S_value = S_value[1,N-1]
		N -= 1
	endif
	if (strsearch(S_value, "\"",N-1,1) == N-1)
		S_value = S_value[0,N-2]
	endif
	String list = keyStrFromBuffer(S_value)

	Variable xalfd, xbetd, dd, ddOffset, NxCCD, NyCCD, dpsx, dpsy, xcent, ycent, xbet, xgam, dia, H0, Hyc, F,Ax,Ay,Az
	xalfd = NumberByKey("xalfd",list,"=")
	xbetd = NumberByKey("xbetd",list,"=")
	dd = NumberByKey("dd",list,"=")
	ddOffset = NumberByKey("ddOffset",list,"=")
	NxCCD = NumberByKey("NxCCD",list,"=")
	NyCCD = NumberByKey("NyCCD",list,"=")
	dpsx = NumberByKey("dpsx",list,"=")
	dpsy = NumberByKey("dpsy",list,"=")
	xcent = NumberByKey("xcent",list,"=")
	ycent = NumberByKey("ycent",list,"=")
	xbet = NumberByKey("xbet",list,"=")
	xgam = NumberByKey("xgam",list,"=")
	dia = NumberByKey("dia",list,"=")
	H0 = NumberByKey("H0",list,"=")
	Hyc = NumberByKey("Hyc",list,"=")
	F = NumberByKey("F",list,"=")
	sscanf StringByKey("wireAxis",list,"="),"{%g,%g,%g}",Ax,Ay,Az
	if (V_flag!=3)
		Ax=NaN  ;  Ay=NaN  ;  Az=NaN 
	endif

	geo.xbet = numtype(xbet) ? geo.xbet : xbet		// set all the geo if they are present
	geo.xgam = numtype(xgam) ? geo.xgam : xgam
	geo.xalfd = numtype(xalfd) ? geo.xalfd : xalfd
	geo.xbetd = numtype(xbetd) ? geo.xbetd : xbetd
	geo.dd = numtype(dd) ? geo.dd : dd
	geo.ddOffset = numtype(ddOffset) ? geo.ddOffset : ddOffset
	geo.NxCCD = numtype(NxCCD) ? geo.NxCCD : NxCCD
	geo.NyCCD = numtype(NyCCD) ? geo.NyCCD : NyCCD
	geo.dpsx = numtype(dpsx) ? geo.dpsx : dpsx
	geo.dpsy = numtype(dpsy) ? geo.dpsy : dpsy
	geo.xcent = numtype(xcent) ? geo.xcent : xcent
	geo.ycent = numtype(ycent) ? geo.ycent : ycent
	geo.wire.dia = numtype(dia) ? geo.wire.dia : dia
	geo.wire.F = numtype(F) ? geo.wire.F : F
	geo.wire.H0 = numtype(H0) ? geo.wire.H0 : H0
	geo.wire.Hyc = numtype(Hyc) ? geo.wire.Hyc : Hyc
	if (numtype(Ax+Ay+Az)==0)
		geo.wire.axis[0] = Ax
		geo.wire.axis[1] = Ay
		geo.wire.axis[2] = Az
	endif
	GeometryUpdateCalc(geo)						// calculate other values
	if (ItemsInList(GetRTStackInfo(0))<3)		// print everything if run from command line
		printf "Loaded geometry information from web on  %s,  %s",Secs2Date(epoch,2),Secs2Time(epoch,1)
		printf "     this geometry created on   %s  %s\r",StringByKey("dateWritten",list,"="),StringByKey("timeWritten",list,"=")
	endif
	return 0
End
//
Static Function/S keyStrFromBuffer(buf)// read in all of the tagged values from a file into a keyword list
	String buf

	buf = ReplaceString("\r",buf,"\n")
	buf = ReplaceString("\t",buf," ")
	buf += "\n"

	Variable i, i0=0,i1
	String line
	String tagName,value,list = ""
	Variable dollar = char2num("$")
	do
		i1 = strsearch(buf, "\n",i0)
		line = buf[i0,i1-1]
		if (char2num(line)==dollar)				// if it starts with a $, probable tag so process it
			i = strsearch(line,"//",0)				// strip off comments
			if (i>=0)
				line = line[0,i-1]
			endif
			for (i=0;char2num(line[i+1])>32;i+=1)// find end of tag, it ends with a space or lower
			endfor
			tagName = line[1,i]
			if (strlen(tagName)<1)
				continue
			endif

			// check if tag is already in list, only take the first instance of a tag, not the last
			if (keyInList(tagName,list,"=",""))			// check if this key is already in list
				continue
			endif

			for (i=i+1;char2num(line[i])<=32;i+=1)	// find first non-white space
			endfor
			value = line[i,Inf]							// value associated with tagName

			for (i=strlen(value)-1;char2num(value[i])<=32 && i>0;i-=1)	// strip off trailing whitespace
			endfor
			value = ReplaceString(";",value[0,i],":")
			list = ReplaceStringByKey(tagName,list,value,"=")
		endif
		i0 = i1 + 1
	while (strlen(line))
	return list
End



Function FillGeometryStructDefault(geo)					//fill the geometry structure with test values
	STRUCT microGeometry &geo								// returns 0 if something set, 0 is nothing done

	String strStruct=StrVarOrDefault(":geoStructStr","")// set to values in current directory
	if (strlen(strStruct)<1)
		strStruct=StrVarOrDefault("root:Packages:geometry:geoStructStr","")	// try the default values
	endif
	if (strlen(strStruct)>1)
		StructGet/S/B=2 geo, strStruct						// found structure information, load into geo
	else
		LoadPackagePreferences/MIS=1 "microGeo","microGeoPrefs",0,geo
		if (V_flag)
			return 1											// did nothing, nothing found
		endif
	endif
	GeometryUpdateCalc(geo)
	return 0
End

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



Function setUseDistortion(use)
	Variable use
	Init_microGeo()
	NVAR useDistortion=root:Packages:geometry:useDistortion
	Variable before = useDistortion							// value on entry
	if (numtype(use)==0)									//  valid input, set the global
		useDistortion = use
	else														// need to prompt user
		Prompt use, "Use Distortion Correction?", popup, "NO Distortion Correction;Use Distortion Correction"
		use = useDistortion+1								// default value
		DoPrompt/Help="3D-Xray Diffraction[Distortion Correction]" "Distortion Correction", use
		if (!V_flag)
			useDistortion = use - 1
		endif
	endif
	if (useDistortion && exists("root:Packages:geometry:xymap")!=1)
		LoadStandardDistortion()
	endif
	if (ItemsInList(GetRTStackInfo(0))<2)
		if (before == useDistortion)
			print SelectString(useDistortion,"Distortion Correction remains OFF","Continuing to Use Distortion Correction")
		else
			printf "¥¥Distortion Correction changed from %s to %s\r",SelectString(before,"OFF","ON"),SelectString(useDistortion,"OFF","ON")
		endif
	endif
	return useDistortion
End


//	xymap,Nx,Ny,cornerX0,cornerY0,cornerX1,cornerY1
Function/S loadCCDCorrectionTable(mName)
	String mName						// name for xymap file
	String CCDxyMapfile = "CCD_distorMay03_corr.dat"
	mName = SelectString(strlen(mName)>0,"xymap",CleanupName(mName,0))

	if (ItemsInList(GetRTStackInfo(0))<2)
		print "Loading the correction table ..."
	endif
	Variable refNum
	Open/M="distortion table"/P=home/R/T="TEXT"/Z=2 refNum as CCDxyMapfile
	if (V_flag)
		return ""
	endif
	String line
	FReadLine refNum, line

	Variable Nx,Ny,cornerX0,cornerY0,cornerX1,cornerY1
	sscanf line, "%d %d %d %d %d %d",Nx,Ny,cornerX0,cornerY0,cornerX1,cornerY1
	if (V_flag!=6)
		Close refNum
		return ""
	endif

	String noteStr=""
	noteStr = ReplaceNumberByKey("cornerX0",noteStr,cornerX0,"=")
	noteStr = ReplaceNumberByKey("cornerY0",noteStr,cornerY0,"=")
	noteStr = ReplaceNumberByKey("cornerX1",noteStr,cornerX1,"=")
	noteStr = ReplaceNumberByKey("cornerY1",noteStr,cornerY1,"=")
	noteStr = ReplaceStringByKey("CCDname",noteStr,"White2084x2084","=")

	Make/N=(Nx,Ny,4)/O $mName
	Wave xymap = $mName
	Note/K xymap
	Note xymap, noteStr

	Variable i,j,x1,y1,x2,y2
	for (j=0;j<Ny;j+=1)
		for (i=0;i<Ny;i+=1)
			FReadLine refNum, line
			sscanf line, "%g %g %g %g",x1,y1,x2,y2
			if (V_flag!=4)
				Close refNum
				return ""
			endif
			xymap[i][j][0] = x1
			xymap[i][j][1] = y1
			xymap[i][j][2] = x2
			xymap[i][j][3] = y2
		endfor
	endfor
	Close refNum
	return GetWavesDataFolder(xymap,2)
End
//
Function LoadStandardDistortion()
	String str= FunctionPath("LoadStandardDistortion")
	if (strlen(str)<2)
		return 1
	endif
	String fileName = ParseFilePath(1, str, ":", 1, 0)+"xymap.itx"
	LoadWave/O/T fileName
	if (V_flag!=1)
		DoAlert 0, "Unable to load xymap.itx"
		return 1
	elseif (!stringmatch(StringFromList(0,S_waveNames),"xymap"))
		DoAlert 0,"Loaded something other than xyamp, see history"
		return 1
	endif

	// move to root:Packages:geometry if not already there
	if (stringmatch(GetDataFolder(1),"root:Packages:geometry:"))
		return 0
	endif
	Duplicate/O xymap root:Packages:geometry:xymap
	KillWaves xymap
	printf "¥¥loaded standard distortion from '%s' into file root:Packages:geometry:xymap\r",fileName
	return 0
End


Function MakeDistortionMap()
	Variable N=20
	Make/N=(N,N)/O totalDistortion
	Make/N=(N,N)/O/C xyDistortion
	SetScale/I x 0,2083,"pixels", totalDistortion
	SetScale/I y 0,2083,"pixels", totalDistortion
	SetScale d 0,0,"pixels", totalDistortion
	Variable x0=DimOffset(totalDistortion,0), dx=DimDelta(totalDistortion,0)
	Variable y0=DimOffset(totalDistortion,1), dy=DimDelta(totalDistortion,1)

	Wave xymap=root:Packages:geometry:xymap
	Variable px,py, i,j
	for (j=0;j<N;j+=1)
		py = y0 + j*dy
		for (i=0;i<N;i+=1)
			px = x0 + i*dx
			xyDistortion[i][j] = peakcorrection2(xymap,px,py)		// returns cmplx(dx,dy)
		endfor
	endfor
	totalDistortion = cabs(xyDistortion[p][q])

	if (strlen(WinList("GraphDistortionMap", "","WIN:1"))>1)
		DoWindow/F GraphDistortionMap
		SetDrawLayer /K UserFront	
	else
		Display/K=1/W=(168,55,798,615)
		DoWindow/C GraphDistortionMap
		AppendImage totalDistortion
		ModifyImage totalDistortion ctab= {0,7.1,Terrain,1}
		ModifyGraph margin(right)=72,width={Aspect,1},mirror=2
		ModifyGraph lblMargin(left)=7,lblMargin(bottom)=5,axOffset(left)=-2.57143,lblLatPos(left)=-44,lblLatPos(bottom)=47
		SetAxis/A/R left
		ColorScale/N=text0/F=0/S=3/A=RC/X=-13.63/Y=0.40 image=totalDistortion
	endif
	xyDistortion *= 20							// scale up to make lines longer
	for (j=0;j<N;j+=1)
		py = y0 + j*dy
		for (i=0;i<N;i+=1)
			px = x0 + i*dx
			SetDrawLayer UserFront
			SetDrawEnv xcoord=bottom,ycoord=left,arrow=1, arrowlen=7.0,arrowfat=0.3	// draw the arrow
			DrawLine px,py,px+real(xyDistortion[i][j]),py+imag(xyDistortion[i][j])
		endfor
	endfor
	KillWaves/Z xyDistortion
End

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

//	Make/O root:Packages:micro:X1correctionWave={0.31,0.32,0.2,0.09,-0.1,-0.3,-0.44,-0.47,-0.42,-0.37,-0.32,-0.19,-0.12,-0.07,0.01,0.17,0.35,0.52,0.67,0.73,0.69,0.53,0.32,0.13,-0.01,-0.13,-0.18,-0.18,-0.18,-0.16,-0.26,-0.25,-0.29,-0.38,-0.41,-0.38,-0.33,-0.18,0.01,0.23,0.42}
//	Make/O root:Packages:micro:Y1correctionWave={0.37,0.38,0.35,0.26,0.1,-0.06,-0.22,-0.31,-0.36,-0.35,-0.62,-0.51,-0.39,-0.28,-0.19,-0.03,0.11,0.25,0.34,0.38,0.38,0.34,0.28,0.17,0.04,-0.09,-0.22,-0.33,-0.4,-0.41,-0.39,-0.36,-0.28,-0.19,-0.08,0.06,0.2,0.3,0.4,0.54,0.66}
//	Make/O root:Packages:micro:Z1correctionWave={-0.09,-0.17,-0.21,-0.24,-0.25,-0.22,-0.14,-0.01,0.16,0.28,0.34,0.38,0.37,0.33,0.3,0.27,0.23,0.1,-0.05,-0.22,-0.18,-0.24,-0.28,-0.29,-0.28,-0.25,-0.18,-0.09,0,0.1,0.28,0.34,0.33,0.32,0.32,0.32,0.31,0.26,0.19,0.09,-0.02}
//	SetScale/P x 0,1,"µm", root:Packages:micro:X1correctionWave,root:Packages:micro:Y1correctionWave,root:Packages:micro:Z1correctionWave
//	SetScale d 0,0,"µm", root:Packages:micro:X1correctionWave,root:Packages:micro:Y1correctionWave,root:Packages:micro:Z1correctionWave
//
Function X1correctedKeyence(X1)		// takes PM500 X1 and returns the "real" X1
	Variable X1					// input the PM500 X1 reading
	Wave deltaW=root:Packages:micro:X1correctionWave
	Variable zWidth = DimDelta(deltaW,0)*(numpnts(deltaW)-1)	// width of the correction wave (µm)
	Variable dX1
	dX1 = mod(X1,zWidth)
	dX1 = dX1 < 0 ? dX1+zWidth : dX1
	return deltaW(dX1)+X1					// return corrected X1, this should be the true X1 value
End
//
Function Y1correctedKeyence(Y1)		// takes PM500 Y1 and returns the "real" Y1
	Variable Y1					// input the PM500 Y1 reading
	Wave deltaW=root:Packages:micro:Y1correctionWave
	Variable zWidth = DimDelta(deltaW,0)*(numpnts(deltaW)-1)	// width of the correction wave (µm)
	Variable dY1
	dY1 = mod(Y1,zWidth)
	dY1 = dY1 < 0 ? dY1+zWidth : dY1
	return deltaW(dY1)+Y1					// return corrected Y1, this should be the true Y1 value
End
//
Function Z1correctedKeyence(Z1)		// takes PM500 Z1 and returns the "real" Z1
	Variable Z1					// input the PM500 Z1 reading
	Wave deltaW=root:Packages:micro:Z1correctionWave
	Variable zWidth = DimDelta(deltaW,0)*(numpnts(deltaW)-1)	// width of the correction wave (µm)
	Variable dZ1
	dZ1 = mod(Z1,zWidth)
	dZ1 = dZ1 < 0 ? dZ1+zWidth : dZ1
	return deltaW(dZ1)+Z1					// return corrected Z1, this should be the true Z1 value
End

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

Static Function/S MakeUnique3Vector(preset)		// make a new 3-vector, and optionally set it to preset[]
	Wave preset							// optional wave, if none call with $""
	String wName = UniqueName("vec3_",1,0)
	Make/N=3/D $wName
	Wave vec=$wName
	if (WaveExists(preset))
		vec = preset
	endif
	return GetWavesDataFolder(vec,2)		// return full path for assingment to a wave
End
//
// Used to make temporary internal waves, a typical use would be:    Wave xxx=$MakeUnique3Vector($"")
// You will probably want to KillWaves/Z xxx, at the end of the routine.
Static Function/S MakeUnique3x3Mat(preset)		// make a new 3x3 matrix, and optionally set it to preset[]
	Wave preset							// optional wave, if none call with $""
	String wName = UniqueName("mat3x3_",1,0)
	Make/N=(3,3)/D $wName
	Wave mat=$wName
	if (WaveExists(preset))
		mat = preset
	endif
	return GetWavesDataFolder(mat,2)		// return full path for assingment to a wave
End

Function angleVec2Vec(a,b)			// return the angle between two vectors (degree)
	Wave a,b
	Variable dot = MatrixDot(a,b) / (norm(a)*norm(b))
	dot = limit(dot,-1,1)			// ensure that the acos will exist
	return acos(dot)*180/PI
End


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


Function Init_microGeo()
	if (NumVarOrDefault("root:Packages:geometry:geoInited",0))
		return 0
	endif
	Variable/G root:Packages:geometry:geoInited=0		// flag that initialization was called successsfully
	NVAR geoInited=root:Packages:geometry:geoInited
	NewDataFolder/O root:Packages						// ensure Packages exists
	NewDataFolder/O root:Packages:geometry			// ensure geometry exists
	NewDataFolder/O root:Packages:micro				// ensure micro exists
	Make/N=3/O/D root:Packages:geometry:pixel2q_ki, root:Packages:geometry:pixel2q_kout
	Make/N=3/O/D root:Packages:geometry:q2pixel_qhat,root:Packages:geometry:q2pixel_ki
	Make/N=3/O/D root:Packages:geometry:q2pixel_kout,root:Packages:geometry:q2pixel_vecB
	Make/N=(3,3)/O/D root:Packages:geometry:q2pixel_mat
	Make/N=3/O/D root:Packages:geometry:pixel2depth_ccd
	Make/N=3/O/D root:Packages:geometry:GeoUpdate_ki
	Make/N=(3,3)/O/D root:Packages:geometry:GeoUpdate_rde

	Make/O root:Packages:micro:X1correctionWave={0.31,0.32,0.2,0.09,-0.1,-0.3,-0.44,-0.47,-0.42,-0.37,-0.32,-0.19,-0.12,-0.07,0.01,0.17,0.35,0.52,0.67,0.73,0.69,0.53,0.32,0.13,-0.01,-0.13,-0.18,-0.18,-0.18,-0.16,-0.26,-0.25,-0.29,-0.38,-0.41,-0.38,-0.33,-0.18,0.01,0.23,0.42}
	Make/O root:Packages:micro:Y1correctionWave={0.37,0.38,0.35,0.26,0.1,-0.06,-0.22,-0.31,-0.36,-0.35,-0.62,-0.51,-0.39,-0.28,-0.19,-0.03,0.11,0.25,0.34,0.38,0.38,0.34,0.28,0.17,0.04,-0.09,-0.22,-0.33,-0.4,-0.41,-0.39,-0.36,-0.28,-0.19,-0.08,0.06,0.2,0.3,0.4,0.54,0.66}
	Make/O root:Packages:micro:Z1correctionWave={-0.09,-0.17,-0.21,-0.24,-0.25,-0.22,-0.14,-0.01,0.16,0.28,0.34,0.38,0.37,0.33,0.3,0.27,0.23,0.1,-0.05,-0.22,-0.18,-0.24,-0.28,-0.29,-0.28,-0.25,-0.18,-0.09,0,0.1,0.28,0.34,0.33,0.32,0.32,0.32,0.31,0.26,0.19,0.09,-0.02}
	SetScale/P x 0,1,"µm", root:Packages:micro:X1correctionWave,root:Packages:micro:Y1correctionWave,root:Packages:micro:Z1correctionWave
	SetScale d 0,0,"µm", root:Packages:micro:X1correctionWave,root:Packages:micro:Y1correctionWave,root:Packages:micro:Z1correctionWave

	if (exists("root:Packages:geometry:useDistortion")!=2)	// create global, it does not exist
		Variable/G root:Packages:geometry:useDistortion=USE_DISTORTION_DEFAULT
	endif
	NVAR useDistortion=root:Packages:geometry:useDistortion
	geoInited = 1
	if (useDistortion && exists("root:Packages:geometry:xymap")!=1)
		geoInited = LoadStandardDistortion() ? 0 : geoInited
	endif
	return 0
End