#pragma rtGlobals=1		// Use modern global access method.
#pragma version=1.2		// removed all of the old stuff
	Strconstant micron ="µm"			// use for Mac, I think that this works for Windows too!
//	Strconstant micron ="micron"		// use for Win
//
//	mu is hex B5 = 181
//	deg is hex A1 = 161
//	Windows deg is hex C1 = 193


// SliceType = 0 gives plane perp to axis 0, Igor calls this x, and I call it X (so graph will be H,F) (or H, Z)
// SliceType = 1 gives plane perp to axis 1, Igor calls this y, and I call it H (so graph will be X,F) (or X,Z)
// SliceType = 2 gives plane perp to axis 2, Igor calls this z, and I call it F or Z (so graph will be X,H)


Menu "SpotTilts"
	"4 Plot XHZ slice of 3d array, new",Make_Graph_slice($"")
	Submenu "   3d surface"
		"Show 3d surface",Make_Graph_3dSurface()
		"autoscale z",ModifySurfer autoscale=1
		"fix z scale",Set3dSurfaceZscale()
	End
	"3 Interpolate 3d grid, new",Interpolate3dTiltsFromXHF($"",NaN)
	"2 Make  XHF_tilt,XHF_Htilt,XHF_Xtilt",Make_XHF_tilt_Arrays()
	"1 Read In Raw Data",readInTheRawDataFile()
End





//	„readInTheRawDataFile()
//	      read in data from  'Macintosh HD:Users:tischler:data:sector 34:Ni on Si:2DAreaDepthScan.txt'
//	      putting data into folder 'root:X2DAreaDepthScan:'
//	  read in 110 wirescans, formed an array of (78,21,21)
//	„Make_Graph_slice($"")
//	„Make_XHF_tilt_Arrays()
//	  X range is [0,20]
//	  H range is [-3.53468,56.5685]
//	  F range is [0,56.5685]
//	  did no clipping of points, since the level=0 could not be found by FindLevel
//	  clip points below 3, removes ~9.2% of the points
//	  clip points below 3, removes ~8.1% of the points
//	  filling the 3-d array took 00:00:30
//	„Interpolate3dTiltsFromXHF($"",NaN)
//	  Interpolating from wave 'root:X2DAreaDepthScan:XHF_tilt', with resolution = 2µm
//	  started at 7:23:59 PM,   should finish at 7:27:06 PM
//	  creating 3-d wave named 'tilt'
//	  X   0  20  11  
//	  H   -4  57  32  
//	  F   0  57  30  
//	  to triangulate the 3-d array took 00:01:28
//	  the whole interpolation process took 00:11:06
//	„|
//	„ModifyImage Slice_Wave ctab= {-10,10,Terrain,1}
//	„ModifyContour Slice_Wave autoLevels={-1,1,11}



Function Make_Graph_slice(sourceWave)
	Wave sourceWave
	if (!WaveExists(sourceWave))
		if (strlen(WaveList("*", ";","DIMS:3" ))<1)
			DoAlert 0, "There are no 3d waves in this data folder, use the 'Data Browser' to move to the right folder"
			Execute "CreateBrowser"
			return 1
		endif
		String sourceWaveName, list = WaveList("*", ";","DIMS:3" )
		sourceWaveName = StringFromList(ItemsInList(list)-1,list)
		Prompt sourceWaveName,"name of 3d wave to show",popup,list	// "Xccd3d"
		DoPrompt "3d wave",sourceWaveName
		if (V_flag)
			return 1
		endif
		Wave sourceWave = $sourceWaveName
	endif
	if (!WaveExists(sourceWave))
		DoAlert 0, "Unable to find 3d wave, '"+sourceWaveName+"'"
		return 1
	endif
	if (WaveDims(sourceWave)!=3)
		DoAlert 0, NameOfWave(sourcewave)+" is "+num2istr(WaveDims(sourceWave))+"-d, not 3-d"
		return 1
	endif

	if (strlen(WinList("Graph_slice", "", "WIN:1")))
		DoWindow /F Graph_slice
		return 1
	endif
	String fldr= GetDataFolder(1)			// this is the folder where everything is created and stored
	Make/N=(10,10)/O Slice_Wave
	String noteStr = note(sourceWave)
	Variable XHZtype=0, XHFtype=0
	XHZtype = (numtype(NumberByKey("Flo", noteStr,"="))>0)
	XHFtype = !XHZtype
	noteStr = ReplaceStringByKey("sourceWave",noteStr,GetWavesDataFolder(sourceWave,2),"=")
	noteStr = ReplaceStringByKey("axisNames",noteStr,SelectString(XHFtype,"X,H,Z","X,H,F"),"=")
	Note/K Slice_Wave
	Note Slice_Wave, noteStr
	Display /W=(528,50,1078,663)/K=1 
	DoWindow /C Graph_slice
	AppendMatrixContour Slice_Wave
	String str = fldr+"SliceValue"
	Variable/G $str							// ensure existance of SliceType
	NVAR SliceValue=$str
	str = fldr+"SliceType"
	Variable/G $str
	NVAR SliceType=$str
	str = fldr+"phiSlice"
	if (exists(str)!=2)
		Variable/G $str=0					// ensure existance of phiSlice
	endif
	str = fldr+"CutResolution"
	if (exists(str)!=2)
		Variable/G $str=1.0				// ensure existance of resolution
	endif
	ModifyContour Slice_Wave autoLevels={*,*,11}, rgbLines=(0,0,0)		// ModifyContour Slice_Wave autoLevels={0,0.001,11}, rgbLines=(0,0,0)
	ModifyContour Slice_Wave labelBkg=1
	AppendImage Slice_Wave
	ModifyImage Slice_Wave ctab= {*,*,RedWhiteBlue,0}					// red is negative, white=0, blue is postive
	if (XHZtype)
		ModifyGraph swapXY=1
	endif
	ModifyGraph margin(top)=35,mirror=1,minor=1
	ModifyGraph/Z lSize('Slice_Wave=0')=3

	Label left SelectString(XHZtype , "Depth (F, \\U)","Z (\\U)")
	Label bottom "X (\\U)"
	SetAxis/A/R left
//	SetAxis/A/R bottom
	SetAxis/A bottom
	Cursor/P/I A Slice_Wave 1,1
	ShowInfo
	TextBox/N=text0/F=0/S=3/A=RB/X=3.06/Y=2.06 "spot tilt (radians)\r"+GetWavesDataFolder(sourceWave,2)
	ColorScale/N=text1/A=RC/X=-4.75/Y=-24.74 image=Slice_Wave, heightPct=30

	Variable lo=ceil(NumberByKey("Hlo",noteStr,"=")), hi=floor(NumberByKey("Hhi",noteStr,"="))
	SetVariable setvarSlice,pos={1,2},size={70,18},proc=SetSliceProc,title="H",fSize=12
	SetVariable setvarSlice,limits={lo,hi,1},value= SliceValue
	SetVariable setvarPhi,pos={137,2},size={71,18},proc=SetPhiCutProc,title="phi”"
	SetVariable setvarPhi,fSize=12,limits={0,180,0},value=phiSlice
	SetVariable setvarResolution,pos={220,2},size={145,18},title="resolution (µm)",fSize=12
	SetVariable setvarResolution,limits={0,inf,0},value= CutResolution,proc=SetResolutionProc
	sprintf str, "PopupMenu XHpopup value=\"%s\"", SelectString(XHFtype,"X;H;Z","X;H;F;phi")
	Execute str										//	PopupMenu XHpopup value="X;H;F"
	PopupMenu XHpopup ,pos={81,1},size={45,20},proc=XHFPopMenuProc, mode=2
	XHFPopMenuProc("",0,"H")					// last of all force an update
