Kspace_JD.ipf

#pragma TextEncoding = "MacRoman"
//file: Kspace		created:  4/25/99  J. Denlinger
// Specialized conversion of BE-polar and hv-Polar intensity maps

#pragma rtGlobals=1		// Use modern global access method.
#pragma version=20.0425
// 4/25/20  shortcut to create K axis labels including unicode Overbar
// v2.53  10/16/15  subtract pol0, tilt0 from xrng, yrng in  PolTiltE2KxyE
// v2.52  2/4/15 fix poloff & Vinner globals in tmp folder
// v2.51  5/25/14 jdd PolarE2KxE:  change default kx step to be 0.01; fix poloff global typo
// v2.5 12/5/13 jdd  add off-normal constKpar parameter to PolPhot2Kxz()
// v2.4  11/19/12 jdd add geometry configuration options to PolarTilt2Kxy() & PolTiltE2KxyE()
// v2.3  5/26/12 jdd add kz2hvNE() routine
// v2.23 10/24/09 jdd swap X & Y in PolElev conversions to match BL7 (BJ Kim)
// v2.22 9/7/09 jdd add duplicate of CT rather than using original
// v2.21  3/5/09  jdd fix bugs in PolElevE2KxyE; add Pol-E-Elev input option 
//  v2.2 , 2/19/09 jdd  Add Kxyz-to-kxyhv for theory conversion to constant hv with intensity = band energy
//  v2.1, 12/13/08  jdd  Add photon momentum (approx) to PhotonE2KzE  only (so far)
//  v2.0, 11/29/08   jdd  Rename 2D comversion routines; Fix up Feng Wang's 3D routines 
//     goal have Proc interface; Function computation with KeyValues
// 2/12/08 jd change output plot fixed aspect ratio to be width follows height  instead of vice versa
// 08/16/04 fw  found can not specify the range in the panel, reason not figured out yet.
// 6/27/02 jdd      corrected BEpolar2Kx to allow KE data with hv=WF=0	


#include "Image_util"		// uses ColorTableStr; replaces Contour/ReGridToContour functionality
//#include "Sph2K"			// added required functions as in this procedure, can use static mode if call from function
//#include "Tilt2K"			// added required functions as static
#include "List_util"
#include "Kspace_calc"

//Function/T ImgRange( dir )
//Proc 	BEPolar2Kpar( img, outn, hv, WF, xrng, yrng, kopt, dopt ) 				: GraphMarquee
//Proc 	hvPolar2Kxz( img, outn, BE, WF, xrng, yrng, kxopt, kzopt, Vo, dopt ) 	: GraphMarquee

Menu "kspace"
	"-"
	"Convert Image on top graph..."
	"   Polar-E -> KxE", PolarE2KxE()
	"   Photon-E -> KzE", PhotonE2KzE()
	"   Polar-Tilt -> Kxy", PolTilt2kxky()
	"   Polar-hv -> Kxz", PolPhot2Kxz()
	"-"
	"Vol: Polar-Tilt-E -> KxyE", PolTiltE_KxyE()		
		help = {"convert Polar-Tilt-Energy 3d volume to Kx-Ky-Energy"}
	"Vol: Polar-hv-E -> KxzE", PolPhotE_KxzE()
		help = {"convert Polar-Photon-Energy 3d volume to Kx-Kz-Energy"}
	"Vol: Kxyz -> Kxyhv", Kxyz_Kxyhv()
		help = {"convert Kxyz (band energy) volume to Kxy-constant hv"}
	"-"
	"Create K axis label", CreateKaxisLabel()
	"Plot Kz points on NE hv", kz2hvNE()
End

Proc CreateKaxisLabel( kptList, kvalList,overbar, fontsz, sfont, basen, mode)
//-------------------
	string kptList =StrVarOrDefault( "root:tmp:KLBL_kptList", "M,G,Z")
	string kvalList =StrVarOrDefault( "root:tmp:KLBL_kvalList", "-1,0,1")
	variable overbar =NumVarOrDefault("root:tmp:KLBL_overbar", 1)
	variable fontsz =NumVarOrDefault("root:tmp:KLBL_fontsz", 4)
	string sfont =StrVarOrDefault("root:tmp:KLBL_sfont", "Arial")
	string basen =StrVarOrDefault( "root:tmp:KLBL_basen", "MGZ")
	variable mode =NumVarOrDefault("root:tmp:KLBL_mode", 2)
	prompt kptList, "K-point list, e.g M,G,X,Z"
	prompt kvalList, "K-value list, e.g -1,0,1,2"
	prompt overbar, "Add over-Bar", popup, "No;Yes"
	prompt fontsz, "Font size", popup, "None;12;14;16;18;20;24"
	Prompt sfont,"Select font", popup, "Times;Arial"
	prompt basen, "base name for _lbl, _k"
	prompt mode, "Output option", popup, "None;Table;Top graph X-axis;Table & Axis"
//	prompt mode, "Encode method", popup, "Unicode;old"

	string curr=GetDataFolder(1)
	NewDataFolder/O/S  root:tmp	
		string/G KLBL_basen=basen, KLBL_kptList=kptList, KLBL_kvalList=kvalList
		string/G KLBL_sfont=sfont
		variable/G KLBL_overbar=overbar, KLBL_fontsz=fontsz, KLBL_mode=mode
	SetDataFolder curr

	string TEXTwn=basen+"_lbl",VALwn=basen+"_k" 
// 
	variable nlbl = ItemsInList( kptList,",")
	
	make/o/n=(nlbl)/T $TEXTwn
	make/o/n=(nlbl)   $VALwn
	
	string textstr, kchar
	variable ii=0
	DO
		if (fontsz>1)
			textstr="\Z"+StringFromList( fontsz-2,"12;14;16;18;20;24")
		endif
		textstr+="\F'"+sfont+"'"

		kchar=StringFromList( ii,kptList,",")
		if (cmpstr(kchar,"G")==0)	// Gamma
			kchar=num2char(915)		// or num2char(0x383)
		endif
		if (overbar==2)    // 0x0305
			if (cmpstr(sfont,"Times")==0)
		 		kchar+="\x-36\y+14\u0305"
			else	//Arial
				kchar+="\x-26\y+06\u0305"
		 	endif
		endif
		textstr+=kchar
		$TEXTwn[ii] = textstr
		$VALwn[ii] = NumFromList( ii,kvalList, ",")
	
		ii+=1
	WHILE (ii<nlbl)

	string cmd
// Set axis on top plot
	if ((mode==3)+(mode==4))
		cmd = "ModifyGraph userticks(bottom)={"+VALwn+","+TEXTwn+"}"
		print cmd
		execute cmd
	endif
	
// Create table for fine editing 
	if ((mode==2)+(mode==4))
		DoWindow/F $(basen+"_labelTable")
		if (V_flag==0)
			edit $TEXTwn, $VALwn
		endif
		DoWindow/C $(basen+"_labelTable")
	endif
End

Proc kz2hvNE( c, sindex, kzlabel, Vo, WF)
//--------------------------
	variable c=NumVarOrDefault("root:NEhvKZ:clattice", 3.456)
	string sindex=StrVarOrDefault("root:NEhvKZ:skzindex", "4,8")
	string kzlabel=StrVarOrDefault("root:NEhvKZ:skzlabel", "G,Z")
	variable Vo=NumVarOrDefault("root:NEhvKZ:InnerPot", 10)
	variable WF=4.5
	prompt c, "vertical lattice constant, c="
	prompt sindex, "kz index range (start, stop)"
	prompt kzlabel, "kz labels, e.g. G, Z"
	prompt Vo, "Inner potential, Vo="
	prompt WF, "Work Function, WF="
	
// version 1 - make all variables and arrays in root:
	NewDataFolder/O/S  root:NEhvKZ	
		variable/G clattice=c, InnerPot=Vo, WorkFct=WF
		string/G skzindex=sindex, skzlabel=kzlabel
		
		variable index0=NumFromList(0, sindex, ","), index1=NumFromList(1, sindex, ",")
		variable np=abs(index1-index0) +1
	SetDataFolder root:

	make/o/n=(np) kzIndex, kzNE, hvNE, kzSYM
	make/o/n=(np)/T kzLBL
	kzindex=index0+p*sign(index1-index0)
	// assume binary repeat
	kzSym=19*(mod(kzIndex,2)==0) + 8*(mod(kzIndex,2)==1)
	kzLBL=StringFromList(mod(kzIndex,2), root:NEhvKZ:skzlabel, ",")
	kzNE := (pi/root:NEhvKZ:clattice)*kzindex
	hvNE := (kzNE/0.5123)^2 - root:NEhvKZ:InnerPot + root:NEhvKZ:WorkFct
	Setscale/P x 0,1E-6, hvNE
	
	CheckDisplayed/W=$WinName(0,1) hvNE
	if (V_flag==0)
		AppendToGraph hvNE
		ModifyGraph mode(hvNE)=3,zmrkNum(hvNE)={kzSYM}
		ModifyGraph rgb(hvNE)=(0,65535,65535), zmrkSize(hvNE)=5
	endif
	

	// make table	
	edit kzindex, kzNE,hvNE, kzLBL,kzSYM
	SetVariable setc title="c",size={70,15},value=:NEhvKZ:clattice,  limits={0.1,inf,0.01}
	SetVariable setV0 title="Vo",size={70,15},value=:NEhvKZ:InnerPot,  limits={-inf,inf,0.5}
	
	//	ControlInfo kwControlBar
	
	
	SetDataFolder root:

End


Proc PolarE2KxE( img, outn, hv, WF, xrng, yrng, kopt, thoff, dopt ) : GraphMarquee
// Formerly BEPolar2Kpar()
//---------------
	string img, outn=StrVarOrDefault( "root:tmp:outkx", "kx")
	variable hv=NumVarOrDefault( "root:tmp:hn", 100), WF=NumVarOrDefault( "root:tmp:Wfct", 4.5)
	variable dopt=1, kopt=NumVarOrDefault( "root:tmp:numkx", 0.01)
	variable thoff=NumVarOrDefault( "root:tmp:poloff", 0)
	string xrng=ImgRange(0,""), yrng=ImgRange(1,"")
	prompt img, "BE vs Polar Image to convert", popup, WaveList("!*ct", ";", "WIN:,DIMS:2")
	prompt outn, "Output image name or extension (kx)"
	prompt hv "Photon energy (eV) [0=KE data]"
	prompt WF "Work Function (eV) [0=KE data]"
	prompt xrng, "Polar range (min,max,inc)"
	prompt yrng, "BE range (min,max,inc)"
	prompt kopt, "kx inc (>0) or Rel. # of Kx points (<0)"
	prompt dopt, "Display new plot", popup, "Yes;No"
	prompt thoff, "Polar offset (deg)"
	NewDataFolder/O root:tmp
	string/G root:tmp:outkx=outn
	variable/G root:tmp:hn=hv, root:tmp:Wfct=WF, root:tmp:numkx=kopt, root:tmp:poloff=thoff
	
	PauseUpdate;Silent 1
	string winnam=WinName(0,1)
	
//Check for marquee; subset selection of image
// obtained from ImgRange() function
	variable x1, x2, xinc, nx, y1, y2, yinc, ny
	if (strlen(xrng)>0)
		x1=ValFromList(xrng, 0,",")
		x2=ValFromList(xrng, 1,",")
		xinc=ValFromList(xrng, 2,",")
	else					//use full range of image
		xinc=DimDelta($img, 0)
		x1=DimOffset($img, 0)	
		x2=x1+(DimSize($img, 0)-1)*xinc
	endif
	nx=round( (x2-x1)/xinc + 1)
		
	if (strlen(yrng)>0)
		y1=ValFromList(yrng, 0,",")
		y2=ValFromList(yrng, 1,",")
//		yinc=ValFromList(yrng, 2,",")
		yinc=DimDelta($img, 1)
	else					//use full range of image
	//print img
		yinc=DimDelta($img, 1)
		y1=DimOffset($img, 1)	
		y2=y1+(DimSize($img, 1)-1)*yinc
	endif
	ny=round( (y2-y1)/yinc + 1)
	
	Duplicate/O/R=(x1,x2)(y1,y2) $img, root:tmp:im
	//print x1, x2, xinc, nx
	print "y:", y1, y2, yinc, ny
	
// Calculate kx range
	variable kx1, kx2, kinc, nk, BEmax=max(y1, y2)
	BEmax=SelectNumber( hv>0, BEmax, min( BEmax, 0) )  // use EF for highest energy
	////print BEmax, hv+BEmax-WF
	kx1=sph2kx_(hv+BEmax-WF, x1-thoff,  0)
	kx1=1E-2*round(kx1*1E2)
	kx2=sph2kx_(hv+BEmax-WF, x2-thoff,  0)
	kx2=1E-2*round(kx2*1E2)
	if (kopt<0)
		nk=round(nx*abs(kopt))
		kinc=(kx2-kx1)/(nk-1)
	else
		kinc=kopt
		nk=abs(round((kx2-kx1)/kinc + 1))
	endif
	print "kx:", kx1, kx2, kinc, nk

// Create new image & theta array	; set BE-Kpar regular grid scaling
	if (cmpstr(outn[0], "k")==0)
		outn=img+outn
	endif	
	Duplicate/o root:tmp:im $outn			//, root:tmp:theta
	Redimension/N=(nk, ny) $outn			//, root:tmp:theta
	Setscale/P x kx1, kinc, "" $outn			//,  root:tmp:theta
	Setscale/P y y1, yinc, "" $outn			//,  root:tmp:theta
	//print "Created: ", outn, "(", nk, "x", ny, ") Æk=", kinc

// Calculate irregular grid theta values	
	//root:tmp:theta=Kp2Polar_( hv+y-WF,  x)
	//$outn=interp2D( $img, root:tmp:theta+thoff , y)
	//	all in one command without intermediate array
	$outn=interp2D( $img, Kp2Polar_( hv+y-WF,  x)+thoff , y)
	
	DoWindow/F $("Kx_"+winnam)
	if ((dopt==1)*(V_flag==0))
		// get color table info from top graph(?)
		string CTstr=ColorTableStr("", img)		//may not work if not proper top graph
		//print outn
		display; appendimage $outn
		Label bottom "Kx (1/)"
		Label left "BE (eV)"
		if (strlen(CTstr)>0) 
			// duplicate CT instead of using old one
			Duplicate/o $(img+"_CT") $(outn+"_CT")
			execute "ModifyImage "+outn+" cindex="+outn+"_CT"
//			execute "ModifyImage "+outn+" "+CTstr
		else
			string ctable=img+"_CT"
			if (exists(ctable)==0)
				//ctable="YellowHot"
				ModifyImage $outn ctab= {*,*,YellowHot,0}
			else
				ModifyImage $outn cindex=$ctable
			endif
		endif
		DoWindow/C $("Kx_"+winnam)
	endif
