#pragma rtGlobals=1		// Use modern global access method.
#pragma IgorVersion = 4.0
#pragma version = 1.0
#include "vector-math"
#include  <TransformAxis>



 Menu "Umwegs"
	"make picture", web("")
	"layout graph", LayoutTopGraph()
	"table of hkl", Table_hkl()
End

Window Table_hkl() : Table
	PauseUpdate; Silent 1		// building window...
	String fldrSav= GetDataFolder(1)
	SetDataFolder root:temp:
	Edit/W=(5,42,197,601) hkl
	ModifyTable width(Point)=32,width(hkl)=38
	ModifyTable font="Helvetica",size=9
	SetDataFolder fldrSav
EndMacro


Window LayoutTopGraph() : Layout
	PauseUpdate; Silent 1		// building window...
	String cmd = "Layout/C=1/W=(5,42,382,521) "+WinName(0,1)+"(39,46,556,754)/O=1/F=0"
	Execute cmd
	Layout_Corner_Labels_Style_()
	DoWindow/C $("Layout"+WinName(0,1))
EndMacro




Function web(name)
	String name="Si"			// string for material
	Variable LAMin=0.8			// minimum lamda (�)
	Variable LAMax=1.2		// maximum lamda (�)
	Variable phimin=0			// minimum phi
	Variable phimax=30		// maximum phi
	if (strlen(name)<1)
		name = "Si"
	endif
	Prompt name, "Material", popup, "Al;Cu;Fe;Ge;Si;W;SrTiO3;Fe2O3;Ho"
	Prompt LAMin, "Smallest wavelength (Angstroms)"
	Prompt LAMax, "Largest wavelength (Angstroms)"
	Prompt phimin, "Minimum phi in (degrees)"
	Prompt phimax, "Maximum phi in (degrees)"
	DoPrompt "for Web", name, LAMin, LAMax, phimin,phimax
	if (V_flag)
		return 1				// Abort "user terminated"
	endif
	if (exists("TransformAxisPanelInitGlobals"))
		TransformAxisPanelInitGlobals()
	endif

	phimin *= pi/180.
	phimax *= pi/180.

//	This program computes the multiple reflections
//	for a diamond structure crystal at many wavelengths
//
//		x/m web
//		Ge
//		1.2
//		3
//		0
//		180
//		2,0,0
//		)
	Variable DIA=227			// Diamond structure
	Variable MAXhkls=8000

	Variable h0=2,k0=2,l0=2	// (hkl) of main reflection
	Variable struct				// defines structure, fcc, bcc, diamond...

	Variable tick0,tick1,tick2
	Variable M1,M2,M3		// (hkl) of reference vector
	Variable hklnum			// actual number of (h'k'l')'s in hkl
	Variable long				// number of points in a particular line
	Variable j					// loop index

	Prompt h0,"h"
	Prompt k0,"k"
	Prompt l0,"l"
	DoPrompt "Give (hkl) for the reflection", h0,k0,l0
	if (V_flag)
		return 1				// Abort "user terminated"
	endif

	if (h0==6 && k0==2 && l0==2)
		M1 = 0								// reference vector for (622)
		M2 = 1								//  is (0 1 -1)
		M3 = -1
	elseif (h0==4 && k0==4 && l0==2)
		M1 = 1								// reference vector for (442)
		M2 = -1								//  is (1 -1 0)
		M3 = 0
	elseif (k0==0 && l0==0)
		M1 = 0								// reference vector for (h00)
		M2 = 1								//  is (0 1 0)
		M3 = 0
	elseif (h0==k0 && k0==l0)
		M1 = 1								// reference vector for (hhh)
		M2 = -1							//  is (1 -1 0)
		M3 = 0
	else
		Prompt M1,"h"
		Prompt M2,"k"
		Prompt M3,"l"
		DoPrompt "Give (h'k'l') perpendicular to (hkl)", M1,M2,M3
		if (V_flag)
			return 1						// Abort "user terminated"
		endif
	endif
	tick0 = ticks

	String fldrSav= GetDataFolder(1)
	NewDataFolder/O/S root:temp