End
Function Make_Graph_3dSurface()
	// Mac V. 3.0.0
	if (exists("CreateSurfer") !=4)		// Do nothing if the Surface Plotter XOP is not available.
		DoAlert 0, "Surface Plotter XOP must be installed"
		return 1
	endif
	if (exists("Slice_WaveAbs")!=1 || exists("Slice_WaveLog")!=1)
		DoAlert 0, "Slice_WaveAbs or  Slice_WaveLog does not exists, try Make_Graph_slice() first"
		return 1
	endif
	Execute "CreateSurfer/K=1"
	MoveWindow 7,53,675,571 		//	left, top, right, bottom
	Execute "ModifySurfer  FactoryDefaults, Update=0"
	Execute "ModifySurfer expandWindow=0"
	Execute "ModifySurfer srcWave=Slice_WaveAbs"
	Execute "ModifySurfer secondWave=Slice_WaveLog"
	Execute "ModifySurfer surfaceColorWave=Slice_WaveLog"
	Execute "ModifySurfer  srcType=1,plotType=3,  setControlView=3"
	Execute "ModifySurfer  theta=28.9,  phi=311.8,  zScale=1,  xStep=1,  yStep=1"
	Execute "ModifySurfer  frame=895,  drawFrame=0, drawBox=1"
	Execute "ModifySurfer  numXTicks=2,  numYTicks=2,  numZTicks=2, drawTicks=5"
	Execute "ModifySurfer  noTicksX=0,  noTicksY=0,  noTicksZ=0"
	Execute "ModifySurfer palette=RedWhiteBlue, reverseColorMap=1"
	Execute "ModifySurfer  xmin=NaN,  xmax=NaN,  ymin=NaN,  ymax=NaN"
	Execute "ModifySurfer  surfaceColorFill=1,  grids=1"
	Execute "ModifySurfer  numContourLevels=12,  contour=1027"
	Execute "ModifySurfer  marker=16,  markerSize=1,  markerColorType=4"
	Execute "ModifySurfer  scatterDepthCue=1, imageType=1, BPP=16"
	//	Execute "ModifySurfer srcWave=Slice_WaveAbs"
	Execute "ModifySurfer  Update=1"
	Execute "ModifySurfer fullUpdate"
End
Function Set3dSurfaceZscale()
	if (exists("GetSurfer") !=4)		// Do nothing if the Surface Plotter XOP is not available.
		DoAlert 0, "Surface Plotter XOP must be installed"
		return 1
	endif
	Execute "GetSurfer surferName"			// store name in S_SurferName
	SVAR S_SurferName=S_SurferName
	if (strlen(S_SurferName)<=1)			// no surface plot is up
		return 1
	endif

	WaveStats/Q/M=1 Slice_WaveAbs
	Execute "ModifySurfer zmin=0,zmax="+num2str(V_max)
	Execute "ModifySurfer autoscale=0"
End

// this updates the slice surface.  It is called each time you want a new level.  To change the kind of slice (x, y, or z) call XHFPopMenuProc
Function SetSliceProc(ctrlName,SliceValue,varStr,varName) : SetVariableControl
	String ctrlName
	Variable SliceValue
	String varStr
	String varName
	if (cmpstr("SliceValue",varName))
		Abort "SetSliceProc only callable from the button"
	endif

	String fldr = GetWavesDataFolder(ImageNameToWaveRef("", "Slice_Wave" ),1)
	Wave Slice_Wave=$(fldr+"Slice_Wave")
	String noteStr = note(Slice_Wave)

	NVAR SliceType = $(fldr+"SliceType")		// 1 is the ususal (H)
	Wave source3d=$StringByKey("sourceWave", noteStr,"=")

	String fldrSav0= GetDataFolder(1)			// needed, Slice3d re-creates Slice_Wave in current folder
	SetDataFolder $fldr

	Make/N=3/O surfNormal_
	switch(SliceType)
		case 0:
			surfNormal_={1,0,0}				// for constant X plane
			break
		case 1:
			surfNormal_={0,1,0}				// for constant H (or Y) plane
			break
		case 2:
			surfNormal_={0,0,1}				// for constant F (or Z) plane
			break
		case 3:
			Variable phi = NumVarOrDefault(fldr+"phiSlice",0)*PI/180
			surfNormal_[0] = cos(phi)			// for a phi cut around F
			surfNormal_[1] = sin(phi)
			surfNormal_[2] = 0
			break
		default:
			Abort "cannot deal with a SliceType of "+num2str(SliceType)
	endswitch

	if (SliceType==3)							// special for a cut at arbitrary phi
		String SliceName
		Variable CutResolution = NumVarOrDefault(fldr+"CutResolution",1.0)
		SetVariable setvarPhi disable=0
		SetVariable setvarResolution disable=0
		if (numtype(CutResolution) || CutResolution<=0)
			DoAlert 0, "You cannot compute with a resolution ² 0"
			return 1
		endif
		SliceName = MakeArbitrarySlice(CutResolution,surfNormal_,SliceValue)	// recalculate a new cut, result goes into Slice_Wave
		SetDataFolder fldrSav0
		Wave Slice_Wave=$SliceName
		Label bottom "transverse (\\U)  at \\F'Symbol'f\\F]0="+num2str(phi*180/PI)+"”"
		if (!WaveExists(Slice_Wave))
			return 1
		endif
	else
		Slice3d(source3d,SliceType,SliceValue)		// recalculate a new cut
		Wave Slice_Wave=$(fldr+"Slice_Wave")
		SetDataFolder fldrSav0
		Variable xdim, ydim
		xdim = SelectNumber(SliceType-1,1,0,0)
		ydim = SelectNumber(SliceType-1,2,2,1)
		SetScale/P x DimOffset(source3d,xdim),DimDelta(source3d,xdim),WaveUnits(source3d,xdim), Slice_Wave
		SetScale/P y DimOffset(source3d,ydim),DimDelta(source3d,ydim),WaveUnits(source3d,ydim), Slice_Wave
		noteStr = ReplaceNumberByKey("SliceValue",noteStr,SliceValue,"=")
		Note/K Slice_Wave
		Note Slice_Wave, noteStr
		SetVariable setvarPhi disable=1
		SetVariable setvarResolution disable=1
	endif

	Duplicate/O Slice_Wave, Slice_WaveAbs, Slice_WaveLog				// these waves needed to make the 3d surface plot
	Slice_WaveAbs = abs(Slice_Wave)
	Slice_WaveLog = log(abs(Slice_Wave)) * sign(Slice_Wave)
	Slice_WaveLog = numtype(Slice_WaveLog) ? 0 : Slice_WaveLog

	String str = StringByKey("RECREATION",ImageInfo("","Slice_Wave",0))		//ctab= {*,*,RedWhiteBlue,0}
	str = StringByKey("ctab",str,"=")
	str = str[2,strlen(str)-2]
	if (strlen(StringFromList(0,str,",")+StringFromList(1,str,","))>27 || strsearch(str,"*,",0)>=0)	// probably auto-scaled
		WaveStats/Q/M=1 Slice_Wave				// these 3 lines make red negative, blue positive
		Variable colRange = max(abs(V_min),abs(V_max))
		colRange = colRange==0 ? sqrt(2)*1e-4 : colRange
		ModifyImage Slice_Wave ctab= {-colRange,colRange,RedWhiteBlue,0}	// red is negative, white=0, blue is postive
	endif
	KillWaves/Z surfNormal_
	if (cmpstr("setvarH",ctrlName))
		return 0
	endif