End

Proc PhotonE2KzE( img, outn, V0, WF, xrng, yrng, kopt, dopt, khvang ) : GraphMarquee
// Formerly BEhv2Kz()
//---------------
// Normal Emission photon-dependence
	string img, outn=StrVarOrDefault( "root:tmp:outkz", "kz")
	variable V0=NumVarOrDefault( "root:tmp:Vinner", 10), WF=NumVarOrDefault( "root:tmp:Wfct", 4.5)
	variable dopt=1, kopt=NumVarOrDefault( "root:tmp:numkx", 0.02), khvang=-1
	string xrng=ImgRange(0,""), yrng=ImgRange(1,"")
	prompt img, "BE vs hv Image to convert", popup, WaveList("!*ct", ";", "WIN:,DIMS:2")
	prompt outn, "Output image name or extension (kz)"
	prompt V0 "Inner Potential (eV)"
	prompt WF "Work Function (eV)"
	prompt xrng, "Photon range (min,max)"
	prompt yrng, "BE range (min,max)"
	prompt kopt, "kz inc (>0) or Rel. # of Kz points (<0)"
	prompt dopt, "Display new plot", popup, "Yes;No"
	prompt khvang, "hv incid. angle (<0, no calc)"
	
	NewDataFolder/O root:tmp
	string/G root:tmp:outkz=outn
	variable/G root:tmp:Vinner=V0, root:tmp:Wfct=WF, root:tmp:numkx=kopt
	
	PauseUpdate;Silent 1
	string winnam=WinName(0,1)
	
	
	// photon momentum option
	variable khv_perp=0, POL=0, ELV=0, d2r=pi/180
	if (khvang>0)
		khv_perp = 5.07E-4*cos(d2r*(khvang+POL))*cos(d2r*ELV)		//*photon energy
	endif
	
//Check for marquee; subset selection of image
// obtained from ImgRange() function
	variable x1, x2, xinc, nx, y1, y2, yinc, ny
	if (strlen(xrng)>0)
		x1=ValFromList(xrng, 0,",")
		x2=ValFromList(xrng, 1,",")
		xinc=ValFromList(xrng, 2,",")
	else					//use full range of image
		xinc=DimDelta($img, 0)
		x1=DimOffset($img, 0)	
		x2=x1+(DimSize($img, 0)-1)*xinc
	endif
	nx=round( (x2-x1)/xinc + 1)
		
	if (strlen(yrng)>0)
		y1=ValFromList(yrng, 0,",")
		y2=ValFromList(yrng, 1,",")
		//yinc=ValFromList(yrng, 2,",")
		yinc=DimDelta($img, 1)
	else					//use full range of image
		yinc=DimDelta($img, 1)
		y1=DimOffset($img, 1)	
		y2=y1+(DimSize($img, 1)-1)*yinc
	endif
	ny=round( (y2-y1)/yinc + 1)
	
	Duplicate/O/R=(x1,x2)(y1,y2) $img, root:tmp:im
	//print x1, x2, xinc, nx
	//print y1, y2, yinc, ny

// Calculate kz range
	variable kz1, kz2, kinc, nk 		//, BEmax=min( max(y1, y2), 0)
	variable BEmax= min(y1,y2)   //maximum negative BE
	variable BEmin = max(y1,y2) //positive BE above EF
//	print BEmin, BEmax
	kz1=sph2kz_(min(x1,x2)+BEmax-WF, 0,  0, V0)  +khv_perp*min(x1,x2)
	kz1=1E-2*round(kz1*1E2)
	kz2=sph2kz_(max(x1,x2)+BEmin-WF, 0,  0, V0)  +khv_perp*max(x1,x2)
	kz2=1E-2*round(kz2*1E2)
	if (kopt<0)
		nk=round(nx*abs(kopt))
		kinc=(kz2-kz1)/(nk-1)
	else
		kinc=kopt
		nk=round((kz2-kz1)/kinc + 1)
	endif
	print "kz:", kz1, kz2, kinc, nk

// Create new image & theta array	; set BE-Kpar regular grid scaling
	if (cmpstr(outn[0], "k")==0)
		outn=img+outn
	endif	
	Duplicate/o root:tmp:im $outn			//, root:tmp:theta
	Redimension/N=(abs(nk), abs(ny)) $outn			//, root:tmp:theta
	Setscale/P x kz1, kinc, "" $outn			//,  root:tmp:theta
	//print "Created: ", outn, "(", nk, "x", ny, ") Æk=", kinc

// Calculate irregular grid theta values	
	//root:tmp:theta=Kp2Polar_( hv+y-WF,  x)
	//$outn=interp2D( $img, root:tmp:theta+thoff , y)
	//	all in one command without intermediate array
		

	
	$outn=interp2D( $img, K2KE_(0, x-khv_perp*(K2KE_(0,x)-V0+WF) )-V0 +WF -y , y)			//assumes scaled
//	$outn=interp2D( $img, K2KE_(0)-V0 +WF -y , y)			//assumes scaled
	// Kp2Polar_  --> Kz2hv:  = 3.81*k^2 -V0 +WF-BE
	
	DoWindow/F $("Kz_"+winnam)
	if ((dopt==1)*(V_flag==0))
		// get color table info from top graph(?)
		string CTstr=ColorTableStr("", img)		//may not work if not proper top graph
		display; appendimage $outn
		Label bottom "Kz (1/)"
		Label left "BE (eV)"
		if (strlen(CTstr)>0) 
			execute "ModifyImage "+outn+" "+CTstr
		else
			string ctable=img+"_CT"
			if (exists(ctable)==0)
				//ctable="YellowHot"
				ModifyImage $outn ctab= {*,*,YellowHot,0}
			else
				ModifyImage $outn cindex=$ctable
			endif
		endif
		DoWindow/C $("Kz_"+winnam)
	endif
End



Function Kcvt_Ang2Kxy( img, inopt, outopt, [mode] )		// mode=#   skip various sections
//================
// Input options:
//		A-ranges:
//			(default) full range of scaled image
//			/X=xmin,xmax;  /Y=ymin,ymax
//		/KE=val
//		/BE=val/WF=val/hv=val   (KE=hv-BE-WF)
//   	/XOFF=val, /YOFF=val   or? /XYOFF=xoff, yoff
//      	/AZOFF=val   (azimuth offset rotation)
//		/GEO = 1,2,3    (1=Det Polar Elev, 2= Sample Poler-tilt+Horiz Slit, 3=Sample Polar-Tilt+Vert Slit
// Output options:
//       /D=kout    (default = *_kxy)
//       K-ranges:
//             (default) minimal span of A->K conversion
//             /KX=kxmin,kxmax,    /KY=kymin, kymax
//             /KINC=kinc  or /KINC=kxinc,kyinc
// optional parameter:
//     mode=#    1=default=execute all steps
//                       2, 3,4 = skip various setup steps
// Usage:
//     
//
	Wave img
	string inopt, outopt
	variable mode
	string imgn=NameOfWave(img)
	
		
	if (ParamIsDefault(mode))
		mode=1
	endif
	print "mode=", mode
	
	NewDataFolder/O/S root:kcvt
//		string/G outkxy=koutn, kincr=kinc, ang_off=angoff
//		variable/G root:kcvt:ken=KE, azim_off=azimoff
		//variable/G root:kcvt:numkx=kxopt, root:kcvt:numky=kyopt
		//variable/G root:kcvt:poloff=pol_off, root:kcvt:elvoff=elv_off
	SetDataFolder root:
	

	// Determine input ranges
	//  Kcvt_InputRange(img, inopt)  //--> variables in folder
	//  
	string xrng=KeyStr("X", inopt), yrng=KeyStr("Y",inopt)
	variable  x1, x2, xinc, nx, xmid
	if (strlen(xrng)>0)
		x1=ValFromList(xrng, 0,",")
		x2=ValFromList(xrng, 1,",")
		xinc=ValFromList(xrng, 2,",")
	else					//default use full range of image
		xinc=DimDelta(img, 0)
		x1=DimOffset(img, 0)	
		x2=x1+(DimSize(img, 0)-1)*xinc
	endif
	nx=round( (x2-x1)/xinc + 1)
	xmid = (x1+x2)/2
	
	variable  y1, y2, yinc, ny, ymid
	if (strlen(yrng)>0)
		y1=ValFromList(yrng, 0,",")
		y2=ValFromList(yrng, 1,",")
		yinc=ValFromList(yrng, 2,",")
	else					//default use full range of image
		yinc=DimDelta(img, 1)
		y1=DimOffset(img, 1)	
		y2=y1+(DimSize(img, 1)-1)*yinc
	endif
	ny=round( (y2-y1)/yinc + 1)
	ymid = (y1+y2)/2
	
	// use x1,x2,xmid (& y) to populate range arrays for plotting purposes,
	//     k-range conversion & k-range evaluation
	make/o/n=9 xrngw, yrngw
	xrngw={x2, x2, xmid, x1, x1, x1, xmid, x2, x2} 
	yrngw={ymid, y2, y2, y2, ymid, y1, y1, y1, ymid}
	
	//Determine constant kinetic energy:
	variable KE=95.5, hv=100, BE=0, WF=4.5
	if (KeySet("KE", inopt))
		KE=KeyVal("KE",inopt)
	else
		hv=KeyVal("hv",inopt)
		// possibly look for photon energy from wave note
		BE=KeyVal("BE",inopt)		// positive values below EF
		WF=KeyVal("WF",inopt)
		KE=hv-BE-WF
	endif
	
	//Determine offset parameters
	variable xoff=0, yoff=0, azoff=0
	if (KeySet("AZOFF",inopt))
		azoff=-KeyVal("XYOFF", inopt)
	endif
	if (KeySet("XYOFF", inopt))
		xoff=str2num(StringFromList(0, KeyStr("XYOFF", ",")) )
		xoff=SelectNumber(numtype(xoff)==2, xoff, 0)		// use default if NAN
		yoff=str2num(StringFromList(1, KeyStr("XYOFF", ",")) )
		yoff=SelectNumber(numtype(yoff)==2, yoff, 0)    //NaN for single value list
	endif
	
	// Geometry formula selection
	variable geometry=3
	if (KeySet("GEO", inopt))
		geometry=KeyVal("GEO",inopt)
	endif
	
	//Convert xy-range (angles) to kxy
	make/o/n=9 kxrngw, kyrngw
	kxrngw=DualA2kx_( KE, xrngw, yrngw, geometry)
	kyrngw=DualA2ky_( KE, xrngw, yrngw, geometry)
	
	//Determine kmin, kmax values & round
	variable kxmin, kxmax, kymin, kymax
	WaveStats/Q kxrngw
	kxmin=1E-2*round(V_min*1E2)
	kxmax=1E-2*round(V_max*1E2)
	WaveStats/Q kyrngw
	kymin=1E-2*round(V_min*1E2)
	kymax=1E-2*round(V_max*1E2)

	// output k-increment, numpnts
	variable kxinc=0.02, kyinc=0.02
	if (KeySet("KINC", outopt))
		kxinc=str2num(StringFromList(0, KeyStr("KINC", ",")) )
		kxinc=SelectNumber(numtype(kxinc)==2, kxinc, 0.02)		// use default if NAN
		kyinc=str2num(StringFromList(1, KeyStr("KINC", ",")) )
		kyinc=SelectNumber(numtype(kyinc)==2, kyinc, kxinc)    //NaN for single value list
	endif
	
	variable nkx, nky	
	nkx=round((kxmax-kxmin)/kxinc + 1)
	nky=round((kymax-kymin)/kyinc + 1)
	//print "kx:", kx1, kx2, kxinc, nkx, "   ky:", ky1, ky2, kyinc, nky
	
	// create scaled output k-array
	string koutn=KeyStr("D", outopt)
	if ((strlen(koutn)==0)+stringmatch(imgn,koutn))
		koutn=imgn+"_kxy"
	endif	
	
//	Make/O/N=(nkx, nky)  $koutn, root:kcvt:polw, root:kcvt:tiltw
//	Setscale/I x kxmin, kxmax, "" $koutn,  root:kcvt:polw, root:kcvt:tiltw
//	Setscale/I y kymin, kymax, "" $koutn,  root:kcvt:polw, root:kcvt:tiltw	
	Make/O/N=(nkx, nky) 	kout, polw, tiltw
	WAVE kout=kout, polw=polw, tiltw=elevw	
	Setscale/I x kxmin, kxmax, "" kout, polw, tiltw
	Setscale/I y kymin, kymax, "" kout, polw, tiltw
	
	// Calculate irregular grid polar& elev values
	// worry about misorientation correction

	variable r2d=180/pi, ktot=0.5123*sqrt(KE)	
	//root:kcvt:polw=r2d*asin(x/sqrt(ktot^2-y^2) )
	//root:kcvt:tiltw=r2d*asin(y/ktot)
	polw=r2d*asin(xrot(x, y, azoff)/sqrt(ktot^2-yrot(x, y, azoff)^2) )
	tiltw=r2d*asin(yrot(x, y, azoff)/ktot)
	//kout=interp2D( img, root:kcvt:theta , y)
	//	all in one command without intermediate array
	kout=interp2D( img, polw+xoff , tiltw+yoff )
//	kout=interp2D( img, r2d*asin( xrot(x, y, azoff)/sqrt(ktot^2-yrot(x, y, azoff)^2) )+xoff , r2d*asin(yrot(x, y, azoff)/ktot)+yoff )
	//$outn=interp2D( $img, r2d*asin( x/sqrt(kout^2-y^2) )+pol0 , r2d*asin(y/kout)+elev0 )
	//$outn=interp2D( $img, Kp2Polar_( K2KE_(x, y)-Vo+BE,  x)+pol0 , K2KE_(x, y)-Vo-BE+WF )

	
End