//	SetDataFolder root:temp

	struct = matter(name)
	Wave recip=recip

	Make/N=3/O Gh				// reciprocal lattice vector for h0,k0,l0
	Make/N=3/O Gm			// reciprocal latt vector for M1,M2,M3
	Make/N=3/O Gn				// perpendicular to Gh and Gm
	Make/N=3/O hklOne		// used to call shoot

	Make/N=2000/O lambdas	// wavelength in Angstroms as array
	Make/N=2000/O phis		// ordered output array in degrees
	Make/N=(MAXhkls,3)/O hkl// list of the (h'k'l')'s

	getG(h0,k0,l0,recip,Gh)						// get Gh
	getG(M1,M2,M3,recip,Gm)						// get Gm
	cross (Gh,Gm,Gn)								// get Gn = Gh X Gm
	cross (Gn,Gh,Gm)								// set Gm = Gn X Gh

	print " calling list..."
	tick1 = ticks
	hklnum = list(LAMin,hkl,h0,k0,l0,Gh)			// get hkl, array of (h'k'l')'s
	print " from list found ",hklnum,"possible reflections,  it took ",(ticks-tick1)/60.15,"sec"
	Redimension/N=(hklnum,3) hkl

	Display /W=(5,42,520,600)
	PauseUpdate
	String str
	sprintf str, "%s (%d %d %d)\rreference  (%d %d %d)", StringByKey("name",note(recip)), h0,k0,l0,M1,M2,M3
//	TextBox/C/N=text0/F=0/S=3/A=LT str
	TextBox/C/N=text0/F=0/S=3/A=MT/Y=1/E str
	String savedDF= GetDataFolder(1)
	NewDataFolder/O/S lines

	String phiName, lambdaName
	Variable iline=0
	print " starting to shoot lines..."
	tick2 = ticks
	for (j=0; j<hklnum; j += 1)				// one line for each h'k'l'
//		printf  "working on (h'k'l') = (%d,%d,%d), no. %d of %d\r", hkl[j][0],hkl[j][1],hkl[j][2],j,hklnum
		hklOne = hkl[j][p]
		long = shoot(lambdas,phis,LAMin,LAMax,phimin,phimax,hklOne)
		if (long>=2)
			phiName = "phi"+num2istr(iline)
			lambdaName = "lambda"+num2istr(iline)
			Execute "Make/O/N=("+num2istr(long)+") "+phiName+", "+ lambdaName
			Wave phiW = $phiName
			Wave lambdaW = $lambdaName
			phiW = phis
			lambdaW = lambdas
			AppendToGraph lambdaW vs phiW
			if (struct==DIA && mod(hkl[j][0],2)==0)
				if ( mod(hkl[j][0]+hkl[j][1]+hkl[j][2],4)!=0 || mod(H0-hkl[j][0]+K0-hkl[j][1]+L0-hkl[j][2],4)!=0 )
					ModifyGraph lstyle($lambdaName)=2
				endif
			endif
			iline += 1
		else
			hkl[j][0] = NaN
			hkl[j][1] = NaN
			hkl[j][2] = NaN
		endif
		if (mod(j,100)==0)
			print "       done with point ",j
		endif
	endfor
	print "to draw all the lines took ",(ticks-tick2)/60.15,"sec"
	removeNaNsRows(hkl)

	ModifyGraph/Z gfSize=18, lowTrip=0.001, mirror=1
	ModifyGraph/Z tick=0, minor=1, standoff=0
	ModifyGraph/Z rgb=(0,0,30000), lsize=0.2
	SetAxis bottom 180.*phimin/pi, 180.*phimax/pi
	SetAxis left LAMin, LAMax
	Label bottom "\\F'Symbol'f\\F]0�"
	Label left "\\F'Symbol'l\\F]0  (�)"
	ModifyGraph lblMargin(left)=4
	SetDataFolder savedDF
	SetDataFolder fldrSav
	print "total execution time = ",(ticks-tick0)/60.15,"sec"
	DoAlert 1, "Add an energy axis on right?"
	DoUpdate
	if (V_flag==1)
		AddEnergyAxis()
	endif
End

Function AddEnergyAxis()
	PauseUpdate
	Make/O hcw={12.39841857,0}
	String theAxis = "left"
	Wave/Z coefWave = hcw
	Variable tickDensity = 3	
	Variable tickSep = 5
	SetupTransformMirrorAxis(WinName(0, 1), "left", "TransAx_ModifiedReciprocal", coefWave, tickDensity, 1, tickSep, 1)
	ModifyGraph mirror(left)=0,mirror(MT_left)=0, mirror(bottom)=1
	ModifyGraph margin(left)=72,margin(bottom)=72,margin(top)=72,margin(right)=72,gfSize=18
	ModifyGraph lblMargin(MT_left)=4, lblPos(MT_left)=68, freePos(MT_left)=0
	Label MT_left "Energy  (keV)"