End

Function XHFPopMenuProc(ctrlName,popNum,popStr) : PopupMenuControl
	String ctrlName
	Variable popNum
	String popStr

	Wave sliceWave = ImageNameToWaveRef("", "Slice_Wave" )
	String noteStr = note(sliceWave)
	String fldr = GetWavesDataFolder(sliceWave,1)
	String str = fldr+"SliceType"
	Variable/G $str
	NVAR SliceType=$str
	SliceType = NumberByKey(popStr,"X:0;H:1;F:2;Z:2;phi:3")		// X=0, H=1, F or Z is 2, phi=3
	if (numtype(SliceType) || SliceType<0 || SliceType>3)
		Abort "in XHFPopMenuProc(), found SliceType="+num2str(SliceType)
	endif

	String xName, yName, axisNames = StringByKey("axisNames",noteStr,"=")+",phi"
	Variable ix,iy
	ix = str2num(StringFromList(SliceType,"1;0;0;3"))
	iy = str2num(StringFromList(SliceType,"2;2;1;2"))
//	ix = SelectNumber(SliceType-1,1,0,0)
//	iy = SelectNumber(SliceType-1,2,2,1)
	xName = StringFromList(ix,axisNames,",")
	yName = StringFromList(iy,axisNames,",")

	String lables = "X:X (\\U);H:H (\\U);Z:Z (\\U);F:Depth (F, \\U);phi:at \\F'Symbol'f\\F]0=x” (\\U)"
	Label bottom StringByKey(xName, lables)
	Label left StringByKey(yName,lables)

	Variable lo=ceil(NumberByKey(popStr+"lo",noteStr,"=")), hi=floor(NumberByKey(popStr+"hi",noteStr,"="))
	SetVariable setvarSlice,limits={lo,hi,1}
	Variable SliceValue=NumVarOrDefault(fldr+"SliceValue",(lo+hi)/2)
//	SetVariable setvarSlice,title=popStr
	SetVariable setvarSlice,title=SelectString(stringmatch(popStr,"phi"),popStr,"r")
	SetSliceProc("",SliceValue,"","SliceValue")		// force a re-calc of the surface wave
End
Function SetPhiCutProc(ctrlName,varNum,varStr,varName) : SetVariableControl
	String ctrlName
	Variable varNum
	String varStr
	String varName
	XHFPopMenuProc("",-1,"phi")
End
Function SetResolutionProc(ctrlName,varNum,varStr,varName) : SetVariableControl
	String ctrlName
	Variable varNum
	String varStr
	String varName

	if (varNum<=0)
		DoAlert 0, "you set the resolution to ² 0, what were you thinking?"
	elseif (abs(varNum)<0.5)
		DoAlert 0, "Caution, this fine a resolution will take a long (really long) time to interpolate!"
	endif
End




Function/S 	MakeArbitrarySlice(resolution,normal,SliceValue)	// recalculate a new cut, result goes into Slice_Wave, comes from M_3DVertexList
	Variable resolution			// resolution, grid spacing of output wave
	Wave normal				// surface normal
	Variable SliceValue			// intercept displaced along normal

	Wave M_3DVertexList=M_3DVertexList
//	Wave source3d = root:X2DAreaDepthScanPeakP_16E728:Htilt
//	if (!WaveExists(M_3DVertexList) || !WaveExists(source3d))
	if (!WaveExists(M_3DVertexList))
//		DoAlert 0, "Error in MakeArbitrarySlice(), M_3DVertexList or source3d does not exist"
		DoAlert 0, "Error in MakeArbitrarySlice(), M_3DVertexList does not exist"
		return ""
	endif
	if (normal[2]!=0)
		Abort "MakeArbitrarySlice(), normal[2]!=0"
	endif

	Wave Slice_Wave=Slice_Wave
	String noteStr = note(Slice_Wave)
	noteStr = ReplaceNumberByKey("SliceValue",noteStr,SliceValue,"=")

	Variable xlo,xhi, hlo,hhi, flo, fhi,  nt,nz
	xlo = floor(NumberByKey("Xlo",noteStr,"="))		// extent of box in x,y
	xhi = ceil(NumberByKey("Xhi",noteStr,"="))
	hlo = floor(NumberByKey("Hlo",noteStr,"="))
	hhi = ceil(NumberByKey("Hhi",noteStr,"="))
	flo = floor(NumberByKey("Flo",noteStr,"="))
	fhi = ceil(NumberByKey("Fhi",noteStr,"="))
	nz = round((fhi-flo)/resolution + 1)

	Variable x0=SliceValue*normal[0], y0=SliceValue*normal[1]
	Variable perpX=-normal[1], perpY=normal[0],len
	len = sqrt(perpX^2 + perpY^2)
	perpX /= len										// the line along the cut is y=y0+t*perpY, x=x0+t*perpX, t is the parametric parameter
	perpY /= len

	Make/N=4/O/D t_temp_							// t where the line in x,y intersects the 4 lines making up the rectangle
	t_temp_[0] = (xlo-x0)/perpX
	t_temp_[1] = (xhi-x0)/perpX
	t_temp_[2] = (hlo-y0)/perpY
	t_temp_[3] = (hhi-y0)/perpY
	Sort t_temp_,t_temp_								// order t1,t2,t3,t4, and pick the middle two for tlo and thi
	Variable tlo,thi										// width of the cut, where it hits box defined by (xlo,ylo), (xhi,yhi)
	tlo = t_temp_[1]
	thi = t_temp_[2]
	KillWaves/Z t_temp_

	// to check that this line crosses inside of box, check if midpoint is internal
	Variable tmid = (thi+tlo)/2, xmid,ymid
	xmid = x0+tmid*perpX
	ymid = y0+tmid*perpY
	if (xmid<=xlo || xmid>=xhi || ymid<=hlo || ymid>=hhi)
		DoAlert 0, "cut lies outside data"
		return ""
	endif

	nt = round((thi-tlo)/resolution + 1)
	Make/N=(nt,nz)/O Slice_Wave
	SetScale/I x tlo,thi,"µm", Slice_Wave
	SetScale/I y flo,fhi,"µm", Slice_Wave
	Note/K Slice_Wave
	Note Slice_Wave, noteStr

	Wave XHFwave=$StringByKey("XHFsource", note(M_3DVertexList),"=")
	if (WaveExists(XHFwave))
		Slice_Wave = Interp3d(XHFwave,x0+x*perpX,y0+x*perpY,y,M_3DVertexList)		// This line takes a long time!
	else
		DoAlert 0, "could not find XHFwave from M_3DVertexList in MakeArbitrarySlice()"
		Slice_Wave=0
	endif
	return GetWavesDataFolder(Slice_Wave,2)