Proc PolTilt2kxky( img, outn, xrng, yrng, kinc, angoff, KE, azimoff,  geometry,dopt ) : GraphMarquee
//---------------
	string img, outn=StrVarOrDefault( "root:kcvt:outkxy", "kxy")
	string kinc=StrVarOrDefault( "root:kcvt:kincr", ".02,.02")
		//variable kxopt=NumVarOrDefault( "root:kcvt:numkx", .1), kyopt=NumVarOrDefault( "root:kcvt:numky", .1)
	string angoff=StrVarOrDefault( "root:kcvt:ang_off", "0,0")
		//variable pol_off=NumVarOrDefault( "root:kcvt:poloff", 0)
		//variable elv_off=NumVarOrDefault( "root:kcvt:elvoff", 0)
	variable KE=NumVarOrDefault( "root:kcvt:ken", 0)
	variable azimoff=NumVarOrDefault( "root:kcvt:azim_off", 0)
	variable geometry=NumVarOrDefault( "root:kcvt:geometry", 3)
	variable dopt=1
	string xrng=ImgRange(0,""), yrng=ImgRange(1,"")
	prompt img, "Elev vs Polar Image to convert", popup, WaveList("!*ct", ";", "WIN:,DIMS:2")
	prompt outn, "Output image name or extension (kxy)"
	prompt xrng, "Polar range (min,max,inc)"
	prompt yrng, "Elev range (min,max,inc)"
	prompt kinc, "k-increments:  kx,ky"
	prompt angoff, "angle offsets:  pol0, elev0"
	prompt KE "Kinetic energy (eV) [hv-BE-WF]"
	prompt azimoff "Azimuth rotation (deg)"
	//prompt kxopt, "kx inc (>0) or Rel. # of Kx points (<0)"
	//prompt kyopt, "ky inc (>0) or Rel. # of Ky points (<0)"
	//prompt pol_off, "Polar offset (deg)"
	//prompt elv_off, "Elevation offset (deg)"
	prompt geometry "Geometry", popup, "Detector Pol-Elev;Sample Ang+Horiz Slit;Sample Ang+Vert Slit"
	prompt dopt, "Display new plot", popup, "Yes;No"
	
	NewDataFolder/O/S root:kcvt
		string/G outkxy=outn, kincr=kinc, ang_off=angoff
		variable/G ken=KE, azim_off=azimoff
		variable/G geometry=geometry
		//variable/G root:kcvt:numkx=kxopt, root:kcvt:numky=kyopt
		//variable/G root:kcvt:poloff=pol_off, root:kcvt:elvoff=elv_off
	SetDataFolder root:
	
	PauseUpdate;Silent 1
	string winnam=WinName(0,1)
	
//Check for marquee; subset selection of image
// obtained from ImgRange() function
	variable x1, x2, xinc, nx, y1, y2, yinc, ny
	if (strlen(xrng)>0)
		x1=ValFromList(xrng, 0,",")
		x2=ValFromList(xrng, 1,",")
		xinc=ValFromList(xrng, 2,",")
	else					//use full range of image
		xinc=DimDelta($img, 0)
		x1=DimOffset($img, 0)	
		x2=x1+(DimSize($img, 0)-1)*xinc
	endif
	nx=round( (x2-x1)/xinc + 1)
		
	if (strlen(yrng)>0)
		y1=ValFromList(yrng, 0,",")
		y2=ValFromList(yrng, 1,",")
		yinc=ValFromList(yrng, 2,",")
	else					//use full range of image
		yinc=DimDelta($img, 1)
		y1=DimOffset($img, 1)	
		y2=y1+(DimSize($img, 1)-1)*yinc
	endif
	ny=round( (y2-y1)/yinc + 1)
	
	Duplicate/O/R=(x1,x2)(y1,y2) $img, root:kcvt:im
	////print x1, x2, xinc, nx
	////print y1, y2, yinc, ny

//Extract angle offsets
	variable pol0, elev0
	pol0=ValFromList( angoff, 0, ",")
	elev0=ValFromList( angoff, 1, ",")
	////print pol0, elev0, azimoff

// Calculate k-grid range   -- add geometry dependence
	variable ang=-azimoff
	variable kx1, kx2, kx_1, kx_2, kxinc, nkx
	variable ky1, ky2, ky_1, ky_2, kyinc, nky
		//convert corners	
	kx1=dualA2kx_(KE, x1-pol0,  0, geometry);  
	kx2=dualA2kx_(KE, x2-pol0,  0, geometry)	   // max at ky=0
	ky1=dualA2ky_(KE, 0, y1-elev0, geometry); 
	ky2=dualA2ky_(KE, 0, y2-elev0, geometry)    //straight lines in kspace, any kx ok
		//sort min/max
	kx_1=min(kx1, kx2); kx_2=max(kx1, kx2)
	ky_1=min(ky1, ky2); ky_2=max(ky1, ky2)
		//rotate
	kx1= min( xrot(kx_1, ky_1, -ang),  xrot(kx_1, ky_2, -ang) )
	kx2=max( xrot(kx_2, ky_1, -ang),  xrot(kx_2, ky_2, -ang) )
	ky1= min( yrot(kx_1, ky_1, -ang),  yrot(kx_2, ky_1, -ang) )
	ky2=max( yrot(kx_1, ky_2, -ang),  yrot(kx_2, ky_2, -ang) )
		//round to simpler numbers (2 decimeal)
	kx1=1E-2*round(kx1*1E2)
	kx2=1E-2*round(kx2*1E2)
	ky1=1E-2*round(ky1*1E2)
	ky2=1E-2*round(ky2*1E2)
	
	kxinc=ValFromList( kinc, 0, ",")
	kyinc=ValFromList( kinc, 1, ",")
	
	nkx=round((kx2-kx1)/kxinc + 1)
	nky=round((ky2-ky1)/kyinc + 1)

	//print "kx:", kx1, kx2, kxinc, nkx, "   ky:", ky1, ky2, kyinc, nky

// Create new image & polar/elev arrays	; set Kx-Ky regular grid scaling
	if (cmpstr(outn[0], "k")==0)
		outn=img+outn
	endif	
	Make/O/N=(nkx, nky) $outn, root:kcvt:polw, root:kcvt:tiltw
	Setscale/I x kx1, kx2, "" $outn,  root:kcvt:polw, root:kcvt:tiltw
	Setscale/I y ky1, ky2, "" $outn,  root:kcvt:polw, root:kcvt:tiltw
	//print "Created: ", outn, "(", nkx, "x", nky, ") Ækxy=", kxinc, kyinc

// Calculate irregular grid polar& elev values
	// worry about misorientation correction

	variable r2d=180/pi, kout=0.5123*sqrt(KE)	
	IF (geometry==1)				// Single channel Detector Polar-Elevation +  fixed sample
	      if (azimoff==0)			      // no azimuth angle correction
			root:kcvt:polw=r2d*asin(x/sqrt(kout^2-y^2) )
			root:kcvt:tiltw=r2d*asin(y/kout)
		else
			root:kcvt:polw=r2d*asin(xrot(x, y, ang)/sqrt(kout^2-yrot(x, y, ang)^2) )
			root:kcvt:tiltw=r2d*asin(yrot(x, y, ang)/kout)
		endif
		$outn=interp2D( $img, root:kcvt:polw + pol0 , root:kcvt:tiltw + elev0)
		// all-in-one command w/o intermediate arrays		
//		$outn=interp2D( $img, r2d*asin( xrot(x, y, ang)/sqrt(kout^2-yrot(x, y, ang)^2) )+pol0 , r2d*asin(yrot(x, y, ang)/kout)+elev0 )	

	ELSE
	IF (geometry==2)			// Sample Polar-Tilt +  Detector Horizonatal Slit  (polw=detector angle)
		 if (azimoff==0)			      // no azimuth angle correction
			root:kcvt:polw=r2d*asin(x/kout)
			root:kcvt:tiltw=r2d*asin(y/sqrt(kout^2-x^2) )
		else						 // with azimuth angle correction to kx-ky
			root:kcvt:polw=r2d*asin(xrot(x, y, ang)/kout)
			root:kcvt:tiltw=r2d*asin(yrot(x, y, ang)/sqrt(kout^2-xrot(x, y, ang)^2) )
		endif
		$outn=interp2D( $img, root:kcvt:polw + pol0 , root:kcvt:tiltw + elev0)
		//	all in one command without intermediate array
//		$outn=interp2D( $img, r2d*asin(xrot(x, y, ang)/kout)+pol0, r2d*asin( yrot(x, y, ang)/sqrt(kout^2-xrot(x, y, ang)^2) )+elev0  )
	
	ELSE
	IF (geometry==3)			// Sample Polar-Tilt +  Detector Vertical Slit  (polw=detector angle)
		 if (azimoff==0)			      // no azimuth angle correction
			root:kcvt:tiltw=r2d*asin(y/kout*cos(elev0/r2d) + sqrt(kout^2-x^2-y^2)/kout*sin(elev0/r2d) )
			root:kcvt:polw=r2d*asin(x/kout/cos(root:kcvt:tiltw/r2d) )
		else						 // with azimuth angle correction to kx-ky
			root:kcvt:tiltw=r2d*asin(yrot(x, y, ang)/kout*cos(elev0/r2d) + sqrt(kout^2-xrot(x, y, ang)^2-yrot(x, y, ang)^2)/kout*sin(elev0/r2d) )
			root:kcvt:polw=r2d*asin(xrot(x, y, ang)/kout/cos(root:kcvt:tiltw/r2d) )
		endif	
//			root:kcvt:tiltw=r2d*asin(y/kout)
		$outn=interp2D( $img, root:kcvt:polw + pol0 , root:kcvt:tiltw  )			// elev0 offset included above
	
	ENDIF
	ENDIF
	ENDIF

//Azimuthally rotate image; ?? about what point; preservation of image scaling?
//	if (azimoff!=0)
//		ImageRotate/A=-azimoff/E=-1/O $outn
//		//$outn*=SelectNumber( $outn==-1, 1, Nan)
//	endif
	
	DoWindow/F $("Kxy_"+winnam)
	if ((dopt==1)*(V_flag==0))
		// get color table info from top graph(?)
		string CTstr=ColorTableStr("", img)		//may not work if not proper top graph

		
		
		display; appendimage $outn

//		ModifyGraph height={Plan,1,left,bottom}
		ModifyGraph width={Plan,1,bottom,left}
		ModifyGraph axThick=0.5,standoff=0
		Label bottom "Kx (1/)"
		Label left "Ky (1/)"
		string titlestr="\JC"+outn+": "+num2str(nkx)+" x "+num2str(nky)
		titlestr+="\rKE="+num2str(KE)+",  NE=("+num2str(pol0)+", "+num2str(elev0)+", "+num2str(azimoff)+")"
		Textbox/N=title/F=0/A=MT/E titlestr
		//worry about changes w/o display
		if (strlen(CTstr)>0) 
			// duplicate CT instead of using old one
			Duplicate/o $(img+"_CT") $(outn+"_CT")
			execute "ModifyImage "+outn+" cindex="+outn+"_CT"
//			execute "ModifyImage "+outn+" "+CTstr
		else
			string ctable=img+"_CT"
			if (exists(ctable)==0)
				//ctable="YellowHot"
				ModifyImage $outn ctab= {*,*,YellowHot,0}
			else
				ModifyImage $outn cindex=$ctable
			endif
		endif
		DoWindow/C $("Kxy_"+winnam)
	endif

End

Function xrot( x, y, ang )
	variable x, y, ang
	ang*=pi/180
	return x*cos(ang)+y*sin(ang)
end

Function yrot( x, y, ang )
	variable x, y, ang
	ang*=pi/180
	return -x*sin(ang)+y*cos(ang)
end

Proc PolPhot2Kxz( img, outn, BE_WF, constKpar, xrng, yrng, kinc, thoff, Vo,dopt ) : GraphMarquee
//   Formerly hvPolar2Kxz()
//-------------------
	string img, outn=StrVarOrDefault( "root:tmp:outkxz", "kxz")
	string BE_WF=StrVarOrDefault( "root:tmp:BE_WFct", "0,4.5")
//	variable BE=NumVarOrDefault( "root:tmp:ben", 0), WF=NumVarOrDefault( "root:tmp:Wfct", 4.5)
	variable constKpar=NumVarOrDefault( "root:tmp:const_Kpar", 0)
	string kinc=StrVarOrDefault( "root:tmp:kxzincr", ".02,.02")
	//variable kxopt=NumVarOrDefault( "root:tmp:numkx", .1), kzopt=NumVarOrDefault( "root:tmp:numkz", .1)
	variable Vo=NumVarOrDefault( "root:tmp:Vinner", 10)	//, dopt=1 negative for varying Vo
	variable thoff=NumVarOrDefault( "root:tmp:poloff", 0)
	string xrng=ImgRange(0,""), yrng=ImgRange(1,"")
	variable dopt=1
	prompt img, "hv vs Polar Image to convert", popup, WaveList("!*ct", ";", "WIN:,DIMS:2")
	prompt outn, "Output image name or extension (kxz)"
	prompt BE_WF "Binding energy, Work Function (eV)"
	prompt constKpar "Constant K_parallel"
	prompt xrng, "Polar range (min,max,inc)"
	prompt yrng, "hv range (min,max,inc)"
	prompt kinc, "k-increments:  kx,kz"
		//prompt kxopt, "kx inc (>0) or Rel. # of Kx points (<0)"
		//prompt kzopt, "kz inc (>0) or Rel. # of Kz points (<0)"
	prompt Vo, "Inner Potential, Vo,negative for varying Vo"   
	prompt thoff, "Polar offset (deg)"
	prompt dopt, "Display new plot", popup, "Yes;No"
	
	
	//variable Vo_base=NumVarOrDefault( "root:tmp:Vo_base", 0)	//, dopt=1      //FW 0507/05
	//variable Vo_slope=NumVarOrDefault( "root:tmp:Vo_slope", 0)	//, dopt=1       //FW 0507/05
	//variable hv_base=NumVarOrDefault( "root:tmp:hv_base", 0)	//, dopt=1       //FW 0507/05
	//prompt Vo_base, "or  Vo="									  //FW 0507/05
	//prompt Vo_slope, "+                     *(hv-"					  //FW 0507/05
	//prompt hv_base, ")"					  //FW 0507/05
	
	NewDataFolder/O root:tmp
	string/G root:tmp:outkxz=outn, root:tmp:kxzincr=kinc
	string/G root:tmp:BE_WFct=BE_WF
//	variable/G root:tmp:ben=BE, root:tmp:Wfct=WF
	variable/G root:tmp:const_Kpar=constKpar
	variable/G root:tmp:poloff=thoff, root:tmp:Vinner=Vo
	//variable/G root:tmp:numkx=kxopt, root:tmp:numkz=kzopt

// Extract BE and WF values from string
	variable BE, WF
	BE = NumFromList(0, BE_WF, ",")
	WF = NumFromList(1, BE_WF, ",")
	print BE, WF

	
	GetVo(Vo)
	variable Vo_base=root:tmp:Vo_base	//, dopt=1      //FW 0507/05
	variable Vo_slope=root:tmp:Vo_slope	//, dopt=1       //FW 0507/05
	variable hv_base=root:tmp:hv_base	//, dopt=1       //FW 0507/05
	
	PauseUpdate;Silent 1
	string winnam=WinName(0,1)
	