end




Function shoot(lambdas,phis,LAMin,LAMax,phimin,phimax,hkl)
//	This routine finds the intersection of reflection hkl with the sphere
//	the wavelengths lambdas(i) and stores results in phis (degrees)
	Wave lambdas
	Wave phis						// ordered output array for input
	Variable LAMin, LAMax
	Variable phimin, phimax
	Wave hkl						// length 3
	Variable long

	Variable re=2.82E-5
	Variable steps=150				// was 600, number of division in lamda
	Variable slope_min=.15			// minimum scaled slope
	Variable slope_max=5			// maximum   "      "

	Wave Gh=::Gh					// reciprocal lattice vector for H0,K0,L0
	Wave Gm=::Gm					// reciprocal lattice vector for M1,M2,M3
	Wave Gn=::Gn					// perpendicular to Gh and GM
	Wave recip=::recip				// recip lattice unit cell column vectors
	Make/N=3/O G					// reciprocal lattice vector for hkl
	Variable LAMDA
	Variable lamda_2,lamda_1,lamda_0
	Variable dLAM0
	Variable phinew
	Variable phi_2,phi_1
	Variable slope					// |slope| of curve (|dlamda/dphi|)
	Variable slope_scale			// slope for a diagonal
	Variable isign					// sign of delta lamda
	Variable enough					// cheks for enough for quadratic
	Variable H0
	Variable ROOT					// number of current root wanted
	Variable STRUCT
	Variable rho					// electrons per Angstrom^3
	rho = NumberByKey("rho",note(recip))

	long = 0
	getG(hkl[0],hkl[1],hkl[2],recip,G)
	slope_scale = (LAMax-LAMin) / (phimax-phimin)
	dLAM0 = (LAMax-LAMin) / steps
	enough = 0
	isign = 1
	slope = slope_min
	LAMDA = LAMin
	lamda_2 = LAMin
	lamda_1 = LAMin
	lamda_0 = LAMin
	phinew = phimax
	phi_1 = phimin
	phi_2 = phimin
	ROOT = 0								// start on first root
	Variable delta, R, RMN
	Variable iflag
	Variable/G phiA,phiB
	do
		delta = LAMDA^2 * re * 4 * rho / pi
		R = 2*pi / LAMDA * (1-DELTA)
		RMN = sqrt( max(R*R - normV(Gh)^2/4, 0) )
		iflag = IAphi(G,Gh,Gm,Gn,R,RMN)

		if (IFLAG>=1)								// get value of current root
			if (iflag==1 || iflag==2)				// found a root
				enough = enough + 1
				phi_2 = phi_1
				phi_1 = phinew
				lamda_2 = lamda_1
				lamda_1 = lamda_0
				lamda_0 = LAMDA
			endif
			IF (IFLAG==1)
				phinew = phiA
			endif
			IF (IFLAG==2)
				phinew = (ROOT==0) ? phiA : phiB
			endif
			IF ((phinew>=phimin) && (phinew<=phimax) && LAMDA<=LAMax)
				lambdas[long] = LAMDA
				phis[long] = 180*phinew/pi
				long += 1
//				print " * lamda = ",lamda,"     phi = ",phinew*180/pi,"     slope =",slope
		       else
//				print "   lamda = ",lamda,"     phi = ",phinew*180/pi,"     slope = ",slope
			endif
		endif

//          if (abs(phinew-phi_1)<1.E-15)
//             slope = slope_max
		if (enough<3)
			slope = slope_min
		else
			slope = quad(lamda_0,lamda_1,lamda_2,phinew,phi_1,phi_2)
//             slope = (lamda_0-lamda_1) / (phinew-phi_1)
			slope = abs(slope/slope_scale)			// rescale to plot diagonal
			slope = max(slope,slope_min)			// minimum slope expansion
			slope = min(slope,slope_max)			// and points not too far apart
		endif

		IF (IFLAG<=0)						// above the peak, go down
			ROOT = 1						// switch to second root
			isign = -1
		elseif (IFLAG==1)					// exactly at the peak, go down
			ROOT = 1						// switch to second root
			isign = -1
		endif

		LAMDA = LAMDA + isign * dLAM0 * slope
		while (!(LAMDA<LAMin || long>=2000))

	KillWaves/Z hklOne
	return long
End