End




Function Interpolate3dTiltsFromXHF(XHF_w,resolution)
	Wave XHF_w		// the XHF wave, probably either root:raw:XHF_tilt, root:raw:XHF_Htilt, or root:raw:XHF_Xtilt
	Variable resolution	// resolution in final array (µm)

	Variable defaultResolution = NumVarOrDefault("CutResolution",1.0)
	defaultResolution = defaultResolution<=0 ? NaN : defaultResolution
	resolution = (resolution<=0) ? NaN : resolution			// force 0 or negative to be invalid

	if (!WaveExists(XHF_w) || numtype(resolution) || resolution<=0)
		if (strlen(WaveList("XHF_*",";","DIMS:2"))<1 && resolution>0)
			DoAlert 0, "There are no 3d waves in this data folder, use the 'Data Browser' to move to the right folder"
			Execute "CreateBrowser"
			return 1
		endif
		String wList=WaveList("XHF_*",";","DIMS:2")
		String XHF_name="XHF_Xtilt"
		Prompt XHF_name,"XHF file to sort into 3d array",popup,wList
		Prompt resolution, "resolution in final array ("+micron+")"
		if (!WaveExists(XHF_w) && (numtype(resolution) ||resolution<=0))
			resolution = numtype(resolution) ? defaultResolution : max(abs(resolution),resolution)
			DoPrompt "pick XHF file and resolution", XHF_name,resolution
		elseif(!WaveExists(XHF_w))
			DoPrompt "pick XHF file", XHF_name
		else
			resolution = numtype(resolution) ? defaultResolution : max(abs(resolution),resolution)
			DoPrompt "set resolution", resolution
		endif
		if (V_flag)
			return 1
		endif
		if (!WaveExists(XHF_w))
			Wave XHF_w = $XHF_name
		endif
	endif
	if (!WaveExists(XHF_w))
		DoAlert 0, "The XHF_ wave cannot be found"
		return 1
	endif
	if (numtype(resolution) || resolution<=0)
		DoAlert 0, "a resolution of "+num2str(resolution)+" is not valid"
		return 1
	endif

	Variable resFactor = 1/(resolution^3)
	printf "Interpolating from wave '%s', with resolution = %g%s\r",GetWavesDataFolder(XHF_w,2),resolution,micron
	String str
	DoAlert 1,"This procedure could take ~"+num2istr(25*resFactor)+"min, do it?"
	if (V_flag>1)
		return 1
	endif
	Variable now = dateTime
	printf "started at %s,   should finish at %s\r",Secs2Time(now,1),Secs2Time(now+25*60*resFactor,1)
	DoUpdate

	Variable N = DimSize(XHF_w,0)
	String noteStr = note(XHF_w)
	String outName = NameOfWave(XHF_w)
	outName = outName[4,Inf]
	printf "creating 3-d wave named '%s'\r",outName

	Variable xlo,xhi, hlo,hhi, flo, fhi,  nx,ny,nz
	xlo = floor(NumberByKey("Xlo",noteStr,"="))
	xhi = ceil(NumberByKey("Xhi",noteStr,"="))
	hlo = floor(NumberByKey("Hlo",noteStr,"="))
	hhi = ceil(NumberByKey("Hhi",noteStr,"="))
	flo = floor(NumberByKey("Flo",noteStr,"="))
	fhi = ceil(NumberByKey("Fhi",noteStr,"="))
	nx = round((xhi-xlo)/resolution + 1)
	ny = round((hhi-hlo)/resolution + 1)
	nz = round((fhi-flo)/resolution + 1)
	//	print "X ",xlo,xhi,nx
	//	print "H ",hlo,hhi,ny
	//	print "F ",flo,fhi,nz

	Variable timer1=startMSTimer, timer2=startMSTimer
	Triangulate3d XHF_w
	Wave M_3DVertexList=M_3DVertexList
//	printf "to triangulate the 3-d array took %s\r",num2sexigesmal(stopMSTimer(timer1)*1e-6,0)
	printf "to triangulate the 3-d array took %s\r",Secs2Time(stopMSTimer(timer1)*1e-6,5,0)
	Note/K M_3DVertexList
	Note M_3DVertexList,"XHFsource="+GetWavesDataFolder(XHF_w,2)
	Interpolate3D /RNGX={xlo,resolution,nx}/RNGY={hlo,resolution,ny}/RNGZ={flo,resolution,nz}/DEST=$outName triangulationWave=M_3DVertexList, srcWave=XHF_w
	Wave tilt = $outName
	tilt = numtype(tilt[p][q][r]) ? 0 : tilt[p][q][r]
	SetScale/I x xlo,xhi,micron, tilt
	SetScale/I y hlo,hhi,micron, tilt
	SetScale/I z flo,fhi,micron, tilt
	Note/K tilt
	Note tilt, noteStr
	Variable/G CutResolution = resolution
	beep
//	printf "the whole interpolation process took %s\r",num2sexigesmal(stopMSTimer(timer2)*1e-6,0)
	printf "the whole interpolation process took %s\r",Secs2Time(stopMSTimer(timer2)*1e-6,5,0)
End