//Check for marquee; subset selection of image
// obtained from ImgRange() function
	variable x1, x2, xinc, nx, y1, y2, yinc, ny
	if (strlen(xrng)>0)
		x1=ValFromList(xrng, 0,",")
		x2=ValFromList(xrng, 1,",")
		xinc=ValFromList(xrng, 2,",")
	else					//use full range of image
		xinc=DimDelta($img, 0)
		x1=DimOffset($img, 0)	
		x2=x1+(DimSize($img, 0)-1)*xinc
	endif
	nx=round( (x2-x1)/xinc + 1)
		
	if (strlen(xrng)>0)
		y1=ValFromList(yrng, 0,",")
		y2=ValFromList(yrng, 1,",")
		yinc=ValFromList(yrng, 2,",")
	else					//use full range of image
		yinc=DimDelta($img, 1)
		y1=DimOffset($img, 1)	
		y2=y1+(DimSize($img, 1)-1)*yinc
	endif
	ny=round( (y2-y1)/yinc + 1)
	
	
	Duplicate/O/R=(x1,x2)(y1,y2) $img, root:tmp:im
	////print x1, x2, xinc, nx
	////print y1, y2, yinc, ny
	
// Calculate kx range
	variable hvmax=max(y1, y2), hvmin=min(y1,y2)
	variable Vomin=Vo_base+Vo_slope*(hvmin-hv_base)        //FW 0507/05
	variable Vomax=Vo_base+Vo_slope*(hvmax-hv_base)        //FW 0507/05
	
	variable kx1, kx2, kxinc, nkx, 
	kx1=min( sph2kx_(hvmax+BE-WF, x1-thoff,  0), sph2kx_(hvmin+BE-WF, x1-thoff,  0) )
	kx2=max( sph2kx_(hvmax+BE-WF, x2-thoff,  0), sph2kx_(hvmin+BE-WF, x2-thoff,  0) )
	kx1=1E-2*round(kx1*1E2)
	kx2=1E-2*round(kx2*1E2)	
	
// Calculate kz range
	variable kz1, kz2, kzinc, nkz, thmax=max(abs(x1-thoff), abs(x2-thoff))
	kz1=sph2kz_(hvmin+BE-WF, thmax,  0, Vomin)		//lower edge			//FW 0507/05
	//kz1=sph2kz_(hvmin+BE-WF, thmax,  0, Vo)		//lower edge
	kz2=sph2kz_(hvmax+BE-WF, 0,  0, Vomax)			// max at normal emission		//FW 0507/05
	//kz2=sph2kz_(hvmax+BE-WF, 0,  0, Vo)			// max at normal emission
	kz1=1E-2*round(kz1*1E2)
	kz2=1E-2*round(kz2*1E2)
	
	kxinc=ValFromList( kinc, 0, ",")
	kzinc=ValFromList( kinc, 1, ",")
	
	nkx=round((kx2-kx1)/kxinc + 1)
	nkz=round((kz2-kz1)/kzinc + 1)

	//print "kx:", kx1, kx2, kxinc, nkx, "   kz:", kz1, kz2, kzinc, nkz

// Create new image & theta array	; set BE-Kpar regular grid scaling
	if (cmpstr(outn[0], "k")==0)
		outn=img+outn
	endif	
	Make/O/N=(nkx, nkz) $outn, root:tmp:theta, root:tmp:hv
	Setscale/I x kx1, kx2, "" $outn,  root:tmp:theta, root:tmp:hv
	Setscale/I y kz1, kz2, "" $outn,  root:tmp:theta, root:tmp:hv
	//print "Created: ", outn, "(", nkx, "x", nkz, ") Ækxz=", kxinc, kzinc

// Calculate irregular grid theta values	

	//root:tmp:theta=Kp2Polar_( 3.81*(x^2+y^2)-Vo+BE,  x)+thoff			//KE=hv+BE-WF,  hv
	//root:tmp:hv=3.81*(x^2+y^2)-Vo-BE+WF								//=(KEin-Vo)-BE+WF
	//root:tmp:hv=(3.81*(x^2+y^2)-Vo_base+V0_slope*hv_base-BE+WF)/(1+V0_slope)						//FW 0507/05
	//root:tmp:theta=Kp2Polar_((3.81*(x^2+y^2)-Vo_base+V0_slope*hv_base-BE+WF)/(1+V0_slope)+BE-WF, x)+thoff			//KE=hv+BE-WF				//FW 0507/05
		
	//For linear hv dependence of Vo, the kx-kz convertion to hv-polar is solvable. For Vo=a+b*(hv-hv_0),
	// hv = (KEin-a+b*hv_0-BE+WF)/(1+b)										
	
	//$outn=interp2D( $img, root:tmp:theta , y)
	//	all in one command without intermediate array
//	$outn=interp2D( $img, Kp2Polar_((3.81*(x^2+y^2)-Vo_base+Vo_slope*hv_base-BE+WF)/(1+Vo_slope)+BE-WF,  x)+thoff , (3.81*(x^2+y^2)-Vo_base+Vo_slope*hv_base-BE+WF)/(1+Vo_slope) )  //FW 0507/05

// Without linear Vo dependence
	//$outn=interp2D($img,Kp2Polar_( 3.81*(x^2+y^2)-Vo+BE,  x)+thoff,3.81*(x^2+y^2)-Vo-BE+WF)
	$outn=interp2D($img,Kp2Polar_( 3.81*(x^2+y^2+constKpar^2)-Vo,  x)+thoff,3.81*(x^2+y^2+constKpar^2)-Vo+BE+WF)
	//$outn=interp2D( $img, Kp2Polar_( K2KE_(x, y)-Vo+BE,  x)+thoff , K2KE_(x, y)-Vo-BE+WF )
	
	
	DoWindow/F $("Kxz_"+winnam)
	if ((dopt==1)*(V_flag==0))
		// get color table info from top graph(?)
		string CTstr=ColorTableStr("", img)		//may not work if not proper top graph
		display; appendimage $outn
		
		Label bottom "Kx (1/)"
		Label left "Kz (1/)"
//		ModifyGraph height={Plan,1,left,bottom}
		ModifyGraph width={Plan,1,bottom,left}
		ModifyGraph axThick=0.5,standoff=0
		
		string titlestr="\JC"+outn+": "+num2str(nkx)+" x "+num2str(nkz)
		titlestr+="\rVo="+num2str(Vomin)+"-"+num2str(Vomax)+", WF="+num2str(WF)+",  NE="+num2str(thoff)
		Textbox/N=title/F=0/A=MT/E titlestr
		if (strlen(CTstr)>0) 
			// duplicate CT instead of using old one
			Duplicate/o $(img+"_CT") $(outn+"_CT")
			execute "ModifyImage "+outn+" cindex="+outn+"_CT"
//			execute "ModifyImage "+outn+" "+CTstr
		else
			string ctable=img+"_CT"
			if (exists(ctable)==0)
				//ctable="YellowHot"
				ModifyImage $outn ctab= {*,*,YellowHot,0}
			else
				ModifyImage $outn cindex=$ctable
			endif
		endif
		DoWindow/C $("Kxz_"+winnam)
	endif

End


Function GetVo(Vo)
//===========
// Feng Wang:  05/07/05
// Not sure why one would want a linearly varying inner potential -JD
	variable Vo
	variable Vo_base,Vo_slope,hv_base
		prompt Vo_base, "Vo_0"									  //FW 0507/05
		prompt Vo_slope, "Vo_slope"					  //FW 0507/05
		prompt hv_base, "hv_0"					  //FW 0507/05
	if(Vo<0)
		Doprompt "Vo(hv)=Vo_0 + Vo_slope*(hv-hv_0)", Vo_base,Vo_slope,hv_base
		variable/G root:tmp:Vo_base=Vo_base	//, dopt=1      //FW 0507/05
		variable/G root:tmp:Vo_slope=Vo_slope	//, dopt=1       //FW 0507/05
		variable/G root:tmp:hv_base=hv_base	//, dopt=1       //FW 0507/05
		Print "Vo="+num2str(Vo_base)+"+"+"("+num2str(Vo_slope)+")*(hv-"+num2str(hv_base)+")"
	else
		variable/G root:tmp:Vo_base=Vo	//, dopt=1      //FW 0507/05
		variable/G root:tmp:Vo_slope=0	//, dopt=1       //FW 0507/05
		variable/G root:tmp:hv_base=0	//, dopt=1       //FW 0507/05
//		Print "Vo Constant="+num2str(Vo)
	endif
End


Function volPolarBE2K(vol,WF,kxstep,newvoln)
//================
wave vol
variable WF,kxstep
string newvoln

variable x1, x2, xinc, nx, y1, y2, yinc, ny,z1,z2,zinc,nz
	//use full range of image
		xinc=DimDelta(vol, 0)
		x1=DimOffset(vol, 0)	
		x2=x1+(DimSize(vol, 0)-1)*xinc
	nx=round( (x2-x1)/xinc + 1)

		yinc=DimDelta(vol, 1)
		y1=DimOffset(vol, 1)	
		y2=y1+(DimSize(vol, 1)-1)*yinc
	ny=round( (y2-y1)/yinc + 1)

		zinc=DimDelta(vol, 2)
		z1=DimOffset(vol, 2)	
		z2=y1+(DimSize(vol, 2)-1)*yinc
		nz=dimsize(vol,2)	
	variable hvmax = max(z1,z2)

	
//	Duplicate/O/R=(x1,x2)(y1,y2) $img, root:tmp:im
	
// Calculate kx range
	variable kx1, kx2, kinc, nk, BEmax=min(y1, y2)
	////print BEmax, hv+BEmax-WF
	kx1=sph2kx_(hvmax-BEmax-WF, x1,  0)
	kx1=1E-2*round(kx1*1E2)
	kx2=sph2kx_(hvmax-BEmax-WF, x2,  0)
	kx2=1E-2*round(kx2*1E2)
	kinc=kxstep
	nk=abs(round((kx2-kx1)/kinc + 1))
	
	//print "kx:", kx1, kx2, kinc, nk

// Create new image & theta array	; set BE-Kpar regular grid scaling

	make/o/N=(nk,ny,nz) $newvoln			//, root:tmp:theta
	wave newvol =  $newvoln
	Setscale/P x kx1, kinc, "" $newvoln			//,  root:tmp:theta
	Setscale/P y y1, yinc, "" $newvoln			//,  root:tmp:theta
	Setscale/P z dimoffset(vol,2),dimdelta(vol,2), $newvoln
	
	//print "Created: ", outn, "(", nk, "x", ny, ") Æk=", kinc

// Calculate irregular grid theta values	
	//root:tmp:theta=Kp2Polar_( hv+y-WF,  x)
	//$outn=interp2D( $img, root:tmp:theta+thoff , y)
	//	all in one command without intermediate array
	newvol=interp3D( vol, Kp2Polar_(z-y-WF,  x) , y,z)

end
 
Proc PolPhotE_KxzE(voln, kvoln, V0,WF,kxinc,kzinc, transp)
//-------------
	string voln, kvoln="KxzE"
	variable V0=10,WF=4.5,kxinc=0.02,kzinc=0.05, transp
	prompt voln, "Polar-hv-E volume to convertt", popup, WaveList("!*ct", ";", "DIMS:3")
	prompt kvoln, "Output volume name or extension (kxzE)"
	prompt V0 "Inner Potential (eV)"
	prompt WF "Work Function (eV)"
	prompt kxinc, "kx-increment"
	prompt kzinc, "kz-increment"
	prompt transp, "Input axes order", popup, "Pol-Phot-E;Pol-E-Phot"
//	prompt kinc, "k-increments:  kx,kz"
	
//	string voln=NameOfWave(vol)
	
	string inopt="/V0="+num2str(V0)+"/WF="+num2str(WF)
	string outopt="/KINC="+num2str(kxinc)+","+num2str(kzinc)
	if (stringmatch(kvoln,"KxzE"))
		kvoln=voln+"_KxzE"
	else
		outopt+="/D="+kvoln
	endif
	
	if (transp==2)
		string cmd1="VolTranspose("+voln+",  \"/XZY/D="+voln+"r\")"		
		print cmd1
		execute cmd1
		voln=voln+"r"
	endif

	string cmd="PolPhotE2KxzE("+voln+",\""+inopt+"\",\""+outopt+"\")"
	print cmd
	execute cmd
End

//Function PolhvE2KxzE(vol,V0,WF,kxstep,kzstep,newvoln)
Function PolPhotE2KxzE(vol,inopt,outopt)
// ==================
//  Input volume:  (Polar-hv-BE)
//       x = Detector Angle (polar)
//       y = Photon Energy (hv)
//       z = Binding energy  (negative below EF)
//   Conversion parameters:
//       V0 = inner potential
//       WF = work function
//  Output info:
//      kxinc, kzinc    (BEstep=same as original)
//      newvoln  --> KxzE  (Kx -Kz - BE)
	wave vol
	string inopt, outopt
	
	//energy correction:
	variable V0=10, WF=4.5
		V0=KeyVal("V0",inopt)
		V0=SelectNumber(numtype(V0)==2, V0, 10)		// use default if NAN
		WF=KeyVal("WF",inopt)
		WF=SelectNumber(numtype(WF)==2, WF, 4.5)		// use default if NAN
	
	// output k-increments
	variable kxinc=0.02, kzinc=0.05
	if (KeySet("KINC", outopt))
		kxinc=str2num(StringFromList(0, KeyStr("KINC", ",")) )
		kxinc=SelectNumber(numtype(kxinc)==2, kxinc, 0.02)		// use default if NAN
		kzinc=str2num(StringFromList(1, KeyStr("KINC", ",")) )
		kzinc=SelectNumber(numtype(kzinc)==2, kzinc, kxinc)    //NaN for single value list
	endif
	
	string voln=NameOfWave(vol), newvoln
	newvoln=KeyStr("D",outopt)
	if ((strlen(newvoln)==0) + stringmatch(newvoln,"KxzE"))
		newvoln=voln+"_KxzE"
	endif
	
	NewDataFolder/O root:tmp
	
	variable nx=dimsize(vol,0)
	variable ny=dimsize(vol,1)
	variable nz=dimsize(vol,2)
	
	variable angmin = dimoffset(vol,0)
	variable angmax =angmin+dimdelta(vol,0)*(nx-1)
	variable y1 = dimoffset(vol,1)
	variable y2 =y1+dimdelta(vol,1)*(ny-1)
	variable e1=dimoffset(vol,2)
	variable e2=e1+dimdelta(vol,2)*(nz-1)
	variable BEmin = min(e1,e2)			// max negative BE
	variable BEmax = max(e1,e2)			// positive BE above EF
	//calculate k grids
	
