#pragma rtGlobals=1		// Use modern global access method.
#pragma IgorVersion = 5.0
#pragma version = 2.0



Function lattice2recip(lattice,recip)	// creates reciprocal lattice from the real-lattice, returns Vc
	Wave lattice			// the three lattice vectors  (�)
	// The 3 lattice vectors are lattice[0,2][0], lattice[0,2][1], lattice[0,2][2]
	// note that the vector indicies are the first index.  The second index identivies a1, a2, or a3
	Wave recip				// (3x3) matrix to hold the 3 reciprocal lattice vectors (1/�)
	// The 3 reciprocal lattice vectors are recip[0,2][0], recip[0,2][1], recip[0,2][2]
	// note that the vector indicies are the first index.  The second index identivies b1, b2, or b3
	// Also, the recip has the factor of 2� in it already.

	if (!WaveExists(lattice) || !WaveExists(recip))
		DoAlert 0, "in 'lattice2recip' either the lattice or the recip does not exist"
		return NaN
	endif
	if (DimSize(lattice,0)!=3 || DimSize(lattice,1)!=3 || DimSize(lattice,2)!=0)
		DoAlert 0, "in 'lattice2recip' the lattice is not 3x3"
		return NaN
	endif
	Redimension/N=(3,3)/D recip					// force the recip to be (3x3) and double precision

	// make three temporary lattice vectors a1,a2,a3
	String wName=UniqueName("latVec",1,0)
	Make/N=3/D $wName
	Wave a1 = $wName
	wName=UniqueName("latVec",1,0)
	Make/N=3/D $wName
	Wave a2 = $wName
	wName=UniqueName("latVec",1,0)
	Make/N=3/D $wName
	Wave a3 = $wName
	a1 = lattice[p][0]
	a2 = lattice[p][1]
	a3 = lattice[p][2]

	Cross a2,a3								// find each of the three reciprocal lattice vectors, and put into recip
	Wave W_Cross=W_Cross
	Variable Vc = MatrixDot(a1,W_Cross)	// Vc = a1 � (a2 x a3),   Volume of the unit cell
	if (Vc<=0)								// error, return bad values
		recip = NaN
		return NaN
	endif
	recip[][0] = W_Cross[p]				// b1 = 2� * (a2 x a3) / Vc
	Cross a3,a1
	recip[][1] = W_Cross[p]				// b2 = 2� * (a3 x a1) / Vc
	Cross a1,a2
	recip[][2] = W_Cross[p]				// b3 = 2� * (a1 x a2) / Vc
	Killwaves/Z a1,a2,a3,W_Cross
	recip *= 2*PI/Vc
	// The reciprocal lattice has been computed
	// The 3 reciprocal lattice vectors are recip[0,2][0], recip[0,2][1], recip[0,2][2]
	// note that the vector indicies are the first index.  The second index identivies b1, b2, or b3
	// Also, the recip has the factor of 2� in it already
	return Vc
End
//Function test_lattice2recip(strain,axis)
//	Variable strain			// strain applied in the x-direction
//	String axis				// axis for the strain "x", "y", or "z"
//	Variable iaxis = WhichListItem(axis,"x;y;z")
//	if (iaxis<0)
//		Abort "axis must \"x\", \"y\", or \"z\""
//	endif
//
//	Make/N=(3,3)/O/D lattice, recip
//	// define the lattice here
//	// The 3 lattice vectors are lattice[0,2][0], lattice[0,2][1], lattice[0,2][2]
//	// note that the vector indicies are the first index.  The second index identivies a1, a2, or a3
//	lattice = (p==q)								// a cubic lattice
//	//	lattice[][1] = {cos(PI/3),sin(PI/3),0}		// This is just a test, it works
//	lattice[][iaxis] = (1+strain)*lattice[p][iaxis]	// add strain to the lattice
//
//	Variable Vc
//	Vc = lattice2recip(lattice,recip)				// create the recip from lattice
//
//	printf "for a strain of %g%% applied in the %s-direction,    Vc = %g\r",strain*100,axis,Vc
//	printWave(lattice)
//	print ""
//	printWave(recip)
//	print ""
//
////	Testing Section
////	Make/O/T hkls={ "0 0 1","0 1 0","1 0 0","1 1 0","1 1 1","1 2 3","3 -2 1","-2 -4 -6"}
//	Make/O/T hkls={ "1 0 0","0 1 0","0 0 1"}
//	Make/O/N=3/D hkl,qv
//	Variable N = numpnts(hkls)
//	String str
//	Variable hh,kk,ll,i
//	printf "  (hkl)\t\tQ-vector\r"
//	for (i=0;i<N;i+=1)
//		sscanf hkls[i], "%g %g %g",hh,kk,ll
//		hkl = {hh,kk,ll}
//		hkl2Qvec(hkl,recip,qv)
//		printf "(%s)\t\t(%g, %g, %g)\r",hkls[i],qv[0],qv[1],qv[2]
//	endfor
//	KillWaves/Z qv,hkl
//End