// This routine makes arrays of the tilts from the raw CCD data.  It makes three arrays, XHF_Xtilt, XHF_Htilt, XHF_tilt
// XHF_Xtilt has four columns and many rows.  Each row contains x,h,f,(tilt in x direction) for a measured spot.
// Likewise XHF_Htilt is similar, but for the tilt in the H direction, and XHF_tilt is for the total tilt.  This routine does no
// interpolating; that is left for another routine, probably testLong()
// Fixed this so the sign of H is correct.  The sample move in H which is like moving the beam in -H
Function Make_XHF_tilt_Arrays()
	Wave XX = XX					// these are the arrays with the information that I need
	Wave HH = HH
	Wave FF = FF
	Wave Xccd = Xccd
	Wave Yccd = Yccd
	Wave Intensity = Intensity
	if (!WaveExists(XX) || !WaveExists(HH) || !WaveExists(FF) || !WaveExists(Xccd) || !WaveExists(Yccd) || !WaveExists(Intensity))
		DoAlert 0, "The waves XX,HH,FF,Xccd,Yccd, & Intensity are not in this folder, use the 'Data Browser' to move to the right folder"
		Execute "CreateBrowser"
		return 1
	endif

	String noteStr = note(XX)
	WaveStats/Q XX
	noteStr = ReplaceNumberByKey("Xlo",noteStr,V_min,"=")
	noteStr = ReplaceNumberByKey("Xhi",noteStr,V_max,"=")
	printf "X range is [%g,%g]\r",V_min,V_max
	WaveStats/Q HH
	noteStr = ReplaceNumberByKey("Hlo",noteStr,V_min,"=")
	noteStr = ReplaceNumberByKey("Hhi",noteStr,V_max,"=")
	printf "H range is [%g,%g]\r",V_min,V_max
	WaveStats/Q FF
	noteStr = ReplaceNumberByKey("Flo",noteStr,V_min,"=")
	noteStr = ReplaceNumberByKey("Fhi",noteStr,V_max,"=")
	printf "F range is [%g,%g]\r",V_min,V_max

	Variable N=numpnts(XX)
	Make/N=(N,4)/O XHF_Xtilt, XHF_Htilt, XHF_tilt
	Note/K XHF_Xtilt
	Note/K XHF_Htilt
	Note/K XHF_tilt
	Note XHF_Xtilt noteStr
	Note XHF_Htilt noteStr
	Note XHF_tilt noteStr
	XHF_Xtilt[][0] = XX[p] ;	XHF_Xtilt[][1] = HH[p];	XHF_Xtilt[][2] = FF[p]
	XHF_Htilt[][0] = XX[p] ;	XHF_Htilt[][1] = HH[p] ;	XHF_Htilt[][2] = FF[p]
	XHF_tilt[][0] = XX[p] ;		XHF_tilt[][1] = HH[p] ;		XHF_tilt[][2] = FF[p]

	Variable tiltX,tiltH,totalTilt		// tilt of the spot from the reference (radians)
	Variable i, zdepth

	Variable timer = startMSTimer
	for (i=0;i<N;i+=1)
		zdepth = sqrt(2)*FF[i]
		spot2tilt(noteStr, Xccd[i],Yccd[i],zdepth,  tiltX,tiltH,totalTilt)	// set tiltX, tiltH,  and totalTilt from ccd position
		XHF_Xtilt[i][3] = tiltX
		XHF_Htilt[i][3] = tiltH
		XHF_tilt[i][3] = totalTilt
	endfor

	// This section removes the points intensity below V_LevelX
	Make/N=100/O Intensity_Hist
	Histogram/B=1 Intensity,Intensity_Hist
	WaveStats/Q Intensity_Hist
	FindLevel /Q/R=(-Inf,V_maxloc) Intensity_Hist, V_maxloc*.01	//	Find where Histogram drops to (V_maxloc*.01)
	Variable j
	if (V_flag==0)
		printf "clip points below %g, removes ~%.1f%% of the points\r",V_LevelX,100*sum(Intensity_Hist,0,V_LevelX)/sum(Intensity_Hist,0,Inf)
		for (i=0,j=0;i<N;i+=1)
			if (Intensity[i]<V_LevelX)
				DeletePoints/M=0 j,1, XHF_Xtilt,XHF_Htilt,XHF_tilt
			else
				j += 1
			endif
		endfor
		printf "clip points below %g, removes ~%.1f%% of the points\r",V_LevelX,100*(N-DimSize(XHF_Xtilt,0))/N
	else
		printf "did no clipping of points, since the level=%g could not be found by FindLevel\r",V_maxloc*.01
		DoAlert 1, "no clipping done, do you want to manually set a threshold on intensity?"
		if (V_flag==1)
			Variable threshold=3
			Prompt threshold, "intensity threshold"
			DoPrompt "threshold",threshold
			if (!V_flag || threshold>=0)
				printf "clip points below %g, removes ~%.1f%% of the points\r",threshold,100*sum(Intensity_Hist,0,threshold)/sum(Intensity_Hist,0,Inf)
				for (i=0,j=0;i<N;i+=1)
					if (Intensity[i]<threshold)
						DeletePoints/M=0 j,1, XHF_Xtilt,XHF_Htilt,XHF_tilt
					else
						j += 1
					endif
				endfor
				printf "clip points below %g, removes ~%.1f%% of the points\r",threshold,100*(N-DimSize(XHF_Xtilt,0))/N
			endif
		endif
	endif
	printf "filling the 3-d array took %s\r",num2sexigesmal(stopMSTimer(timer)/1e6,0)
	printf "filling the 3-d array took %s\r",Secs2Time(stopMSTimer(timer)/1e6,5,0)
	KillWaves/Z Intensity_Hist
End


//	This fuction converts the spot position on the CCD (Xccd,Yccd) to tilts perpendicular to the 
//	surface normal (assumed to be -F).  All tilts are returned as radians
//	it gets all of the other geometry informaiton from the values in list.  
//	It is only called by Make_XHF_tilt_Arrays()
Static Function spot2tilt(list,dXccd,dYccd,zdepth,tiltX,tiltH,totalTilt)
	String list
	Variable dXccd,dYccd				// distance on CCD from un-tilted spot (pixels)
	Variable zdepth						// depth in z-direction below surface for this spot
	Variable &tiltX,&tiltH,&totalTilt	// tilt of the spot from the reference (radians)
										// tilts are toward X=(100) and the -H=(0-1-1)
//	   CCDX0=1104.24;CCDY0=986.833;psizeX=23.9971;psizeY=24.073;
//	   x0=210.757;y0=59498.8;z0=-2094.17;
//	   ki=0.0, 0.0, 1.0;Zlo=3;Zhi=80;Xlo=-0;Xhi=20;Hlo=-0;Hhi=7.07000017166138;

	// make the waves I need to work with
	String wName
	wName = UniqueName("kincident",1,0)
	Make/N=3/D $wName
	Wave ki = $wName							// direction of incident beam
	wName = UniqueName("kfinal",1,0)
	Make/N=3/D $wName
	Wave kf = $wName							// direction of exit beam (to the spot)
	wName = UniqueName("qvec",1,0)
	Make/N=3/D $wName
	Wave q0 = $wName							// direction of the q vector for the un-tilted spot
	wName = UniqueName("qvec",1,0)
	Make/N=3/D $wName
	Wave qv = $wName							// direction of the q vector for the tilted spot

	// extract all variables from 'list'
	Variable x0=NumberByKey("x0",list,"=")	// (x,y,z) position of un-tilted spot (µm)
	Variable y0=NumberByKey("y0",list,"=")	// this xyz is in the beamline coord system
	Variable z0=NumberByKey("z0",list,"=")
	Variable psizeX=NumberByKey("psizeX",list,"=")	// size of pixel in x-direction on CCD (µm)
	Variable psizeY=NumberByKey("psizeY",list,"=")	// size of pixel in z-direction on CCD (µm)
	String kiStr = StringByKey("ki",list,"=")

	ki[0] = str2num(StringFromList(0,kiStr,","))
	ki[1] = str2num(StringFromList(1,kiStr,","))
	ki[2] = str2num(StringFromList(2,kiStr,","))
	normalize(ki)

	kf = {x0,y0,z0-zdepth}						// kf = {x0,y0,z0} in beamline coordinate system
	normalize(kf)
	q0 = kf-ki									// calculte q of the un-tilted spot
	normalize(q0)

	kf = {x0+(dXccd*psizeX),y0,z0-zdepth+(dYccd*psizeY)}	// remember Yccd is in the z direction
	normalize(kf)
	qv = kf-ki									// calculate q of the tilted spot
	normalize(qv)

	wName = UniqueName("qvec",1,0)			// make some more vectors, qx, qh, & dq
	Make/N=3/D $wName
	Wave qx = $wName
	wName = UniqueName("qvec",1,0)
	Make/N=3/D $wName
	Wave qh = $wName
	wName = UniqueName("dq",1,0)
	Make/N=3/D $wName
	Wave dq = $wName

	// resolve all tilts in terms of the outward-surface normal (the -F direction), or in terms of xyz it is (0,1,-1)
	// the two perpendicular directions are (1,0,0) and (0,-1,-1),    and (1,0,0) x (0,-1,-1) =  (0,1,-1)
	qx = {1,0,0}								// direction of the tilt of the q vector toward +Xhat direction
	qh = {0,-1,-1}								// direction of the tilt of the q vector toward -Hhat direction
	qh /= sqrt(2)
	dq = qv-q0									// total deltaQ beyween qv and q0

	tiltX = asin(MatrixDot(dq,qx))				// tilt of q towad qx direction, toward (1,0,0)
	tiltH = asin(MatrixDot(dq,qh))				// tilt of q towad qh direction, toward (0,-1,-1)
	totalTilt = acos(MatrixDot(qv,q0))			// total tilt of qv away from q0, always positive, (remember that |qv|==1)