Function list(lambda,hkl,h0,k0,l0,Gh)
//		This routine finds all possible h'k'l' for a given wavelength and
//		and returns them in the array hkl.
	Variable lambda			// wavelength (�)
	Wave hkl				// list of the hkl's (maxx,3) to be filled
	Variable H0, K0, L0		// hkl of main reflection
	Wave Gh				// reciprocal lattice vector for H0,K0,L0

	Wave recip=recip			// recip lattice unit cell column vectors
	Make/N=3/O G_temp_		// reciprocal lattice vector for one hkl

	Variable h,k,l
	Variable struct				// Defines the lattice, FCC, BCC, DIA,...
	Variable rho				// electrons per Angstrom^3

	rho = NumberByKey("rho",note(recip))
	struct = NumberByKey("struct",note(recip))
	Variable hklnum=DimSize(hkl, 0)					// maximum number of hkl's
	Variable delta = lambda^2 * 2.82E-5 * 4. * rho / pi
	Variable R = 2*pi / lambda * (1-delta)				// R is radius of Ewald sphere with index correction
	Variable HMAX = 2*R/getG(0,1,0,recip,G_temp_)	// compare to length of (010)

	Variable RMN = sqrt(R^2 - normV(Gh)^2/4)
	if (numtype(RMN))
		return 0
	endif
	Variable AHP										// magnitude of G_temp_
	Variable hklLimit = round(HMAX)
	Variable index = 0
	Variable/G phiA,phiB

	for (h=-hklLimit; h<=hklLimit; h+=1)
		for (k=-hklLimit; k<=hklLimit; k+=1)
			for (l=-hklLimit; l<=hklLimit; l+=1)
				if (index < hklnum)
					AHP = getG(h,k,l,recip,G_temp_)
					if ( (AHP<=HMAX && AHP>0) && ((h!=0 || k!=0 || l!=0) && (h!=H0 || k!=K0 || l!=L0)) && (allow(h,k,l,struct) && allow(H0-h,K0-k,L0-l,struct)) )
						if (IAphi(G_temp_,Gh,Gm,Gn,R,RMN)>=1)
							hkl[index][0] = h
							hkl[index][1] = k
							hkl[index][2] = l
							index += 1
						endif
					endif
				endif
			endfor
		endfor
	endfor
	KillWaves/Z G_temp_
	return index
End



Function IAphi(G,Gh,Gm,Gn,R,RMN)
//	This routine calculates the phi value of a given hkl.  If necessary,
//	it also determies how many times hkl intersects the Ewald sphere.
	Wave G	 				// reciprocal lattice vector of h'k'l'
	Wave Gh				// reciprocal lattice vector for H0,K0,L0
	Wave Gm				// reciprocal lattice vector for M1,M2,M3
	Wave Gn				// perpendicular to GH and GM
	Variable R				// radius of Ewald sphere
	Variable RMN

	Variable/G phiA, phiB	// global to pass back the answer

	Variable HN, HM
	Variable HH
	Variable Rcut

	Variable iflag = 2
	HH = MatrixDot(G,GH) / normV(GH)
	Rcut = sqrt(R^2 - (HH-normV(GH)/2)^2)
	if (numtype(Rcut))
		return 0
	endif
	Variable RLIMU = Rcut + RMN
	Variable RLIML = abs(Rcut-RMN)
	HM = MatrixDot(G,Gm) / normV(Gm)
	HN = MatrixDot(G,Gn) / normV(Gn)
	Variable HMN = sqrt(normV(G)^2-HH^2)
	HMN = numtype(HMN) ? 0 : HMN
	if (HMN>RLIMU || HMN<RLIML)
		return 0
	endif
	Variable cosB = (HMN == 0) ? 0 : (HMN^2 + RMN^2 - Rcut^2) / (2.*HMN*RMN)
	if (abs(cosB) > 1.)
		return 0
	endif
	Variable beta = acos(cosB)
	Variable eta = (HM==0 && HN==0) ? 0 : atan2(HN,HM)
	iflag = (abs(cosB)==1 || HMN==0) ? 1 : iflag
	phiA = range(eta + beta)
	phiB = range(eta - beta)
//	print "phiA=",phiA,"    phiB=",phiB
	return iflag
End

Function range(x)
	Variable x
	Variable range
	range = mod(x,2*pi)
	range = (range <= -pi) ? range + 2*pi : range
	range = (range > pi) ? range - 2*pi : range
	return range
End