// Calculate kx range
	variable hvmax=max(y1, y2), hvmin=min(y1,y2)
	variable kx1, kx2, nkx
	kx1=min( sph2kx_(hvmax-BEmax-WF, angmin,  0), sph2kx_(hvmin-BEmax-WF, angmin,  0) )
	kx2=max( sph2kx_(hvmax-BEmax-WF, angmax,  0), sph2kx_(hvmin-BEmax-WF, angmax,  0) )
	kx1=1E-2*round(kx1*1E2)
	kx2=1E-2*round(kx2*1E2)	
	
// Calculate kz range
	variable kz1, kz2, nkz, thmax=max(abs(angmin), abs(angmax))
	kz1=sph2kz_(hvmin+BEmin-WF, thmax,  0, V0)		//lower edge
	kz2=sph2kz_(hvmax+BEmax-WF, 0,  0, V0)			// max at normal emission
	kz1=1E-2*round(kz1*1E2)
	kz2=1E-2*round(kz2*1E2)
	print "BEmin.max=",BEmin, BEmax, "Kz1,2=",kx1,kz2
	
	nkx=round((kx2-kx1)/kxinc + 1)
	nkz=round((kz2-kz1)/kzinc + 1)
	
	// Create new image & theta array	; set BE-Kpar regular grid scaling
	Make/O/N=(nkx, nkz,nz) $newvoln//, root:tmp:theta, root:tmp:hv
	Setscale/I x kx1, kx2, $newvoln//,  root:tmp:theta, root:tmp:hv
	Setscale/I y kz1, kz2, $newvoln//,  root:tmp:theta, root:tmp:hv
	Setscale/P z dimoffset(vol,2),dimdelta(vol,2), $newvoln 			//energy
	//print "Created: ", outn, "(", nkx, "x", nkz, ") Ækxz=", kxinc, kzinc
	wave newvol = $newvoln

// Calculate irregular grid theta values	
	//root:tmp:theta=Kp2Polar_( 3.81*(x^2+y^2)-Vo-z,  x)			//KE=hv+BE-WF,  hv
	//root:tmp:hv=3.81*(x^2+y^2)-Vo-BE+WF								//=(KEin-Vo)-BE+WF
	newvol=interp3D( vol, Kp2Polar_( 3.81*(x^2+y^2)-V0-z,  x) , 3.81*(x^2+y^2)-V0-z+WF,z)
	
	//$outn=interp2D( $img, Kp2Polar_( K2KE_(x, y)-Vo+BE,  x)+thoff , K2KE_(x, y)-Vo-BE+WF )
	
End

Proc PolTiltE_KxyE(voln, kvoln, hv, WF,kinc,angoff, geometry, transp)
//-------------
	string voln, kvoln="KxyE"
	variable hv=100, WF=4.5,transp
	string kinc=StrVarOrDefault( "root:kcvt:kincr", ".02,.02")
	string angoff=StrVarOrDefault( "root:kcvt:ang_off", "0,0")
	variable geometry=NumVarOrDefault( "root:kcvt:geometry", 3)
	prompt voln, "Polar-Elevation-E volume to convertt", popup, WaveList("!*ct", ";", "DIMS:3")
	prompt kvoln, "Output volume name or extension (kxyE)"
	prompt hv "Photon Energy (eV), (KE=hv-WF+BE)"
	prompt WF "Work Function (eV)"
	prompt kinc, "k-increments:  kx,ky"
	prompt angoff, "angle offsets:  pol0, elev0"
	prompt geometry "Geometry", popup, "Detector Polar-Elev;Sample Ang+Horiz Slit;Sample Ang+Vert Slit"
	prompt transp, "Input axes order", popup, "Pol-Elev-E;Pol-E-Elev"
//	prompt kinc, "k-increments:  kx,kz"
	
//	string voln=NameOfWave(vol)
	NewDataFolder/O/S root:kcvt
		string/G outkxyE=kvoln, kincr=kinc, ang_off=angoff
//		variable/G ken=KE, azim_off=azimoff
		variable/G geometry=geometry
//		variable/G root:kcvt:numkx=kxopt, root:kcvt:numky=kyopt
		//variable/G root:kcvt:poloff=pol_off, root:kcvt:elvoff=elv_off
	SetDataFolder root:
	
	string inopt="/hv="+num2str(hv)+"/WF="+num2str(WF)
	inopt+="/AOFF="+angoff
	inopt+="/GEO="+num2str(geometry)
	
	string outopt="/KINC="+	kinc					//num2str(kxinc)+","+num2str(kyinc)
	if (stringmatch(kvoln,"KxyE"))
		kvoln=voln+"_KxyE"
	else
		outopt+="/D="+kvoln
	endif
	
	if (transp==2)
		string cmd1="VolTranspose("+voln+",  \"/XZY/D="+voln+"r\")"		
		print cmd1
		execute cmd1
		voln=voln+"r"
	endif

	string cmd
	cmd="PolTiltE2KxyE("+voln+",\""+inopt+"\",\""+outopt+"\")"
//	cmd="PolTiltE2KxyE_T("+voln+",\""+inopt+"\",\""+outopt+"\")"
	print cmd
	execute cmd
End

Function PolTiltE2KxyE(vol,inopt,outopt)
// ==================
//  Input volume:  (Polar-hv-BE)
//       x = Detector Angle (polar)
//       y = Elevation Angle (tilt, beta)
//       z = Binding energy  (negative below EF)
//   Conversion parameters:
//       KE = kinetic energy AT EF   (or KE_EF=hv-WF)
//       hv = photon energy
//       WF = work function
//       AOFF = angle offsets (pol0, tilt0)
//	     GEO = 1,2,3    (1=Det Polar Elev, 2= Sample Poler-tilt+Horiz Slit, 3=Sample Polar-Tilt+Vert Slit
//  Output info:
//      kxinc, kyinc    (BEstep=same as original)
//      newvoln  --> KxyE  (Kx -Ky - BE)
	wave vol
	string inopt, outopt
	
	//energy correction:
	variable hv=10, WF=4.5
		hv=KeyVal("hv",inopt)
		hv=SelectNumber(numtype(hv)==2, hv, 100)		// use default if NAN
		WF=KeyVal("WF",inopt)
		WF=SelectNumber(numtype(WF)==2, WF, 4.5)		// use default if NAN
	variable KE_EF=hv-WF
	
	//angle offsets:
	variable pol0=0, tilt0=0
	string tmpstr=KeyStr("AOFF", inopt)
	if (KeySet("AOFF", inopt))
		pol0=str2num(StringFromList(0, tmpstr, ",")) 
		tilt0=str2num(StringFromList(1, tmpstr, ",")) 
	endif
	print KeyStr("AOFF", inopt), pol0, tilt0
	
	// output k-increments
	variable kxinc=0.02, kyinc=0.02
	if (KeySet("KINC", outopt))
		kxinc=str2num(StringFromList(0, KeyStr("KINC", outopt), ",")) 
		print kxinc, numtype(kxinc), KeyStr("KINC", ",")
		kxinc=SelectNumber(numtype(kxinc)==2, kxinc, 0.02)		// use default if NAN
		kyinc=str2num(StringFromList(1, KeyStr("KINC", outopt), ",")) 
		kyinc=SelectNumber(numtype(kyinc)==2, kyinc, kxinc)    //NaN for single value list
	endif
	
	string voln=NameOfWave(vol), newvoln
	newvoln=KeyStr("D",outopt)
	if ((strlen(newvoln)==0) + stringmatch(newvoln,"KxyE"))
		newvoln=voln+"_KxyE"
	endif
	
	//NewDataFolder/O root:tmp
	
	
	// Determine input ranges
	//  Kcvt_InputRange(img, inopt)  //--> variables in folder
	//  
	string xrng=KeyStr("X", inopt), yrng=KeyStr("Y",inopt)
	variable  x1, x2, xinc, nx, xmid
	if (strlen(xrng)>0)
		x1=ValFromList(xrng, 0,",")		//specified subrange
		x2=ValFromList(xrng, 1,",")
//		xinc=ValFromList(xrng, 2,",")
		xinc=DimDelta(vol, 0)
	else					//default use full range of image
		xinc=DimDelta(vol, 0)
		x1=DimOffset(vol, 0)	
		x2=x1+(DimSize(vol, 0)-1)*xinc
	endif
	nx=round( (x2-x1)/xinc + 1)
	xmid = (x1+x2)/2
//	print "x:", x1, x2, xmid, xinc, nx
	
	variable  y1, y2, yinc, ny, ymid
	if (strlen(yrng)>0)
		y1=ValFromList(yrng, 0,",")
		y2=ValFromList(yrng, 1,",")
//		yinc=ValFromList(yrng, 2,",")
		yinc=DimDelta(vol, 1)
	else					//default use full range of image
		yinc=DimDelta(vol, 1)
		y1=DimOffset(vol, 1)	
		y2=y1+(DimSize(vol, 1)-1)*yinc
	endif
	ny=round( (y2-y1)/yinc + 1)
	ymid = (y1+y2)/2
//	print "y:", y1, y2, ymid, yinc, ny
	
	variable ne=dimsize(vol,2)
	variable einc=dimdelta(vol,2)
	variable e1=dimoffset(vol,2)
	variable e2=e1+einc*(ne-1)
	variable BEmin = min(e1,e2)			// max negative BE
	variable BEmax = max(e1,e2)			// positive BE above EF
	
	// use x1,x2,xmid (& y) to populate range arrays for plotting purposes,
	//     k-range conversion & k-range evaluation
	make/o/n=9 xrngw, yrngw
	xrngw={x2, x2, xmid, x1, x1, x1, xmid, x2, x2} 
	yrngw={ymid, y2, y2, y2, ymid, y1, y1, y1, ymid}
// Add pol0 and tilt0 offsets
	xrngw-=pol0
	yrngw-=tilt0
	
	// Geometry formula selection
	variable geometry=2
	if (KeySet("GEO", inopt))
		geometry=KeyVal("GEO",inopt)
	endif
	
	//Convert xy-range (angles) to kxy
	// pol0 & tilt0 offests NOT included in this range conversion
	make/o/n=9 kxrngw, kyrngw
	kxrngw=DualA2kx_( KE_EF, xrngw, yrngw, geometry)
	kyrngw=DualA2ky_( KE_EF, xrngw, yrngw, geometry)
	
	//Determine kmin, kmax values & round
	variable kxmin, kxmax, kymin, kymax
	WaveStats/Q kxrngw
	kxmin=1E-2*round(V_min*1E2)
	kxmax=1E-2*round(V_max*1E2)
	WaveStats/Q kyrngw
	kymin=1E-2*round(V_min*1E2)
	kymax=1E-2*round(V_max*1E2)
	
	variable nkx, nky	
	nkx=round((kxmax-kxmin)/kxinc + 1)
	nky=round((kymax-kymin)/kyinc + 1)
	print "kx:", kxmin, kxmax, kxinc, nkx, "   ky:", kymin, kymax, kyinc, nky
	
	// Create new image & theta array	; set BE-Kpar regular grid scaling
	Make/O/N=(nkx, nky,ne) $newvoln	, root:kcvt:polEw, root:kcvt:tiltEw
	WAVE newvol = $newvoln
	WAVE polEw=root:kcvt:polEw, tiltEw=root:kcvt:tiltEw
	Setscale/I x kxmin, kxmax, "" $newvoln, polEw, tiltEw
	Setscale/I y kymin, kymax, "" $newvoln, polEw, tiltEw
	Setscale/P z e1, einc, $newvoln, polEw, tiltEw
	//print "Created: ", outn, "(", nkx, "x", nky, "x", nky, ") ÆkxzE=", kxinc, kzinc, nz



	variable r2d=180/pi		//, ktot=0.5123*sqrt(KE_EF)	
	make/o/n=(ne) root:kcvt:koutw
	WAVE koutw=root:kcvt:koutw
	Setscale/P x e1, einc, koutw
	koutw=0.5123*sqrt(hv-WF+x)			//x=negative BE

	variable tref=StartMSTimer	
	
	//root:kcvt:polw=r2d*asin(x/sqrt(ktot^2-y^2) )
	//root:kcvt:tiltw=r2d*asin(y/ktot)
//	polw=r2d*asin(xrot(x, y, azoff)/sqrt(ktot^2-yrot(x, y, azoff)^2) )
//	tiltw=r2d*asin(yrot(x, y, azoff)/ktot)
// Calculate irregular grid theta values	
	//root:tmp:theta=Kp2Polar_( 3.81*(x^2+y^2)-Vo-z,  x)			//KE=hv+BE-WF,  hv
	//root:tmp:hv=3.81*(x^2+y^2)-Vo-BE+WF								//=(KEin-Vo)-BE+WF
//	newvol=interp3D( vol, Kp2Polar_( 3.81*(x^2+y^2)-V0-z,  x) , 3.81*(x^2+y^2)-V0-z+WF,z)


//	newvol=interp3D( vol, r2d*asin( x/sqrt(koutw(z)^2-y^2) ) , r2d*asin(y/koutw(z)) , z )
	//SWAP POLAR & ELEVATION ANGLES
//	newvol=interp3D( vol,  r2d*asin(x/koutw(z)), r2d*asin( y/sqrt(koutw(z)^2-x^2) ) ,  z )		
//	$outn=interp2D( $img, r2d*asin( x/sqrt(kout^2-y^2) )+pol0 , r2d*asin(y/kout)+elev0 )
	
//		kout=interp2D( img, polw+xoff , tiltw+yoff )
//	kout=interp2D( img, r2d*asin( xrot(x, y, azoff)/sqrt(ktot^2-yrot(x, y, azoff)^2) )+xoff , r2d*asin(yrot(x, y, azoff)/ktot)+yoff )
	//$outn=interp2D( $img, r2d*asin( x/sqrt(kout^2-y^2) )+pol0 , r2d*asin(y/kout)+elev0 )
	//$outn=interp2D( $img, Kp2Polar_( K2KE_(x, y)-Vo+BE,  x)+pol0 , K2KE_(x, y)-Vo-BE+WF )
	
	//$outn=interp2D( $img, Kp2Polar_( K2KE_(x, y)-Vo+BE,  x)+thoff , K2KE_(x, y)-Vo-BE+WF )

//	print stopMStimer(tref)/1E6, " sec"
//	abort
	
	IF (geometry==1)				// Single channel Detector Polar-Elevation +  fixed sample

		polEw=r2d*asin( x/sqrt(koutw(z)^2-y^2) )
		tiltEw=r2d*asin( y/koutw(z) )

		newvol=interp3D( vol, polEw + pol0 , tiltEw + tilt0, z )
		// all-in-one command w/o intermediate arrays		