//	tiltX = 2*asin(MatrixDot(dq,qx)/2)		// tilt of q towad qx direction, toward (1,0,0)
//	tiltH = 2*asin(MatrixDot(dq,qh)/2)		// tilt of q towad qh direction, toward (0,-1,-1)
//	totalTilt = 2*asin(norm(dq)/2)			// total tilt of q away from q0, always positive, (remember that |qv|==1)

	Killwaves/Z ki,kf,q0,qv,dq,qx,qh
	return 0
End







Function/T readInTheRawDataFile()
	String/G root:rawDataFullPathName
//	if (exists("root:rawDataFullPathName")!=2)
//	endif
	SVAR rawDataFullPathName=root:rawDataFullPathName
	Variable refNum

	// see if rawDataFullPathName exists?
	Open/R/Z=1 refNum rawDataFullPathName
	if (refNum>0)
		Close refNum
	endif

	if (!V_flag)								// V_flag==0 if file exists
		DoAlert 1, "use: "+rawDataFullPathName
		V_flag = (V_flag==2)
	endif
	if (V_flag)
		rawDataFullPathName = ""
	endif
	Open/R/M="Wirescan peak positions from Wenge"/P=home/Z=2 refNum rawDataFullPathName
	if (refNum<=0)
		Abort "unable to re-open the data file"
	endif
	rawDataFullPathName = S_fileName		// save file name for future use

	String line								// determine if this file is a new or old type
	FReadLine refNum, line
	FSetPos refNum,0
	Close refNum
	String noteStr=""
	if (char2num(line[0])==char2num("$"))
		noteStr = readInNewRawDataFile(rawDataFullPathName)	//	new
	else
		Abort "You have chosen an old style file, they are not supported"
	endif
	if (strlen(noteStr)<1)
		DoAlert 0, "Unable to properly read "+rawDataFullPathName
		return ""
	endif
	return noteStr
End



Function/T readInNewRawDataFile(FullFileName)
	String FullFileName

	String noteStr = readSmarterHeader(FullFileName)
	if (strlen(noteStr)<1)
		return ""
	endif

	String fldrSav0= GetDataFolder(1)
	Variable refNum
	Open /R/Z=1 refNum as FullFileName
	if (strlen(S_fileName)<1 || V_flag)
		return ""
	endif
	String pathSep
	if (stringmatch(IgorInfo(2),"Macintosh"))
		pathSep = ":"
	elseif (stringmatch(IgorInfo(2),"Windows"))
		pathSep = "\\"
	else
		Abort "This is neither Mac nor Win, what is this computer?"
	endif
	String fldrName = ParseFilePath(3,FullFileName, pathSep, 0, 0)
	fldrName = CleanupName(fldrName, 0)
	if (stringmatch(fldrName,"raw"))
		fldrName = UniqueName(fldrName,11,0)
	endif

	if (CheckName(fldrName,11))
		DoAlert 2, "This data folder already exists, delete the folder and re-load the data?"
		if (V_flag==1)
			KillDataFolder $fldrName
		elseif (V_flag==2)
			fldrName = UniqueName(fldrName,11,0)
		else
			Close refNum
			return ""
		endif
	endif
	NewDataFolder /O/S $fldrName