Function matter(name)
//	Routine to get the reciprocal lattice and structure (struct) for the
//	program LAUE.  If none is given the program requests it.
//	the name, struct, rho are passed back as a note in recip
	String name					// Label for plot, name of material

	Variable FCC=255			// Face Centered Cubic
	Variable DIA=227			// Diamond structure
	Variable BCC=229			// Body Centered Cubic
	Variable SC=221			// Simple Cubic
	Variable RHOM=167		// Rhombehedral Cr2o3 structure
	Variable HEX=194			// Hexagonal Structure
	Variable PEROV=221		// Perovskite, CaTiO3 structure

	Make/N=(3,3)/O recip		// The reciprocal lattice, column vectors
	Make/N=3/O d0,d1,d2		//   "    direct           "            "           "

	Variable struct				// Defines the lattice, FCC, BCC, DIA,...
	Variable rho				// electrons per Angstrom^3
	String longName			// long version of name
	Variable ao=-1				// Lattice const for a cubic lattice
	Variable norm				// running sum, used to normalize direct
	Variable SINthe				// Sin(theta), theta is just some angle
	Variable COSthe				// Cos(theta),    "     "    "        "     "
	Variable Z					// total number of electrons per cell
	Variable Vc					// volume of unit cell

	d0 = {1,0,0}
	d1 = {0,1,0}
	d2 = {0,0,1}

	if (stringmatch(name,"Al"))					// Al     Aluminum
		ao = 4.049
		struct = FCC
		longName = "Aluminum"
		Z = 4 * 13
	elseif (stringmatch(name,"Cu"))				// Cu     Copper
		ao = 3.62
		struct = FCC
		longName = "Copper"
		Z = 4 * 29
	elseif (stringmatch(name,"Fe"))				// Fe     alpha-Iron
		ao = 2.867
		struct = BCC
		longName = "alpha-Iron"
		Z = 2 * 26
	elseif (stringmatch(name,"Ge"))				// Ge     Germanium
		ao = 5.63
		struct = DIA
		longName = "Germanium"
		Z = 8 * 32
	elseif (stringmatch(name,"Si"))				// Si     Silicon
		ao = 5.43072
		struct = DIA
		longName = "Silicon"
		Z = 8 * 14
	elseif (stringmatch(name,"W"))				// W      Tungsten
		ao = 3.163
		struct = BCC
		longName = "Tungsten"
		Z = 2 * 74
	elseif (stringmatch(name,"SrTiO3"))			// Perovskite
		ao = 3.9051
		struct = PEROV
		longName = "SrTiO3"
		Z = 4 * ( 38 + 22 + 3*8 )
	elseif (stringmatch(name,"Fe2O3"))			// Fe2O3  Hematite
		ao = 5.42
		struct = RHOM
		longName = "Hematite"
		Z = 2 * ( 2*26 + 3*8 )
		SINthe = SQRT(2./3.*.4297642)			// 0.4297642 = 1-COS(ALPHA)
		COSthe = SQRT(1. - 2./3.*.4297642)
		d0 = {SINthe, 0, COSthe}
		d1 = {-SINthe/2., sqrt(3.)*SINthe/2., COSthe}
		d2 = {-SINthe/2., -SQRT(3.)*SINthe/2., COSthe}
	elseif (stringmatch(name,"Ho"))				// Holmium
		struct = HEX
		longName = "Holmium"
		Z = 2 * 67
		d0 = {3.564, 0, 0}							// a-axis
		d1 = {cos(2.*pi/3.), sin(2.*pi/3.), 0.}		// b-axis
		d1 *= 3.564
		d2 = {0, 0, 5.631}							// c-axis
	else
		Abort "matter -- ERROR, Material not listed"
	endif

	if (ao>0.)										// For a single lattice const
		d0 *= ao / sqrt(MatrixDot(d0,d0 ))			// make d0 ao long
		d1 *= ao / sqrt(MatrixDot(d1,d1 ))			// make d1 ao long
		d2 *= ao / sqrt(MatrixDot(d2,d2 ))			// make d2 ao long
	endif

	Make/N=3/O vector__							// vector used to get Vc
	cross(d1,d2,vector__)							// Vc = a1 . a2 X a3
	Vc = MatrixDot(d0, vector__ )