//		$outn=interp2D( $img, r2d*asin( xrot(x, y, ang)/sqrt(kout^2-yrot(x, y, ang)^2) )+pol0 , r2d*asin(yrot(x, y, ang)/kout)+elev0 )	

	ELSE
	IF (geometry==2)			// Sample Polar-Tilt +  Detector Horizonatal Slit  (polw=detector angle)

		polEw=r2d*asin(x/koutw(z))
		tiltEw=r2d*asin(y/sqrt(koutw(z)^2-x^2) )

		newvol=interp3D( vol, polEw + pol0 , tiltEw + tilt0, z)
		//	all in one command without intermediate array
//		$outn=interp2D( $img, r2d*asin(xrot(x, y, ang)/kout)+pol0, r2d*asin( yrot(x, y, ang)/sqrt(kout^2-xrot(x, y, ang)^2) )+elev0  )
	
	ELSE
	IF (geometry==3)			// Sample Polar-Tilt +  Detector Vertical Slit  (polw=detector angle)

		tiltEw=r2d*asin(y/koutw(z)*cos(tilt0/r2d) + sqrt(koutw(z)^2-x^2-y^2)/koutw(z)*sin(tilt0/r2d) )
		polEw=r2d*asin(x/koutw(z)/cos(tiltEw/r2d) )

		newvol=interp3D( vol, polEw + pol0 , tiltEw, z )			// tilt0 offset included above
	
	ENDIF
	ENDIF
	ENDIF
	
	
	print stopMStimer(tref)/1E6, " sec"
	
End



Function PolElevE2KxyE_T(vol,inopt,outopt)
// ==================
//  _T = Triangulation method
//  Input volume:  (Polar-hv-BE)
//       x = Detector Angle (polar)
//       y = Elevation Angle (tilt, beta)
//       z = Binding energy  (negative below EF)
//   Conversion parameters:
//       KE = kinetic energy AT EF   (or KE_EF=hv-WF)
//       hv = photon energy
//       WF = work function
//  Output info:
//      kxinc, kyinc    (BEstep=same as original)
//      newvoln  --> KxyE  (Kx -Ky - BE)
	wave vol
	string inopt, outopt
	
	//energy correction:
	variable hv=10, WF=4.5
		hv=KeyVal("hv",inopt)
		hv=SelectNumber(numtype(hv)==2, hv, 100)		// use default if NAN
		WF=KeyVal("WF",inopt)
		WF=SelectNumber(numtype(WF)==2, WF, 4.5)		// use default if NAN
	variable KE_EF=hv-WF
	
	// output k-increments
	variable kxinc=0.02, kyinc=0.02
	if (KeySet("KINC", outopt))
		kxinc=str2num(StringFromList(0, KeyStr("KINC", outopt), ",")) 
		print kxinc, numtype(kxinc), KeyStr("KINC", ",")
		kxinc=SelectNumber(numtype(kxinc)==2, kxinc, 0.02)		// use default if NAN
		kyinc=str2num(StringFromList(1, KeyStr("KINC", outopt), ",")) 
		kyinc=SelectNumber(numtype(kyinc)==2, kyinc, kxinc)    //NaN for single value list
	endif
	
	string voln=NameOfWave(vol), newvoln
	newvoln=KeyStr("D",outopt)
	if ((strlen(newvoln)==0) + stringmatch(newvoln,"KxyE"))
		newvoln=voln+"_KxyE"
	endif
	
	//NewDataFolder/O root:tmp
	
	
	// Determine input ranges
	//  Kcvt_InputRange(img, inopt)  //--> variables in folder
	//  
	string xrng=KeyStr("X", inopt), yrng=KeyStr("Y",inopt)
	variable  x1, x2, xinc, nx, xmid
	if (strlen(xrng)>0)
		x1=ValFromList(xrng, 0,",")		//specified subrange
		x2=ValFromList(xrng, 1,",")
//		xinc=ValFromList(xrng, 2,",")
		xinc=DimDelta(vol, 0)
	else					//default use full range of image
		xinc=DimDelta(vol, 0)
		x1=DimOffset(vol, 0)	
		x2=x1+(DimSize(vol, 0)-1)*xinc
	endif
	nx=round( (x2-x1)/xinc + 1)
	xmid = (x1+x2)/2
//	print "x:", x1, x2, xmid, xinc, nx
	
	variable  y1, y2, yinc, ny, ymid
	if (strlen(yrng)>0)
		y1=ValFromList(yrng, 0,",")
		y2=ValFromList(yrng, 1,",")
//		yinc=ValFromList(yrng, 2,",")
		yinc=DimDelta(vol, 1)
	else					//default use full range of image
		yinc=DimDelta(vol, 1)
		y1=DimOffset(vol, 1)	
		y2=y1+(DimSize(vol, 1)-1)*yinc
	endif
	ny=round( (y2-y1)/yinc + 1)
	ymid = (y1+y2)/2
//	print "y:", y1, y2, ymid, yinc, ny
	
	variable ne=dimsize(vol,2)
	variable einc=dimdelta(vol,2)
	variable e1=dimoffset(vol,2)
	variable e2=e1+einc*(ne-1)
	variable BEmin = min(e1,e2)			// max negative BE
	variable BEmax = max(e1,e2)			// positive BE above EF
	
	// use x1,x2,xmid (& y) to populate range arrays for plotting purposes,
	//     k-range conversion & k-range evaluation
	make/o/n=9 xrngw, yrngw
	xrngw={x2, x2, xmid, x1, x1, x1, xmid, x2, x2} 
	yrngw={ymid, y2, y2, y2, ymid, y1, y1, y1, ymid}
	
	// Geometry formula selection
	variable geometry=3
	if (KeySet("GEO", inopt))
		geometry=KeyVal("GEO",inopt)
	endif
	
	//Convert xy-range (angles) to kxy
	make/o/n=9 kxrngw, kyrngw
	kxrngw=DualA2kx_( KE_EF, xrngw, yrngw, geometry)
	kyrngw=DualA2ky_( KE_EF, xrngw, yrngw, geometry)
	
	//Determine kmin, kmax values & round
	variable kxmin, kxmax, kymin, kymax
	WaveStats/Q kxrngw
	kxmin=1E-2*round(V_min*1E2)
	kxmax=1E-2*round(V_max*1E2)
	WaveStats/Q kyrngw
	kymin=1E-2*round(V_min*1E2)
	kymax=1E-2*round(V_max*1E2)
	
	variable nkx, nky	
	nkx=round((kxmax-kxmin)/kxinc + 1)
	nky=round((kymax-kymin)/kyinc + 1)
	print "kx:", kxmin, kxmax, kxinc, nkx, "   ky:", kymin, kymax, kyinc, nky
	
	// Create new image & theta array	; set BE-Kpar regular grid scaling
	Make/O/N=(nkx, nky,ne) $newvoln				//, root:tmp:theta, root:tmp:elev
	Setscale/I x kxmin, kxmax, "" $newvoln		//, polw, tiltw
	Setscale/I y kymin, kymax, "" $newvoln		//, polw, tiltw
	Setscale/P z e1, einc, $newvoln 					//energy
	//print "Created: ", outn, "(", nkx, "x", nky, "x", nky, ") ÆkxzE=", kxinc, kzinc, nz
	wave newvol = $newvoln


	variable r2d=180/pi, d2r=pi/180		//, ktot=0.5123*sqrt(KE_EF)	
	make/o/n=(ne) koutw
	WAVE koutw=koutw
	Setscale/P x e1, einc, koutw
	koutw=0.5123*sqrt(hv-WF+x)			//x=negative BE
	
	
	//Create 4 vector for interpolation
	duplicate/o vol $(voln+"4v")
	WAVE vol4v = $(voln+"4v")
	Redimension/n=(-1,-1,-1, 4) vol4v
//	vol4v[][][][0]=x1+xinc*p			//x = polar = delta - theta
//	vol4v[][][][1]=y1+yinc*q			//y= elev  = beta  
//  Put different k-conversions here

	variable beamline=4
	if (beamline==7)
//  *** BL7:  Polar-Tilt Sample Goniometer + Scient Multi-angle detection w/ Horiz. Slit
//      kx = k_out*sin( delta - theta)		// delta=det. ang., theta=fixed angle = 0 (typically)
//     ky = k_out*cos(delta-theta)*sin(beta)    // beta = sample tilt/flip angle
		vol4v[][][][0]=koutw[e1+einc*r]*sin( d2r* (x1+xinc*p) )			//kx
		vol4v[][][][1]=koutw[e1+einc*r]*cos( d2r* (x1+xinc*p) )*sin(d2r* (y1+yinc*q) )			//ky
	else	
//  *** BL403:  Polar-Tilt Sample Goniometer + Scient Multi-angle detection w/ Vertical Slit
//     kx = k_out*cos( delta ) sin( theta)		// x: theta=sample polar angle, y: delta=det. ang
//     ky = k_out*(  cos(delta) cos(theta)*sin(beta) + sin(delta)cos(beta)	 )    //beta=fixed tilt angle=0
	variable beta_ang=0
		vol4v[][][][0]=koutw[e1+einc*r]*cos( d2r* (y1+yinc*p) )*sin( d2r* (x1+xinc*p) )			//kx
		vol4v[][][][1]=koutw[e1+einc*r]*( cos( d2r* (y1+yinc*p) )*sin(d2r* (x1+xinc*q) )	*sin(d2r* beta_ang ) + sin(d2r* (y1+yinc*q) )*cos(d2r* beta_ang ) )			//ky
	endif	
	vol4v[][][][2]=e1+einc*r		// no transformation of BE
 	vol4v[][][][3]=vol[p][q][r]	
	redimension/n=(nx*ny*ne*4) vol4v
	redimension/n=(nx*ny*ne, 4) vol4v
	
	vol4v+=gnoise(0.0001)
	
	variable tref=StartMSTimer
	
	Triangulate3D/OUT=1 vol4v
	print "Triangulate: ", stopMStimer(tref)/1E6, " sec"
	tref=StartMSTimer
	
	variable method=2
	if (method==1)
		newvol = interp3D( Data4col, x, y, z, M_3DVertexList)
	else
		 Interpolate3D /RNGX={kxmin, kxinc, nkx}/RNGY={kymin, kyinc, nky}/RNGZ={e1, einc, ne}/DEST=$newvoln triangulationWave=M_3DVertexList, srcWave=vol4v
	
	endif
	//root:kcvt:polw=r2d*asin(x/sqrt(ktot^2-y^2) )
	//root:kcvt:tiltw=r2d*asin(y/ktot)
//	polw=r2d*asin(xrot(x, y, azoff)/sqrt(ktot^2-yrot(x, y, azoff)^2) )
//	tiltw=r2d*asin(yrot(x, y, azoff)/ktot)
// Calculate irregular grid theta values	
	//root:tmp:theta=Kp2Polar_( 3.81*(x^2+y^2)-Vo-z,  x)			//KE=hv+BE-WF,  hv
	//root:tmp:hv=3.81*(x^2+y^2)-Vo-BE+WF								//=(KEin-Vo)-BE+WF
//	newvol=interp3D( vol, Kp2Polar_( 3.81*(x^2+y^2)-V0-z,  x) , 3.81*(x^2+y^2)-V0-z+WF,z)


//	newvol=interp3D( vol, r2d*asin( x/sqrt(koutw(z)^2-y^2) ) , r2d*asin(y/koutw(z)) , z )
	//SWAP POLAR & ELEVATION ANGLES
//	newvol=interp3D( vol,  r2d*asin(x/koutw(z)), r2d*asin( y/sqrt(koutw(z)^2-x^2) ) ,  z )		
//	$outn=interp2D( $img, r2d*asin( x/sqrt(kout^2-y^2) )+pol0 , r2d*asin(y/kout)+elev0 )
	
//		kout=interp2D( img, polw+xoff , tiltw+yoff )
//	kout=interp2D( img, r2d*asin( xrot(x, y, azoff)/sqrt(ktot^2-yrot(x, y, azoff)^2) )+xoff , r2d*asin(yrot(x, y, azoff)/ktot)+yoff )
	//$outn=interp2D( $img, r2d*asin( x/sqrt(kout^2-y^2) )+pol0 , r2d*asin(y/kout)+elev0 )
	//$outn=interp2D( $img, Kp2Polar_( K2KE_(x, y)-Vo+BE,  x)+pol0 , K2KE_(x, y)-Vo-BE+WF )
	
	//$outn=interp2D( $img, Kp2Polar_( K2KE_(x, y)-Vo+BE,  x)+thoff , K2KE_(x, y)-Vo-BE+WF )
	print "Interpolate3d:", stopMStimer(tref)/1E6, " sec"
	
End

	
	
	
Function PolElevhv2Kxyhv( vol,kxinc,kyinc,azimoff)
//============
// Does 2D k-conversion in z-loop, e.g. stack of photon dependent FS maps
//
//  Input volume:  (Polar-hv-BE)
//       x = Detector Angle (polar)
//       y = Elevation angle (tilt, beta)
//       z = Photon  (stack of FS's)
//   Conversion parameters:
//       V0 = inner potential
//       WF = work function
//  Output info:
//      kxinc, kyinc   (z_hv=same as original)
//      newvoln  --> Kxzhv  (Kx -Kz - hv)


	wave vol
	variable kxinc,kyinc
	variable azimoff
	
	variable xsize=dimsize(vol,0),ysize=dimsize(vol,1),zsize=dimsize(vol,2)
	variable xoffset=dimoffset(vol,0),yoffset=dimoffset(vol,1),zoffset=dimoffset(vol,2)
	variable xdelta=dimdelta(vol,0),ydelta=dimdelta(vol,1),zdelta=dimdelta(vol,2)
	
	variable hv,hvmax,hvmin
	variable WF=4.5
	if(zdelta>0)
	hvmin=zoffset
	hvmax=zoffset+zdelta*(zsize-1)
	else
	hvmax=zoffset
	hvmin=zoffset+zdelta*(zsize-1)	
	endif
	
	// Geometry formula selection
	variable geometry=2
//	if (KeySet("GEO", inopt))
//		geometry=KeyVal("GEO",inopt)
//	endif
	