//	FSetPos refNum,0
//	FStatus refNum			// get V_filePos
	if (strlen(findTagInFile(refNum,"wirescan"))<1)			// found the tag
		Close refNum
		SetDataFolder fldrSav0
		return ""												// unable to find $wirescan
	endif

	LoadWave/N=wire/L={0,0,0,0,4}/G/Q FullFileName
	Variable NwireScans = floor(V_flag/4)						// number of wire scans

	Make/N=(NwireScans)/O Xs,Hs
	Xs = NaN
	Hs = NaN
	Variable i
	FSetPos refNum,0
	for (i=0;i<NwireScans;i+=1)
		if (strlen(findTagInFile(refNum,"wirescan"))<1)		// found the tag
			Close refNum
			SetDataFolder fldrSav0
			DoAlert 0, "Unable to find enough $wirescan tags"
			return ""											// unable to find $wirescan
		endif
		Xs[i] = str2num(findTagInFile(refNum,"X"))			// X of PM500 (mm)
		Hs[i] = str2num(findTagInFile(refNum,"H"))			// H of PM500 (mm)
	endfor
	Close refNum
	printf "    read in data from  '%s'\r",FullFileName
	printf "    putting data into folder '%s'\r",GetDataFolder(1)

	WaveStats/Q Xs
	if (V_numNans+V_numINFs)
		SetDataFolder fldrSav0
		DoAlert 0, "Unable to read enough valid $X tags"
		return ""
	endif
	WaveStats/Q Hs
	if (V_numNans+V_numINFs)
		SetDataFolder fldrSav0
		DoAlert 0, "Unable to read enough valid $H tags"
		return ""
	endif

	// rename waves to useful values
	for (i=0;i<NwireScans;i+=1)
		Rename $("wire"+num2istr(4*i+0)) $("z"+num2istr(i))
		Rename $("wire"+num2istr(4*i+1)) $("xccd"+num2istr(i))
		Rename $("wire"+num2istr(4*i+2)) $("yccd"+num2istr(i))
		Rename $("wire"+num2istr(4*i+3)) $("intens"+num2istr(i))
	endfor

	// determine dimensionality from Xs[] and Hs[]
	Variable threshold = 0.18e-3								// (mm) a motion of the PM500 greater than this is intentional (less is jitter)
	Variable Xoff,dX, dimX=1									// for setting scaling in X-dimension
	Variable Hoff,dH, dimH=1
	Variable Zoff,dZ, dimZ=1
	FindScalingFromVec(Xs,threshold,Xoff,dX,dimX)
	FindScalingFromVec(Hs,threshold,Hoff,dH,dimH)

	Variable N = 0												// total number of Z points
	for (i=0;i<NwireScans;i+=1)
		N += numpnts($("z"+num2istr(i)))
	endfor
	Make/N=(N) Zs
	Zs = NaN
	Variable j0=0,j1
	for (i=0;i<NwireScans;i+=1)
		Wave zwave = $("z"+num2istr(i))
		j1 = j0+numpnts(zwave)-1
		if (j1<j0)
			continue
		endif
		Zs[j0,j1] = zwave[p-j0]
		j0 += j1-j0+1
	endfor
	FindScalingFromVec(Zs,threshold,Zoff,dZ,dimZ)
	KillWaves/Z zs

	printf "read in %d wirescans, formed an array of (%d,%d,%d)\r",NwireScans,dimZ,dimX,dimX
	Make/N=(dimZ,dimX,dimH)/O Xccd3d,Yccd3d,Intensity3d
	SetScale/P x Zoff,dZ,micron,Xccd3d,Yccd3d,Intensity3d
	SetScale/P y Xoff,dX,micron,Xccd3d,Yccd3d,Intensity3d
	SetScale/P z Hoff,dH,micron,Xccd3d,Yccd3d,Intensity3d
	noteStr = ReplaceNumberByKey("Zlo",noteStr,Zoff,"=")
	noteStr = ReplaceNumberByKey("Zhi",noteStr,Zoff+(dimZ-1)*dZ,"=")
	noteStr = ReplaceNumberByKey("Xlo",noteStr,Xoff,"=")
	noteStr = ReplaceNumberByKey("Xhi",noteStr,Xoff+(dimX-1)*dX,"=")
	noteStr = ReplaceNumberByKey("Hlo",noteStr,Hoff,"=")
	noteStr = ReplaceNumberByKey("Hhi",noteStr,Hoff+(dimH-1)*dH,"=")
	Note/K Xccd3d  ;  Note Xccd3d, noteStr
	Note/K Yccd3d  ;  Note Yccd3d, noteStr
	Note/K Intensity3d  ;  Note Intensity3d, noteStr
	Xccd3d = NaN
	Yccd3d = NaN
	Intensity3d = NaN

	// now fill the arrays Xccd3d & Yccd3d
	Variable j,iz,ix,ih
	for (i=0;i<NwireScans;i+=1)
		Wave zwave = $("z"+num2istr(i))
		Wave ccdx = $("xccd"+num2istr(i))
		Wave ccdy = $("yccd"+num2istr(i))
		Wave intens = $("intens"+num2istr(i))

		// X is Xs[i], H is Hs[i], and Z comes from zwave
		ih = (Hs[i]-Hoff)/dH
		ix = (Xs[i]-Xoff)/dX
		N = numpnts(zwave)
		for (j=0;j<N;j+=1)
			iz = (zwave[j]-Zoff)/dZ
			Xccd3d[iz][ix][ih] = ccdx[j]
			Yccd3d[iz][ix][ih] = ccdy[j]
			Intensity3d[iz][ix][ih] = intens[j]
		endfor
		KillWaves/Z zwave,ccdx,ccdy,intens
	endfor

	Variable ccdX0,ccdY0
	ccdX0 = NumberByKey("CCDX0",noteStr,"=")
	ccdY0 = NumberByKey("CCDY0",noteStr,"=")
	if (numtype(ccdX0+ccdY0))
		SetDataFolder fldrSav0
		DoAlert 0, "cannot find valid CCDX0 or CCDY0"
		return ""
	endif
	Xccd3d -= ccdX0
	Yccd3d -= ccdY0

	// finally create the long 1-d arrays that can be used for interpolating
	WaveStats/Q Intensity3d
	Make/N=(V_npnts) XX,HH,FF,Xccd,Yccd,Intensity
	Variable xval,hval,fval
	i = 0
	for (ih=0;ih<dimH;ih+=1)
		hval = ih*dH+Hoff
		for (ix=0;ix<dimX;ix+=1)
			xval = ix*dX+Xoff
			for (iz=0;iz<dimZ;iz+=1)
				if (numtype(Intensity3d[iz][ix][ih])==0 && Intensity3d[iz][ix][ih]>0)
					Xccd[i] = Xccd3d[iz][ix][ih]		// save valid point
					Yccd[i] = Yccd3d[iz][ix][ih]
					Intensity[i] = Intensity3d[iz][ix][ih]
					fval = (iz*dZ+Zoff)/sqrt(2)
					XX[i] = xval
					HH[i] = -hval + fval
					FF[i] = fval
					i += 1
				endif
			endfor
		endfor
	endfor
	Note/K Intensity  ;  Note Intensity, noteStr
	Note/K XX  ;  Note XX, noteStr
	Note/K HH  ;  Note HH, noteStr
	Note/K FF  ;  Note FF, noteStr
	Note/K Xccd  ;  Note Xccd, noteStr
	Note/K Yccd  ;  Note Yccd, noteStr
	SetScale d 0,0,micron, XX,HH,FF
	SetDataFolder fldrSav0
	return noteStr
	// SetDataFolder ::
End



Function/T readSmarterHeader(fullFileName)
	String fullFileName
	Variable refNum
	Open /R/M="Wirescan peak positions from Wenge"/P=home/Z=1 refNum as fullFileName
	if (strlen(S_fileName)<1 || V_flag)
		return ""
	endif
	String noteStr = ""

	noteStr = fetchAndAddNumber2List(refNum,"CCDX0",noteStr,1)	// x-position of un-tilted spot (pixels)
	noteStr = fetchAndAddNumber2List(refNum,"CCDY0",noteStr,1)	// y-position of un-tilted spot (pixels)
	noteStr = fetchAndAddNumber2List(refNum,"psizeX",noteStr,1)	// pixel size (in x-direction) (mm)
	noteStr = fetchAndAddNumber2List(refNum,"psizeY",noteStr,1)	// pixel size (in y-direction) (mm)
	noteStr = fetchAndAddNumber2List(refNum,"x0",noteStr,1)		// x-coord of CCD center, toward door (mm)
	noteStr = fetchAndAddNumber2List(refNum,"y0",noteStr,1)		// y-coord of CCD center, up (mm)
	noteStr = fetchAndAddNumber2List(refNum,"z0",noteStr,1)		// z-coord of CCD center, along beam (mm)
	noteStr = fetchAndAddSring2List(refNum,"ki",noteStr)			// direction of incident beam direction {X(toward door), Y(up) Z(along beam)}
	Close refNum

	// change x0,y0,z0 from Wenge coord system to beam-line coord sytem
	Variable x0,y0,z0
	x0 = str2num(StringByKey("x0",noteStr,"="))
	y0 = str2num(StringByKey("y0",noteStr,"="))
	z0 = str2num(StringByKey("z0",noteStr,"="))
	noteStr = ReplaceNumberByKey("x0", noteStr, x0,"=")			// coord system swap done here (Wenge->beamline)
	noteStr = ReplaceNumberByKey("y0", noteStr, z0,"=")
	noteStr = ReplaceNumberByKey("z0", noteStr, -y0,"=")