//	Now calculate the reciprocal lattice, NOTE: recip(i,j) has a 2pi in it!
	Make/N=3/O r0, r1, r2
	cross(d1,d2,r0)
	cross(d2,d0,r1)
	cross(d0,d1,r2)
	recip[][2] = r2[p]
	recip[][1] = r1[p]
	recip[][0] = r0[p]
	recip *= 2*pi/Vc
	KillWaves/Z r0,r1,r2
	KillWaves/Z vector__, d0,d1,d2

	rho = Z / Vc										// number of electrons per Ang^3
	if (Vc<=0.)
		Abort " MATTER -- ERROR, Vc is negative or zero"
	endif
	String str
	sprintf str, "name:%s;struct:%d;rho:%g;", longName, struct, rho
	Note recip, str
	return struct
end



Function allow(h,k,l,struct)	//	This function is true if the (hkl) is an allowed reflection for struct
	Variable h,k,l				// Crystal index (as in hkl)
	Variable struct				// Defines the lattice, FCC, BCC, DIA,...

	Variable FCC=255			// Face Centered Cubic
	Variable DIA=227			// Diamond structure
	Variable BCC=229			// Body Centered Cubic
	Variable SC=221			// Simple Cubic
	Variable RHOM=167		// Rhombehedral Cr2o3 structure
	Variable HEX=194			// Hexagonal Structure
	Variable PEROV=221		// Perovskite, CaTiO3 structure

	Variable mod2 = mod(abs(H),2) + mod(abs(K),2) + mod(abs(L),2)

	if (h==0 && k==0 && l==0)
		return 1
	elseif (struct==FCC && (mod2==3 || mod2==0))									// FCC
		return 1
	elseif (struct==BCC && !mod(H+K+L,2))											// BCC
		return 1
	elseif (struct==DIA && ((mod2==3 || (mod2==0 && mod(H+K+L,4)==0))))		// Diamond
		return 1
	elseif (struct==DIA && abs(H)==2 && abs(K)==2 && abs(L)==2)					// special, 222 allowed
		return 1
	elseif (struct==SC)																	// Simple cubic
		return 1
	elseif(struct==RHOM && !(H==K && mod(L,2) || H==L && mod(K,2) || K==L && mod(H,2) )  )	// Cr2O3
		return 1
	elseif (struct==HEX && (!mod(L,2) || mod(L,2) && mod(H+2*K,3)!=0))			// Hexagonal
		return 1
	elseif (struct==PEROV)																// Perovskite
		return 1
	else																					// An extinction
		return 0
	endif
end



Function getG(h,k,l,recip,Gvec)		// fills Gvec and returns |Gvec|
	Variable h,k,l		// integral (hkl)
	Wave recip			// reciprocal lattice vectors
	Wave Gvec			// desired reciprocal lattice vector
	Make/O M_product = {h,k,l}
	MatrixMultiply recip, M_product
	Gvec = M_product
	KillWaves/Z M_product
	return sqrt(Gvec[0]^2+Gvec[1]^2+Gvec[2]^2)
End


Function quad(lamda_0,lamda_1,lamda_2,phi_0,phi_1,phi_2)
	Variable lamda_0,lamda_1,lamda_2
	Variable phi_0,phi_1,phi_2

	Variable l0, l1			// rescaled lambdas
	Variable p0, p1			//    "     phis
	Variable a, b			// coefficients
	Variable thresh=1.e-20

	if (abs(phi_0-phi_2)<thresh || abs(phi_1-phi_2)<thresh)
		return 0
	endif
	l0 = ( lamda_0-lamda_2 ) / ( phi_0-phi_2 )
	l1 = ( lamda_1-lamda_2 ) / ( phi_1-phi_2 )
	p0 = ( phi_0^2-phi_2^2 ) / ( phi_0-phi_2 )
	p1 = ( phi_1^2-phi_2^2 ) / ( phi_1-phi_2 )

	if (abs(p0-p1)<thresh)
		return 0
	endif
	a = (l0-l1) / (p0-p1)
	b = l0 - a*p0
	return (2*a*phi_0 + b)
end




Function removeNaNsRows(mat)
//	used to remove NaN's from the hkl matrix
	Wave mat
	if (Dimsize(mat,1)<=0 || Dimsize(mat,2)>0)
		Abort "removeNaNs only works on 2-d data"
	endif

	Variable N=Dimsize(mat,0)
	Variable M=Dimsize(mat,1)
	Variable i=0
	do
		if (numtype(mat[i][0]))
			N  -= 1
			mat[i,N-1][0,M] = mat[p+1][q]
			Redimension/N=(N,M) mat
		else
			i += 1
		endif
	while (i<N)
End