// Calculate k-grid range
	variable ang=-azimoff
	variable KE=hvmax-WF
	
	variable kx1, kx2, kx_1, kx_2,nkx
	variable ky1, ky2, ky_1, ky_2,nky
		//convert corners
	kx1=dualA2kx_(KE, xoffset,  0, geometry);  
	kx2=dualA2kx_(KE, xoffset+xdelta*(xsize-1),  0, geometry)	   // max at ky=0
	ky1=dualA2ky_(KE, 0, yoffset, geometry);
	ky2=dualA2ky_(KE, 0, ydelta*(ysize-1), geometry)    //straight lines in kspace, any kx ok
		//sort min/max
	kx_1=min(kx1, kx2); kx_2=max(kx1, kx2)
	ky_1=min(ky1, ky2); ky_2=max(ky1, ky2)
		//rotate
	kx1= min( xrot(kx_1, ky_1, -ang),  xrot(kx_1, ky_2, -ang) )
	kx2=max( xrot(kx_2, ky_1, -ang),  xrot(kx_2, ky_2, -ang) )
	ky1= min( yrot(kx_1, ky_1, -ang),  yrot(kx_2, ky_1, -ang) )
	ky2=max( yrot(kx_1, ky_2, -ang),  yrot(kx_2, ky_2, -ang) )
		//round to simpler numbers (2 decimeal)
	kx1=1E-2*round(kx1*1E2)
	kx2=1E-2*round(kx2*1E2)
	ky1=1E-2*round(ky1*1E2)
	ky2=1E-2*round(ky2*1E2)
	
	nkx=round((kx2-kx1)/kxinc + 1)
	nky=round((ky2-ky1)/kyinc + 1)

	//print "kx:", kx1, kx2, kxinc, nkx, "   ky:", ky1, ky2, kyinc, nky

// Create new image & polar/elev arrays	; set Kx-Ky regular grid scaling
	Make/O/N=(nkx, nky,zsize) volkxymap, root:tmp:polw, root:tmp:tiltw
	Setscale/I x kx1, kx2, "" volkxymap,  root:tmp:polw, root:tmp:tiltw
	Setscale/I y ky1, ky2, "" volkxymap,  root:tmp:polw, root:tmp:tiltw
	Setscale/p z,zoffset,zdelta,volkxymap,  root:tmp:polw, root:tmp:tiltw 
	//print "Created: ", outn, "(", nkx, "x", nky, ") Ækxy=", kxinc, kyinc

// Calculate irregular grid polar& elev values
	// worry about misorientation correction

	variable r2d=180/pi
	wave polw=root:tmp:polw
	wave tiltw=root:tmp:tiltw
	make/o/n=(zsize) kout=0.5123*sqrt(zoffset+(p-1)*zdelta-WF)
	//root:tmp:polw=r2d*asin(x/sqrt(kout^2-y^2) )
	//root:tmp:tiltw=r2d*asin(y/kout)
	polw=r2d*asin(xrot(x, y, ang)/sqrt(kout[r]^2-yrot(x, y, ang)^2) )
	tiltw=r2d*asin(yrot(x, y, ang)/kout[r])
	//$outn=interp2D( $img, root:tmp:theta , y)
	//	all in one command without intermediate array
	make/o/n=(xsize,ysize) root:tmp:tempimage
	make/o/n=(nkx,nky) root:tmp:tempkimage
	wave tempimage=root:tmp:tempimage
	wave tempkimage=root:tmp:tempkimage
	Setscale/I x kx1, kx2, tempkimage
	Setscale/I y ky1, ky2, tempkimage
	Setscale/p x xoffset, xdelta, tempimage
	Setscale/p y yoffset, ydelta, tempimage
	
	variable ii
	For(ii=0;ii<zsize;ii+=1)			// 2D k-conversion; loop over z
		tempimage = vol[p][q][ii]
		tempkimage=interp2D( tempimage, r2d*asin( xrot(x, y, ang)/sqrt(kout[ii]^2-yrot(x, y, ang)^2) ) , r2d*asin(yrot(x, y, ang)/kout[ii]))
		volkxymap[][][ii]=tempkimage(x)(y)
	endfor

End


Function PolElevhv2Kxyhv3( vol,V0,kxinc,kyinc,kzinc,azimoff)
//=================
//  Input volume:  (Polar-Elev-BE)
//       x = Detector Angle (polar)
//       y = Elevation angle (tilt, beta)
//       z = Binding Energy
//   Conversion parameters:
//       V0 = inner potential
//       WF = work function
//      ?? Geometry
//  Output info:
//      kxinc, kyinc   (z_hv=same as original)
//      newvoln  --> Kxzhv  (Kx -Kz - hv)
	wave vol
	variable kxinc,kyinc,kzinc
	variable V0
	variable azimoff
	
	variable xsize=dimsize(vol,0),ysize=dimsize(vol,1),zsize=dimsize(vol,2)
	variable xoffset=dimoffset(vol,0),yoffset=dimoffset(vol,1),zoffset=dimoffset(vol,2)
	variable xdelta=dimdelta(vol,0),ydelta=dimdelta(vol,1),zdelta=dimdelta(vol,2)
	variable xmax=max(abs(xoffset),abs(xoffset+xdelta*(xsize-1)))
	variable ymax=max(abs(yoffset),abs(yoffset+ydelta*(ysize-1)))
	variable hv,hvmax,hvmin
	variable WF=4.5
	if(zdelta>0)
	hvmin=zoffset
	hvmax=zoffset+zdelta*(zsize-1)
	else
	hvmax=zoffset
	hvmin=zoffset+zdelta*(zsize-1)	
	endif
	print "x:", xoffset,xdelta,xsize
	print "y:",yoffset,ydelta,ysize
	print "z:",zoffset,zdelta,zsize
	print "hvmin" ,hvmin
	print "hvmax", hvmax
	
// Calculate k-grid range
	variable ang=-azimoff
	variable KE=hvmax-WF
	
	variable kx1, kx2, kx_1, kx_2,nkx
	variable ky1, ky2, ky_1, ky_2,nky
	variable kz1,kz2,nkz
	
	variable geometry=2
		//convert corners
	kx1=dualA2kx_(KE, xoffset,  0, geometry);  
	kx2=dualA2kx_(KE, xoffset+xdelta*(xsize-1),  0, geometry)	   // max at ky=0
	ky1=dualA2ky_(KE, 0, yoffset, geometry); 
	ky2=dualA2ky_(KE, 0, yoffset+ydelta*(ysize-1), geometry)    //straight lines in kspace, any kx ok
	kz2=sph2kz_(hvmax-WF, 0, 0, V0 );			// Photon stack
	kz1=sph2kz_( hvmin-WF, xmax, ymax,V0)
		//sort min/max
	kx_1=min(kx1, kx2); kx_2=max(kx1, kx2)
	ky_1=min(ky1, ky2); ky_2=max(ky1, ky2)
		//rotate
	print "---------------------"
	print kx_1,kx_2,ky_1,ky_2
	
	kx1= min( xrot(kx_1, ky_1, -ang),  xrot(kx_1, ky_2, -ang) )
	kx2=max( xrot(kx_2, ky_1, -ang),  xrot(kx_2, ky_2, -ang) )
	ky1= min( yrot(kx_1, ky_1, -ang),  yrot(kx_2, ky_1, -ang) )
	ky2=max( yrot(kx_1, ky_2, -ang),  yrot(kx_2, ky_2, -ang) )
		//round to simpler numbers (2 decimeal)
	print "----------------------------------------------------"
	print kx1,kx2,ky1,ky2
	 
	kx1=1E-2*round(kx1*1E2)
	kx2=1E-2*round(kx2*1E2)
	ky1=1E-2*round(ky1*1E2)
	ky2=1E-2*round(ky2*1E2)
	kz1=1E-2*round(kz1*1E2)
	kz2=1E-2*round(kz2*1E2)
		
	nkx=round((kx2-kx1)/kxinc + 1)
	nky=round((ky2-ky1)/kyinc + 1)
	nkZ=round((kZ2-kZ1)/kZinc + 1)

	variable r2d=180/pi
		
	print "kx:", kx1, kx2, kxinc, nkx, "   ky:", ky1, ky2, kyinc, nky, "kz",kz1,kz2,kzinc
	print r2d*asin( xrot(kx1, ky1, ang)/sqrt(xrot(kx1,ky1,ang)^2+kz1^2-V0))
	print r2d*asin( xrot(kx1, ky1, ang)/sqrt(xrot(kx1,ky1,ang)^2+kz2^2-V0))
	print r2d*asin( xrot(kx1, ky2, ang)/sqrt(xrot(kx1,ky2,ang)^2+kz1^2-V0))
	print r2d*asin( xrot(kx2, ky1, ang)/sqrt(xrot(kx2,ky1,ang)^2+kz1^2-V0))
	print r2d*asin( xrot(kx2, ky2, ang)/sqrt(xrot(kx2,ky2,ang)^2+kz1^2-V0))
	print r2d*asin( xrot(kx1, ky2, ang)/sqrt(xrot(kx1,ky2,ang)^2+kz2^2-V0))
	print r2d*asin( xrot(kx2, ky1, ang)/sqrt(xrot(kx2,ky1,ang)^2+kz2^2-V0))
	print r2d*asin( xrot(kx2, ky2, ang)/sqrt(xrot(kx2,ky2,ang)^2+kz2^2-V0))
	
	print"---------------"
	print 3.81*(xrot(kx1,ky1,ang)^2+yrot(kx1,ky1,ang)^2+kz1^2)-V0+WF
	print 3.81*(xrot(kx2,ky2,ang)^2+yrot(kx2,ky2,ang)^2+kz1^2)-V0+WF
	print 3.81*(xrot(kx1,ky2,ang)^2+yrot(kx1,ky2,ang)^2+kz1^2)-V0+WF
	print 3.81*(xrot(kx2,ky1,ang)^2+yrot(kx2,ky1,ang)^2+kz1^2)-V0+WF
// Create new image & polar/elev arrays	; set Kx-Ky regular grid scaling
	Make/O/N=(nkx, nky,nkz) volkxymap, root:tmp:polw, root:tmp:tiltw
	Setscale/I x kx1, kx2, "" volkxymap,  root:tmp:polw, root:tmp:tiltw
	Setscale/I y ky1, ky2, "" volkxymap,  root:tmp:polw, root:tmp:tiltw
	Setscale/I z,kz1,kz2,volkxymap
	 
	//print "Created: ", outn, "(", nkx, "x", nky, ") Ækxy=", kxinc, kyinc

// Calculate irregular grid polar& elev values
	// worry about misorientation correction

	//wave polw=root:tmp:polw
	//wave tiltw=root:tmp:tiltw
	//make/o/n=(zsize) kout=0.5123*sqrt(zoffset+(p-1)*zdelta)
	//root:tmp:polw=r2d*asin(x/sqrt(kout^2-y^2) )
	//root:tmp:tiltw=r2d*asin(y/kout)
	//polw=r2d*asin(xrot(x, y, ang)/sqrt(kout[r]^2-yrot(x, y, ang)^2) )
	//tiltw=r2d*asin(yrot(x, y, ang)/kout[r])
	//$outn=interp2D( $img, root:tmp:theta , y)
	//	all in one command without intermediate array
	//make/o/n=(xsize,ysize) root:tmp:tempimage
	//make/o/n=(nkx,nky) root:tmp:tempkimage
	//wave tempimage=root:tmp:tempimage
	//wave tempkimage=root:tmp:tempkimage
	//Setscale/I x kx1, kx2, tempkimage
	//Setscale/I y ky1, ky2, tempkimage
	//Setscale/p x xoffset, xdelta, tempimage
	//Setscale/p y yoffset, ydelta, tempimage
	
	volkxymap=interp3D( vol, r2d*asin( xrot(x, y, ang)/sqrt(xrot(x,y,ang)^2+z^2-V0) ) , r2d*asin(yrot(x, y, ang)/sqrt(xrot(x,y,ang)^2+yrot(x,y,ang)^2+z^2-V0)),3.81*(xrot(x,y,ang)^2+yrot(x,y,ang)^2+z^2)-V0+WF)
	

End


Proc KzPolar2Kxz( img, outn, BE, WF, xrng, yrng, kinc, thoff, Vo,dopt ) : GraphMarquee
//---------------
	string img, outn=StrVarOrDefault( "root:tmp:outkxz", "kxz")
	variable BE=NumVarOrDefault( "root:tmp:ben", 0), WF=NumVarOrDefault( "root:tmp:Wfct", 4.5)
	string kinc=StrVarOrDefault( "root:tmp:kxzincr", ".05,.05")
	//variable kxopt=NumVarOrDefault( "root:tmp:numkx", .1), kzopt=NumVarOrDefault( "root:tmp:numkz", .1)
	variable Vo=NumVarOrDefault( "root:tmp:Vinner", 10)	//, dopt=1 negative for varying Vo
	variable thoff=NumVarOrDefault( "root:tmp:poloff", 0)
	string xrng=ImgRange(0,""), yrng=ImgRange(1,"")
	variable dopt=1
	prompt img, "Kz(hv) vs Polar Image to convert", popup, WaveList("!*ct", ";", "WIN:,DIMS:2")
	prompt outn, "Output image name or extension (kxz)"
	prompt BE "Binding energy (eV)"
	prompt WF "Work Function (eV)"
	prompt xrng, "Polar range (min,max,inc)"
	prompt yrng, "hv range (min,max,inc)"
	prompt kinc, "k-increments:  kx,kz"
		//prompt kxopt, "kx inc (>0) or Rel. # of Kx points (<0)"
		//prompt kzopt, "kz inc (>0) or Rel. # of Kz points (<0)"
	prompt Vo, "Inner Potential, Vo,negative for varying Vo = 0"   
	prompt thoff, "Polar offset (deg)"
	prompt dopt, "Display new plot", popup, "Yes;No"
	
	
	//variable Vo_base=NumVarOrDefault( "root:tmp:Vo_base", 0)	//, dopt=1      //FW 0507/05
	//variable Vo_slope=NumVarOrDefault( "root:tmp:Vo_slope", 0)	//, dopt=1       //FW 0507/05
	//variable hv_base=NumVarOrDefault( "root:tmp:hv_base", 0)	//, dopt=1       //FW 0507/05
	//prompt Vo_base, "or  Vo="									  //FW 0507/05
	//prompt Vo_slope, "+                     *(hv-"					  //FW 0507/05
	//prompt hv_base, ")"					  //FW 0507/05
	
	NewDataFolder/O root:tmp
	string/G root:tmp:outkxz=outn, root:tmp:kxzincr=kinc
	variable/G root:tmp:ben=BE, root:tmp:Wfct=WF, root:tmp:Vinner=Vo
	//variable/G root:tmp:numkx=kxopt, root:tmp:numkz=kzopt
	
	GetVo(Vo)
	variable Vo_base=root:tmp:Vo_base	//, dopt=1      //FW 0507/05
	variable Vo_slope=root:tmp:Vo_slope	//, dopt=1       //FW 0507/05
	variable hv_base=root:tmp:hv_base	//, dopt=1       //FW 0507/05
	
	PauseUpdate;Silent 1
	string winnam=WinName(0,1)
	