Function hkl2Qvec(hkl,recip,qv)	// for an hkl and and a reciprocal lattice, computs the Q of hkl, by Q = recip x hkl
	Wave hkl
	Wave recip
	Wave qv
//	if (!WaveExists(hkl) || !WaveExists(recip) || !WaveExists(qv) || DimSize(qv,0)!=3 || DimSize(qv,1)>0)
	if (!WaveExists(hkl) || !WaveExists(recip) || !WaveExists(qv))
		return 1
	endif
	MatrixMultiply recip,hkl
	Wave M_product=M_product
	qv = M_product
	KillWaves/Z M_product
	return 0
End




Function angleBetweenVectors(a,b)	// returns angle (rad) between vectors a and b
	Wave a,b
	if (DimSize(a,1)>0 || DimSize(b,1)>0)	// 1-d vectors only
		return NaN
	endif
	Variable cosine, angle
	cosine = MatrixDot(a,b)/(norm(a)*norm(b))
//	cosine = max(-1,min(1,dot))			// ensure that truncation does not |cosine| > 1
	if (cosine>=1)
		angle = 0
	elseif (cosine<=-1)
		angle = PI
	else
		angle = acos(cosine)
	endif
	if (!strlen(StringFromList(1,GetRTStackInfo(0))))
		printf "angle between '%s' and '%s' is %g�\r",NameOfWave(a),NameOfWave(b),angle*180/PI
	endif
	return angle
End





Function allowFCC(H)	// flags whether H=(hkl) is an allowed FCC reflection (not mixed)
	Wave H
	if ( (abs(H[0])+abs(H[1])+abs(H[2]))==0 )	// the (000) is not allowed
		return 0
	endif
	Variable sum = mod(abs(H[0]),2)+mod(abs(H[1]),2)+mod(abs(H[2]),2)
	sum = round(sum)
	return (sum==0 || sum==3)
End


Function allowBCC(H)	// flags whether H=(hkl) is an allowed BCC reflection (sum is even)
	Wave H
	return !mod(sum(H,-inf,inf),2)
End


//  This function suplanted by new built in Igor function,  MatrixDot(a,b)
//
//Function dot(a,b)		// vector dot product
//	Wave a,b
//	Variable i=0, n, sum
//	n = numpnts(a)
//	if (n!=numpnts(b))
//		Abort "in dot product, waves must be same "
//	endif
//
//	do
//		sum += a[i]*b[i]
//		i += 1
//	while (i<n)
//	return sum
//End


Function normV(a)	// returns the magnitude of vector a (a is not changed)
	Wave a
	return norm(a)
//	Variable i=0, n, sum=0
//	n = numpnts(a)
//
//	do
//		sum += a[i]*a[i]
//		i += 1
//	while (i<n)
//	return sqrt(sum)
End


Function normalize(a)	// normalize a and return the initial magnitude
	Wave a
	Variable norm_a = norm(a)
	if (norm_a==0)
		return 0
	endif
	a /= norm_a
	return norm_a
//	Variable i=0, n, norm=0
//	n = numpnts(a)
//
//	do
//		norm += a[i]*a[i]
//		i += 1
//	while (i<n)
//	norm = sqrt(norm)
//	if (norm<=0)
//		return 0
//	endif
//
//	i = 0
//	do
//		a[i] /= norm
//		i += 1
//	while (i<n)
//	return norm
End

// supplanted by the Cross operation in Igor 5
//Function cross(a,b,c)	// vector cross product
//	Wave a,b,c
//	c[2] = a[0]*b[1] - a[1]*b[0]
//	c[0] = a[1]*b[2] - a[2]*b[1]
//	c[1] = a[2]*b[0] - a[0]*b[2]
//End


//Function printvec(w)		// print a vector to screen
//	Wave w
//	Variable i=0, n
//	n = numpnts(w)
//	printf "%s = {", NameOfWave(w)
//	do
//		if (i<(n-1))
//			printf "%g, ", w[i]
//		else
//			printf "%g}\r",w[i]
//		endif
//		i += 1
//	while (i<n)
//End