//	if (stringmatch(StringByKey("ki",noteStr,"="),"0.0, 0.0, 1.0" ))
//		noteStr = ReplaceStringByKey("ki", noteStr, "0.0, -1.0, 0.0","=")	// all in Wenge coord system
//	endif
	return noteStr
End

Function FindScalingFromVec(vec,threshold,first,stepSize,dimN)
	Wave vec
	Variable threshold			// a change greater than this is intentional (less is jitter)
	Variable &first				// use 	SetScale/P x first,stepSize,"",waveName
	Variable &stepSize
	Variable &dimN

	Variable last,i, Ndim=1
	Duplicate/O vec testdim
	Variable N = numpnts(testdim)
	SetScale/P x 0,1,"", testdim
	Sort testDim, testDim

	Variable delta
	stepSize=Inf
	for (i=0;i<(N-1);i+=1)
		delta = testDim[i+1]-testDim[i]
		if (delta<stepSize && delta>threshold)
			stepSize = delta
		endif
	endfor
	// stepSize is now approximately right

	Variable Nj, stepSum=0
	for (i=0;i<(N-1);i+=1)
		delta = testDim[i+1]-testDim[i]
		if (round(delta/stepSize)==1)
			stepSum += delta
			Nj += 1
		endif
	endfor
	stepSize = stepSum/Nj				// this is now a good average step size

	i = floor(BinarySearchInterp(testDim, testDim[0]+threshold))
	first = faverage(testDim,0,i)
	if (ItemsInList(GetRTStackInfo(0))<=2)
		print stepSize, first
	endif

	i = ceil(BinarySearchInterp(testDim, testDim[N-1]-threshold))
	dimN = round((faverage(testDim,i,N-1)-first)/stepSize + 1)
	KillWaves/Z testdim
	return 0
End



//	$CCDX0		1104.24		//x-position (pixels) of undeformed peak
//	$CCDY0		986.833		//y-position (pixels) of undeformed peak
//	$PsizeX		23.9971		//pixel size (in x-direction) (micron)
//	$PsizeY		24.0730		//pixel size (in y-direction) (micron)
//	$X0		210.757		//x-coord(micron) of un-deformed peak, toward door
//	$Y0		2094.17		//y-coord(micron) of un-deformed peak, downstream
//	$Z0		59498.8		//z-coord(micron) of un-deformed peak, upward
//	$ki			0.0, 0.0, 1.0


//	X-tilt0=1132.61     Y-tilt0=781.16
//	x=1096.87    y=938.30    z=1249.7833 of tilt0
//	incident beam direction  0.0, 0.0, 1.0 (X Y Z) 
//	H=trunc(XH/11)*2*sqrt(2) X=mod(XH,11)*2.0


Proc AboutOLDDataFileHeader()
	Abort "cannot run this!"


the scan was X= 5336µm -> 5356µm,  2 mm step size
H = H0-8.3 -> Ho+20µm,  2.828µm setp size  (2.828 ~ 2*sqrt(2))


   tag						definition
Xccd0				x on the CCD of the un tilted spot
Yccd0				y on the CCD of the un tilted spot
(x0,y0,z0)			absolute (xyz) position of the un-tilted spot on the CCD (but this of course is not what Wenge has in file)
					in the file x0,y0 is the point on the CCD from CCD corner and z0 is height
(k0x,k0y,k0z)		the direction of the incident x-ray beam


// the first 5 lines from the file

X-tilt0=1132.8089     Y-tilt=781.15955
x=1096.87    y=938.30    z=1249.7833 of tilt0
incident beam direction  0.0, 0.0, 1.0 (X Y Z) 
H=trunc(XH/11)*2*sqrt(2) X=mod(XH,11)*2.0
       XH      Z        X-tilt         Y-tilt          Intensity

„Duplicate xh HHfirst, XX
„HHfirst=trunc(XH/11)*2*sqrt(2)
„XX=mod(XH,11)*2.0

XH is now a superfluous wave

„Duplicate XX HH FF
„FF = ZZ/sqrt(2)
„HH = HHfirst + ZZ/sqrt(2)
„KillWaves XH,ZZ

Now our coordinate system in the sample and for displaying is (X,H,F)
End




// ==================================================================
// =============== Start of utility section, very general stuff =====================



Function/T fetchAndAddNumber2List(refNum,tagStr,list,scaleFactor)
	Variable refNum
	String tagStr
	String list
	Variable scaleFactor			// a scale factor to apply to numbers, if NaN, use 1
	scaleFactor = numtype(scaleFactor) ? 1 : scaleFactor

	Variable value
	FSetPos refNum,0
	value = str2num(findTagInFile(refNum,tagStr))*scaleFactor
	if (numtype(value)!=2)
		list = ReplaceNumberByKey(tagStr,list,value,"=")
	endif
	return list
End
Function/T fetchAndAddSring2List(refNum,tagStr,list)
	Variable refNum
	String tagStr
	String list

	String str
	FSetPos refNum,0
	str = findTagInFile(refNum,tagStr)
	if (strlen(str)>1 || (strlen(str)==1 && char2num(str[0])!=32))
		list = ReplaceStringByKey(tagStr,list,str,"=")
	endif
	return list
End

Function/T findTagInFile(fileRef,tagStr)	// returns the 'value' string associated with the tag
	Variable fileRef				// ref to an opened file
	String tagStr				// the tag to search for

	String line
	tagStr = "$"+tagStr
	Variable i,c, tagLen=strlen(tagStr)
	do
		FReadLine fileRef, line
		i = strsearch(line,tagStr,0,2)
		i = (i==0 && char2num(line[tagLen])>32) ? -1 : i		// in case tag is a substring of line
	while (i && strlen(line)>0)

	if (strlen(line)<1)								// end of file
		return ""
	endif

	// find index to first useful character
	i = strlen(tagStr)-1
	do
		i += 1
		c = char2num(line[i])
	while (c<=32 || c==char2num("="))
	line = line[i,Inf]

	if (char2num(line[0])==char2num("'"))		// a quoted string
		i = strsearch(line,"'",1)
		if (i<0)
			Abort "encountered a partially quoted string in tagged data file"
		endif
		line = line[1,i-1]
		return line
	endif

	// find end of value
	i = strsearch(line,"//",0)						// remove possible comment
	line = SelectString(i,line," ",line[0,i-1])
	// i<0		no comment
	// i=0		only a comment, so use 1 space
	// i>0		comment found, remove it

	// remove trailing white space
	i = strlen(line)
	do
		i -= 1
		c = char2num(line[i])
	while(c<=32 && i>0)
	line = line[0,i]
	line = SelectString(strlen(line)==0,line," ")	// empty string means EOF, so make contents at least one space
	return line
End



Static Function normalize(a)	// normalize a and return the initial magnitude
	Wave a
	Variable norm_a = norm(a)
	a /= norm_a
	return norm_a
End