//Check for marquee; subset selection of image
// obtained from ImgRange() function
	variable x1, x2, xinc, nx, y1, y2, yinc, ny
	if (strlen(xrng)>0)
		x1=ValFromList(xrng, 0,",")
		x2=ValFromList(xrng, 1,",")
		xinc=ValFromList(xrng, 2,",")
	else					//use full range of image
		xinc=DimDelta($img, 0)
		x1=DimOffset($img, 0)	
		x2=x1+(DimSize($img, 0)-1)*xinc
	endif
	nx=round( (x2-x1)/xinc + 1)
		
	if (strlen(xrng)>0)
		y1=ValFromList(yrng, 0,",")
		y2=ValFromList(yrng, 1,",")
		yinc=ValFromList(yrng, 2,",")
	else					//use full range of image
		yinc=DimDelta($img, 1)
		y1=DimOffset($img, 1)	
		y2=y1+(DimSize($img, 1)-1)*yinc
	endif
	ny=round( (y2-y1)/yinc + 1)
	
	
	Duplicate/O/R=(x1,x2)(y1,y2) $img, root:tmp:im
	////print x1, x2, xinc, nx
	////print y1, y2, yinc, ny
	
// Calculate kx range
	//Convert y1,y2 values to hv with Vo=WF=0
	variable hvy1, hvy2
		hvy1=sqrt(y1/0.523)^2
		hvy1=sqrt(y2/0.523)^2
//	variable hvmax=max(y1, y2), hvmin=min(y1,y2)
		variable hvmax=max(hvy1, hvy2), hvmin=min(hvy1,hvy2)
		print y1, y2, hvmin, hvmax
	variable Vomin=Vo_base+Vo_slope*(hvmin-hv_base)        //FW 0507/05
	variable Vomax=Vo_base+Vo_slope*(hvmax-hv_base)        //FW 0507/05
	
	variable kx1, kx2, kxinc, nkx
	kx1=min( sph2kx_(hvmax+BE-WF, x1-thoff,  0), sph2kx_(hvmin+BE-WF, x1-thoff,  0) )
	kx2=max( sph2kx_(hvmax+BE-WF, x2-thoff,  0), sph2kx_(hvmin+BE-WF, x2-thoff,  0) )
	kx1=1E-2*round(kx1*1E2)
	kx2=1E-2*round(kx2*1E2)	
	
// Calculate kz range
	variable kz1, kz2, kzinc, nkz, thmax=max(abs(x1-thoff), abs(x2-thoff))
	kz1=sph2kz_(hvmin+BE-WF, thmax,  0, Vomin)		//lower edge			//FW 0507/05
	//kz1=sph2kz_(hvmin+BE-WF, thmax,  0, Vo)		//lower edge
	kz2=sph2kz_(hvmax+BE-WF, 0,  0, Vomax)			// max at normal emission		//FW 0507/05
	//kz2=sph2kz_(hvmax+BE-WF, 0,  0, Vo)			// max at normal emission
	kz1=1E-2*round(kz1*1E2)
	kz2=1E-2*round(kz2*1E2)
	
	kxinc=ValFromList( kinc, 0, ",")
	kzinc=ValFromList( kinc, 1, ",")
	
//Manually fix
	kx1=-2
	kx2=2
	kz1=y1
	kz2=y2
	
	nkx=round((kx2-kx1)/kxinc + 1)
	nkz=round((kz2-kz1)/kzinc + 1)
	

	print "kx:", kx1, kx2, kxinc, nkx, "   kz:", kz1, kz2, kzinc, nkz

// Create new image & theta array	; set BE-Kpar regular grid scaling
	if (cmpstr(outn[0], "k")==0)
		outn=img+outn
	endif	
	Make/O/N=(nkx, nkz) $outn, root:tmp:theta, root:tmp:hv
	Setscale/I x kx1, kx2, "" $outn,  root:tmp:theta, root:tmp:hv
	Setscale/I y kz1, kz2, "" $outn,  root:tmp:theta, root:tmp:hv
	print "Created: ", outn, "(", nkx, "x", nkz, ") Ækxz=", kxinc, kzinc

// Calculate irregular grid theta values	

	//root:tmp:theta=Kp2Polar_( 3.81*(x^2+y^2)-Vo+BE,  x)+thoff			//KE=hv+BE-WF,  hv
	//root:tmp:hv=3.81*(x^2+y^2)-Vo-BE+WF								//=(KEin-Vo)-BE+WF
	//root:tmp:hv=(3.81*(x^2+y^2)-Vo_base+V0_slope*hv_base-BE+WF)/(1+V0_slope)						//FW 0507/05
	//root:tmp:theta=Kp2Polar_((3.81*(x^2+y^2)-Vo_base+V0_slope*hv_base-BE+WF)/(1+V0_slope)+BE-WF, x)+thoff			//KE=hv+BE-WF				//FW 0507/05
		
	//For linear hv dependence of Vo, the kx-kz convertion to hv-polar is solvable. For Vo=a+b*(hv-hv_0),
	// hv = (KEin-a+b*hv_0-BE+WF)/(1+b)										
	
	//$outn=interp2D( $img, root:tmp:theta , y)
	//	all in one command without intermediate array
//	$outn=interp2D( $img, Kp2Polar_((3.81*(x^2+y^2)-Vo_base+Vo_slope*hv_base-BE+WF)/(1+Vo_slope)+BE-WF,  x)+thoff , (3.81*(x^2+y^2)-Vo_base+Vo_slope*hv_base-BE+WF)/(1+Vo_slope) )  //FW 0507/05
	$outn=interp2D( $img, Kp2Polar_((3.81*(x^2+y^2)-Vo_base+Vo_slope*hv_base-BE+WF)/(1+Vo_slope)+BE-WF,  x)+thoff , y )  //FW 0507/05
	//$outn=interp2D($img,Kp2Polar_( 3.81*(x^2+y^2)-Vo+BE,  x)+thoff,3.81*(x^2+y^2)-Vo-BE+WF)
	//$outn=interp2D( $img, Kp2Polar_( K2KE_(x, y)-Vo+BE,  x)+thoff , K2KE_(x, y)-Vo-BE+WF )
	
	
	DoWindow/F $("Kxz_"+winnam)
	if ((dopt==1)*(V_flag==0))
		// get color table info from top graph(?)
		string CTstr=ColorTableStr("", img)		//may not work if not proper top graph
		display; appendimage $outn
		
		Label bottom "Kx (1/)"
		Label left "Kz (1/)"
//		ModifyGraph height={Plan,1,left,bottom}
		ModifyGraph width={Plan,1,bottom,left}
		ModifyGraph axThick=0.5,standoff=0
		
		string titlestr="\JC"+outn+": "+num2str(nkx)+" x "+num2str(nkz)
		titlestr+="\rVo="+num2str(Vomin)+"-"+num2str(Vomax)+", WF="+num2str(WF)+",  NE="+num2str(thoff)
		Textbox/N=title/F=0/A=MT/E titlestr
		if (strlen(CTstr)>0) 
			execute "ModifyImage "+outn+" "+CTstr
		else
			string ctable=img+"_CT"
			if (exists(ctable)==0)
				//ctable="YellowHot"
				ModifyImage $outn ctab= {*,*,YellowHot,0}
			else
				ModifyImage $outn cindex=$ctable
			endif
		endif
		DoWindow/C $("Kxz_"+winnam)
	endif

End

// ** Low-level functions from Tilt2K.ipf and SPh2K.ipf
// cannot define as "static" because used by Procs  (--> rewrite as functions)

Function K2KE_(kp, kz )			// KEin=3.81*(kp^2+kz^2) = KEout+Vo
//==============
	variable kp, kz
	return 3.81*(kp^2+kz^2) 	//- Vo		//KEout=3.81*kin^2-Vo
End

Function Kp2Polar_( KEout, Kp )		// alternate  (kout, kp)
//==============
	variable KEout, Kp
	//variable Kp=sqrt(kx^2+ky^2)
	variable r2d=180/pi	, Kout=0.5123*sqrt(KEout)
	return r2d*asin( Kp/Kout )
End


Function Sph2Kx_( KE, TH, PHI )
//==============
	variable KE, TH, PHI
	variable d2r=pi/180
	return 0.5123*sqrt(KE)*sin(d2r*TH)*cos(d2r*PHI)
End

Function Sph2Kz_( KE, TH, PHI, Vo )
//==============
	variable KE, TH, PHI, Vo
	variable  d2r=pi/180
	return  0.5123*sqrt( KE*(cos(d2r*TH))^2+Vo )
End

Function DualA2Kx_( KE, POL, ELV, geometry )
//==============
	variable KE, POL, ELV, geometry
	variable d2r=pi/180
	variable val
	if (geometry==1)				// Detector Polar-Elevation
		val = 0.5123*sqrt(KE)*sin(d2r*POL)*cos(d2r*ELV)
	else						// Detector Polar-Tilt
		val = 0.5123*sqrt(KE)*sin(d2r*POL)
	endif
	return val
End

Function DualA2Ky_( KE, POL, ELV, geometry )
//==============
	variable KE, POL, ELV, geometry	
	variable d2r=pi/180
	variable val
	if (geometry==1)				// Detector Polar-Elevation
		val = 0.5123*sqrt(KE)*sin(d2r*ELV)
	else						// Detector Polar-Tilt
		val = 0.5123*sqrt(KE)*cos(d2r*POL)*sin(d2r*ELV)
	endif
	return 0.5123*sqrt(KE)*sin(d2r*ELV)
End


Proc Kxyz_Kxyhv(voln, kvoln, V0,WF,kxinc,hvinc)
//-------------
	string voln, kvoln="Kxyhv"
	variable V0=10,WF=4.5,kxinc=0.02,hvinc=1
	prompt voln, "Kx-Ky-Kz volume to convertt", popup, WaveList("!*ct", ";", "DIMS:3")
	prompt kvoln, "Output volume name or extension (Kxyhv)"
	prompt V0 "Inner Potential (eV)"
	prompt WF "Work Function (eV)"
	prompt kxinc, "kx,ky-increment (0=same)"
	prompt hvinc, "hv-increment (0=auto)"
	
//	string voln=NameOfWave(vol)
	
	string inopt="/V0="+num2str(V0)+"/WF="+num2str(WF)
	string outopt="/KINC="+num2str(kxinc)+"/hv="+num2str(hvinc)
	if (stringmatch(kvoln,"Kxyhv"))
		//kvoln=voln+"_Kxyhv"
		//do not pass output destination /D option, will auto name as voln+"_Kxyhv"
	else
		outopt+="/D="+kvoln
	endif

	string cmd="Kxyz2Kxyhv("+voln+",\""+inopt+"\",\""+outopt+"\")"
	print cmd
	execute cmd
End

Function Kxyz2Kxyhv(vol,inopt,outopt)
// ==================
//  Input volume:  (Kx-Kx-Kz), 
//       Intensity = band energy (e.g. theory)
//   Conversion parameters:
//       V0 = inner potential
//       WF = work function
//  Output info:
//      KINC = kxinc  (default = same as input volume)
//      HV = hvinc  (default = auto-determine from Kz range & steps)
//      newvoln  --> Kxyhv  (Kx -Ky - hv)
	wave vol
	string inopt, outopt
	
	//energy correction:
	variable V0=10, WF=4.5
		V0=KeyVal("V0",inopt)
		V0=SelectNumber(numtype(V0)==2, V0, 10)		// use default if NAN
		WF=KeyVal("WF",inopt)
		WF=SelectNumber(numtype(WF)==2, WF, 4.5)		// use default if NAN
	
	// output k-increments
	variable kxinc=0.02, hvinc=1
	if (KeySet("KINC", outopt))
		kxinc=str2num(StringFromList(0, KeyStr("KINC", ",")) )
		kxinc=SelectNumber(numtype(kxinc)==2, kxinc, 0.02)		// use default if NAN
	endif
	if (KeySet("hv", outopt))

		hvinc=str2num(StringFromList(1, KeyStr("hv", ",")) )
		hvinc=SelectNumber(numtype(hvinc)==2, hvinc, 0)    //0=>auto-calculate
	endif
	
	string voln=NameOfWave(vol), newvoln
	newvoln=KeyStr("D",outopt)
	if ((strlen(newvoln)==0) + stringmatch(newvoln,"Kxzhv"))
		if (stringmatch(voln, "*Kxyz")==1)
			newvoln=voln[0,strlen(voln)-2]+"hv"		//replace just last character "z"
		else
			newvoln=voln+"_Kxyhv"
		endif
	endif
	
	NewDataFolder/O root:tmp
	
	variable nx=dimsize(vol,0)
	variable ny=dimsize(vol,1)
	variable nz=dimsize(vol,2)
	
	variable kx0 = dimoffset(vol,0)
	variable kx1 =kx1+dimdelta(vol,0)*(nx-1)
	variable ky0 = dimoffset(vol,1)
	variable ky1 =ky0+dimdelta(vol,1)*(ny-1)
	variable kz0=dimoffset(vol,2)
	variable kz1=kz0+dimdelta(vol,2)*(nz-1)

// Calculate hv range from normal emission
      variable hvmin = K2KE_(0, min(kz0,kz1) )-V0+WF		// assume BE=0, EF
	variable hvmax = K2KE_(0, max(kz0,kz1) )-V0+WF	
	if (hvinc==0)
		hvinc=(hvmax-hvmin)/(nz-1)
	endif
	if (hvinc<1)
		hvinc=round(hvinc*10)*0.1 
		hvmin=round(hvmin*10)*0.1
		hvmax=round(hvmax*10)*0.1
	else
		hvinc=round(hvinc )
		hvmin=round(hvmin)
		hvmax=round(hvmax)
	endif
	variable nhv = round((hvmax-hvmin)/hvinc + 1)
	print hvmin, hvmax, hvinc
	
	// Create new image & theta array	; set BE-Kpar regular grid scaling
	Make/O/N=(nx, ny,nhv) $newvoln//, root:tmp:theta, root:tmp:hv
	Setscale/I x kx0, kx1, $newvoln//,  root:tmp:theta, root:tmp:hv
	Setscale/I y ky0, ky1, $newvoln//,  root:tmp:theta, root:tmp:hv
	Setscale/P z hvmin, hvinc, $newvoln 			//photon energy
	//print "Created: ", outn, "(", nkx, "x", nkz, ") Ækxz=", kxinc, kzinc
	wave newvol = $newvoln

// Calculate with hv, V0, WF converson to kz	
//	newvol=interp3D( vol, kx, ky, kz=fct(hv, V0, WF) )
	newvol=interp3D( vol, x, y, sqrt( Sph2Kz_( z-WF, 0, 0, V0 )^2-x^2-y^2) )		// kz^2 = kin^2-kp^2,  kin
	
End