#pragma rtGlobals=1		// Use modern global access method.

// to add to packages, make this addition to the "Analysis->Packages" menu
//
//Menu "Analysis"
//	Submenu "Packages"
//		"hex stacking faults in GaN",Execute/P "INSERTINCLUDE  \"StackingFaults\"";Execute/P "COMPILEPROCEDURES ";Execute/P "initStackingFaults()"
//		help = {"look at scattering from stacking faults in GaN"}
//	End
//End

//	############################################
//							Start of Stacking faults package
//
#include  "Xray", version>=2.2
#pragma IgorVersion = 5.03
#pragma version = 1.0
#pragma ModuleName= StackingFaults

Constant aoGaN = 0.31894							// lattice constants (nm) for GaN
Constant boGaN = 0.31894
Constant coGaN = 0.518614
Constant alphaGaN=90, betaGaN=90, gamGaN=120	// angles for hexagonal

Constant keV = 8.048								// energy for the experiment, CuKa

Menu "Stacking"
	"Calculate an L scan",Lscan($"",NaN,NaN,NaN, NaN,NaN,NaN)
	"-"
	"Make a New Stacking Sequence", MakeNewStacking("",NaN)
	"Edit an Existing Stacking Sequence", EditExistingStacking($"")
	"-"
	"Reset Lattice to Hexagonal GaN",setLattices(aoGaN,boGaN,coGaN,alphaGaN,betaGaN,gamGaN)
	// "Reset Lattice [general]",setLattices(NaN,NaN,NaN,NaN,NaN, NaN)
	"Re-Initialize 'StackingFaults'",initStackingFaults()
End



Function/T Lscan(stacking,reps,h,k,lo,hi,Nintervals)
	Wave/T stacking							// wave with one {A,B,C} style stacking
	Variable reps								// number of repeats of stacking[]
	Variable h,k								// the constant h,k
	Variable lo,hi,Nintervals					// range of the L scan, [lo,hi], with Nintervals+1 points

	String stackName
	if (!WaveExists(stacking) || WaveType(stacking)!=0 || !(reps>0) || numtype(h+k+lo+hi) || !(hi>=lo) || !(Nintervals>=0))
		stackName = SelectString(!WaveExists(stacking) || WaveType(stacking)!=0,NameOfWave(stacking),"")
		reps = !(reps>0) ? 50 : reps
		h = numtype(h) ? 1 : h
		k = numtype(k) ? 0 : k
		lo = !(lo>=0) ? 0.2 : lo
		hi = !(hi>=lo) ? lo+3.8 : hi
		Nintervals = !(Nintervals>=0) ? 0 : Nintervals
		Prompt stackName,"stacking wave",popup,WaveList("*",";","DIMS:1,TEXT:1")
		Prompt reps,"number of times stacking[] repeats"
		Prompt h,"h"
		Prompt k,"k"
		Prompt lo, "low end of L-scan"
		Prompt hi, "high end of L-scan"
		Prompt Nintervals, "number of intervals (1 less than number of points), (=0 to auto choose)"
		DoPrompt "set up L scan",stackName,reps,h,k,lo,hi,Nintervals
		if (V_flag)
			return ""
		endif
		Wave/T stacking = $stackName
		if (Nintervals<=0)							// auto choose Nintervals
			Nintervals = (hi-lo) * numpnts(stacking) * reps * 1.3
		endif
		printf "Lscan(%s,%d,%g,%g, %g,%g,%d)\r",NameOfWave(stacking),reps,h,k,lo,hi,Nintervals
	endif
	if (!WaveExists(stacking))
		DoAlert 0, "stacking wave does not exist in amplitude_hkl()"
		return ""
	endif
	Variable N=numpnts(stacking)
	if (!(N>=1) || WaveType(stacking)!=0)
		DoAlert 0, "the wave '"+NameOfWave(stacking)+"'is not valid, it must be a A,B,C type wave, and at least 1 long"
	endif
	if (!ValidStacking(stacking))
		DoAlert 0,"the wave '"+NameOfWave(stacking)+"' is not a valid stacking"
		return ""
	endif
	if (Nintervals<=0)								// auto choose Nintervals
		Nintervals = (hi-lo) * numpnts(stacking) * reps * 1.3
	endif

	stackName = UniqueName("stackingNumbers",1,0)
	Make/N=(reps*N) $stackName
	Wave wstack=$stackName

	Variable i, iN, A=char2num("A")
	for (i=0;i<reps;i+=1)							// fill in wstack
		iN = i*N
		wstack[iN,iN+N-1] = char2num(UpperStr(stacking[p-iN])) - A
	endfor

	String LscanName=CleanupName(NameOfWave(stacking)+"Lscan",0)
	Make/N=(Nintervals+1)/O $LscanName
	Wave scan = $LscanName
	SetScale/I x lo,hi,"L", scan
	scan = magsqr(amplitude_hkl(h,k,x,wstack))

	KillWaves/Z wstack
	Smooth 3, scan
	return GetWavesDataFolder(scan,2)
End


Static Function/C amplitude_hkl(h,k,l,stacking)		// calculate the scattering amplitude at (hkl), from the stacking seq
	Variable h,k,l									// (hkl)
	Wave stacking									// 0,1,2 style stacking

	if (!WaveExists(stacking))
		DoAlert 0, "stacking wave does not exist in amplitude_hkl()"
		return cmplx(NaN,NaN)
	endif
	Variable N=numpnts(stacking)
	if (!(N>=1) || WaveType(stacking)==0)
		DoAlert 0, "the wave '"+NameOfWave(stacking)+"'is not valid, it must be a 0,1,2 type wave, and at least 1 long"
	endif

	Variable/C amp=cmplx(1,0)
	amp = phaseOfOneStacking(h,k,l,stacking)
	amp *= StructureFactorOfBasis_GaN(h,k,l,keV)

	KillWaves/Z wstack
	return amp
End
//
Static Function/C phaseOfOneStacking(h,k,l,stacking)// returns the phase factor for a particular stacking[]
	Variable h,k,l									// (hkl) used to get kvector
	Wave stacking									// stacking sequence, here it is {0,1,2}, not {A,B,C}

	Variable i, N=numpnts(stacking)				// loop index and size of crystalite along each edge
	Variable/C one=cmplx(1,0)						// convienent short cut
	Wave a=a, b=b, c=c								// the real lattice vectors
	Wave ast=ast, bst=bst, cst=cst					// the reciprocal lattice vectors
	Wave kvec=root:Packages:StackingFaults:kvec	// kvector (contains the 2¹)
	kvec = 2*PI*(h*ast[p] + k*bst[p] + l*cst[p])

	Variable/C basePhase							// phase factor for for all atoms in base plane
	Variable/C pa,pb,pc2,ph						// a phase factor
	pa = one * exp(cmplx(0,MatrixDot(kvec,a)))	// phase factor = exp(i k¥a)
	pb = one * exp(cmplx(0,MatrixDot(kvec,b)))	// phase factor = exp(i k¥b)
//	pc = one * exp(cmplx(0,MatrixDot(kvec,c)))	// phase factor = exp(i k¥c)
	pc2 = one * exp(cmplx(0,MatrixDot(kvec,c)/2))// phase factor = exp(i k¥c/2), only want half phae along c

	Variable/C sumA, sumB							// sum(i=0, N-1, pa^i) and   sum(j=0,N-1,pb) 
	if (pa==one)
		sumA = N									// in case pa==1, then you cannot use this formula for a geometric sum
	else
		sumA = (pa^N - one)/(pa-one)				// sum(i=0, N-1, pa^i)
	endif
	if (pb==one)
		sumB = N
	else
		sumB = (pb^N - one)/(pb-one)
	endif
	basePhase = sumA * sumB						// sum(i=0, N-1, pa^i) * sum(j=0,N-1,pb)

	Wave rvec=root:Packages:StackingFaults:rvec
	Variable/C Aph,Bph,Cph							// phase factor to get A, B, or C position
	Aph = one
	rvec = a*1/3 + b*2/3							// vector displacement from A site to B site
	Bph = one * exp(cmplx(0,MatrixDot(kvec,rvec)))// exp[i*k¥r]
	rvec = -a*1/3 - b*2/3						// vector displacement from A site to C site
	Cph = one * exp(cmplx(0,MatrixDot(kvec,rvec)))// exp[i*k¥r]

	Variable/C phase=cmplx(0,0)					// the total phase factor which will be returned
	Variable/C pci									// exp(i k¥c) ^i
	// sum up the phase factors vertically
	for (i=0;i<N;i+=1)								// cannot ues the geometric sum here because of different Aph,Bph,Cph
		pci = pc2^i
		if (stacking[i]==0)
			phase += basePhase*Aph*pci			// 0 -> A
		elseif (stacking[i]==1)
			phase += basePhase*Bph*pci			// 1 -> B
		elseif (stacking[i]==2)
			phase += basePhase*Cph*pci			// 2 -> C
		else
			Abort "bad value in "+NameOfWave(stacking)+"[], all must be {0,1,2}"
		endif
	endfor
	return phase
End
//
// conventional hexagonal cell is atoms at (0,0,0) and (1/3, 2/3, 1/2)
// base plane has atoms (actually the basis) at all (i*a + j*b + k*c)
//
Static Function/C StructureFactorOfBasis_GaN(h,k,l,keV)// returns structure factor for this bais at (hkl),  fGa + fN*exp(i*k¥r)
	Variable h,k,l									// (hkl) used to get kvector
	Variable keV									// beam energy in keV

	Wave ast=ast, bst=bst, cst=cst					// the reciprocal lattice vectors
	Wave c=c										// real space lattice vector, only need c[], not a[] or b[]
	Wave kvec = root:Packages:StackingFaults:kvec	// kvector (contains the 2¹)
	kvec = 2*PI*(h*ast[p] + k*bst[p] + l*cst[p])
	Variable Q = norm(kvec)						// length of kvec is the Q

	Variable/C fGa = Get_f("Ga",Q,keV)				// atomic structure factor of Gallium
	Variable/C fN = Get_f("N",Q,keV)				// atomic structure factor of Nitrogen
	Wave rN = root:Packages:StackingFaults:rN		// location of N atom, Ga is at origin
	rN = c/4
	Variable/C Fbasis = fGa + fN*exp(cmplx(0,MatrixDot(kvec,rN)))	// structure factor of basis at this k
	return Fbasis
End



Function xxx()
	Wave a=a,b=b,c=c
	Wave lattice=lattice
	Make/N=3/D/O rvec

	rvec = a*1/3 + b*2/3
	rvec = abs(rvec[p]) < 1e-16 ? 0 : rvec[p]
	printWave(rvec)

	rvec = -a*1/3 - b*2/3
	rvec = abs(rvec[p]) < 1e-16 ? 0 : rvec[p]
	printWave(rvec)

	rvec = {-1/3, -2/3, 0}
	MatrixOp/O rout = lattice x rvec
	rout = abs(rout[p]) < 1e-16 ? 0 : rout[p]
	printWave(rout)

	rvec = {1/3, 2/3, 0}
	MatrixOp/O rout = lattice x rvec
	rout = abs(rout[p]) < 1e-16 ? 0 : rout[p]
	printWave(rout)
	KillWaves/Z rvec
End




Function setLattices(ao,bo,co,alpha,beta,gam)	// set both the real lattice, and reciprocal lattice
	Variable ao,bo,co							// three lattice constants (nm)
	Variable alpha,beta,gam						// and three angles (degrees)
	Make/N=3/O/D a,b,c

	if (numtype(ao+bo+co+alpha+beta+gam) || ao<=0 || bo<=0 || co<=0 || alpha<=0 || beta<=0 || gam<=0)
		ao = !(ao>0) ? aoGaN : ao				// default to GaN values
		bo = !(bo>0) ? boGaN : bo
		co = !(co>0) ? coGaN : co
		alpha = !(alpha>0) ? alphaGaN : alpha
		beta = !(beta>0) ? betaGaN : beta
		gam = !(gam>0) ? gamGaN : gam
		Prompt ao, "ao (nm)"
		Prompt bo, "bo (nm)"
		Prompt co, "co (nm)"
		Prompt alpha, "alpha ¡"
		Prompt beta, "beta ¡"
		Prompt gam, "gamma ¡"
		DoPrompt "lattice,"ao,alpha,bo,beta,co,gam
		if (V_flag)
			return NaN
		endif
	endif
	if (numtype(ao+bo+co+alpha+beta+gam) || ao<=0 || bo<=0 || co<=0  || alpha<=0 || beta<=0 || gam<=0)
		return NaN
	endif

	if (ao==bo && ao==co && alpha==90 && beta==90 && gam==90)
		print "this lattice is cubic"
	elseif (ao==bo && alpha==90 && beta==90 && gam==120)
		print "this lattice is hexagonal"
	elseif (ao==bo && ao==co && alpha==beta && beta==gam)
		print "this lattice is rhombohedral"
	elseif (ao==bo && alpha==90 && beta==90 && gam==90)
		print "this lattice is tetragonal"
	elseif (alpha==90 && beta==90 && gam==90)
		print "this lattice is orthorhombic"
	endif

	a = {1,0,0}
	b = {cos(gam*PI/180),sin(gam*PI/180),0}

	Make/N=(3,3)/O/D tempMat
	tempMat[][0] = a[p]
	tempMat[][1] = b[p]
	tempMat[][2] = {0,0,1}
	Make/N=3/O/D tempVec
	tempVec = {cos(beta*PI/180),cos(alpha*PI/180),1}
	tempVec = abs(tempVec[p]) < 1e-16 ? 0 : tempVec[p]
	KillWaves/Z tempOut
	MatrixOp tempOut = Inv(tempMat) x tempVec
	c = tempOut
	normalize(c)
	a *= ao
	b *= bo
	c *= co

	Variable Vc = abc2recip(a,b,c)	// calc recip-lattice and set volume of unit cell
	printf "volume of unit cell = %g (nm^3)\r",Vc
	KillWaves/Z tempMat,tempVec,tempOut

	Make/N=(3,3)/O/D lattice,recipLattice
	lattice[][0] = a[p]
	lattice[][1] = b[p]
	lattice[][2] = c[p]

	Wave ast=ast, bst=bst, cst=cst
	recipLattice[][0] = ast[p]
	recipLattice[][1] = bst[p]
	recipLattice[][2] = cst[p]

	print "created the real space lattice, and the reciprocal space lattice"
	printf "  setLattices(%g, %g, %g,    %g, %g, %g)		// in (nm) and (¡)\r",ao,bo,co,alpha,beta,gam
	return Vc
End
//
Static Function abc2recip(a,b,c)	// given the three lattice vectors, make three recip-lattice vectors, returns Vc
	Wave a,b,c						// three lattice vectors

	Variable Vc	 					// a ¥ (b x c), the triple product, volume of unit cell
	Cross b, c
	Wave W_Cross=W_Cross
	Vc = MatrixDot(W_Cross,a)

	Duplicate/O a,ast,bst,cst
	ast = W_Cross/Vc				// (b x c) / ( a ¥ (b x c) )
	Cross c,a
	bst = W_Cross/Vc				// (c x a) / ( a ¥ (b x c) )
	Cross a,b
	cst = W_Cross/Vc				// (a x b) / ( a ¥ (b x c) )
	KillWaves W_Cross
	return Vc
End


Function/S MakeNewStacking(wName,length)
	String wName				// name to hold the stacking
	Variable length

	Variable N
	if (exists(wName)==1)		// this stacking exists, ask if user wants to edit it
		if (WaveType($wName))
			DoAlert 0, "the wave '"+wName+"' already exists, and it is not a text wave, exiting"
			return ""
		endif
		DoAlert 1, "the stacking '"+wName+"' already exists, edit it?"
		N = numpnts($wName)
		if (V_flag)
			EditExistingStacking($wName)
		endif
		return wName
	endif
	if (CheckName(wName,1) || !(length>0))
		length = !(length>0) ? 9 : length
		Prompt wName,"name to use for stacking wave"
		Prompt length, "number of elements in stacking"
		DoPrompt "stacking wave",wName,length
		if (V_flag)
			return ""
		endif
		length = round(length)
	endif
	if (CheckName(wName,1) || !(length>0))
		if (CheckName(wName,1))
			DoAlert 0, "you selected a new name that already exists, exiting"
		else
			DoAlert 0,"length of "+num2str(length)+" is invalid"
		endif
		return ""
	endif
	Make/N=(length)/T $wName
	wName = EditExistingStacking($wName)
	if (ItemsInList(GetRTStackInfo(0))<2)
		printf "made a new stacking named '%s', that is %d long\r",wName,numpnts($wName)
	endif
	return GetWavesDataFolder($wName,2)
End
//
Function/T EditExistingStacking(stacking)
	Wave/T stacking

	if (!WaveExists(stacking))						// no wave given, ask user to choose from list of existing waves
		String stackName
		Prompt stackName,"stacking wave",popup,WaveList("*",";","DIMS:1,TEXT:1")
		DoPrompt "pick stacking wave",stackName
		if (V_flag)
			return ""
		endif
		Wave/T stacking = $stackName
	endif
	String tName = UniqueName("stackTable",7,0)
	Edit/W=(5,44,258,722)/K=1/N=$tName stacking as "Kill This Table to Continue"
	ModifyTable width(Point)=44,width(stacking)=110
	PauseForUser $tName
	if (!ValidStacking(stacking))
		DoAlert 0, "the stacking '"+NameOfWave(stacking)+"' is not valid, check it"
	endif
	if (numpnts(stacking)<1)
		return ""
	endif
	if (ItemsInList(GetRTStackInfo(0))<2)
		printf "edited '%s', that is %d long\r",NameOfWave(stacking),numpnts(stacking)
	endif
	return GetWavesDataFolder(stacking,2)
End
//
Static Function ValidStacking(stacking)	// returns 1 if stacking[] is and a valid A,B,C type sequence
	Wave/T stacking

	Variable i, N=numpnts(stacking)
	String letter, last=""
	for (i=0;i<N;i+=1)
		letter = UpperStr(stacking[i])
		if (strlen(letter)!=1)
			return 0					// only a single letter is allowed
		elseif (WhichListItem(letter,"A;B;C")<0)
			return 0					// letter must be A, B,or C
		elseif (stringmatch(last,letter))
			return 0					// a repeat
		endif
	endfor

	if (stringmatch(stacking[0],stacking[N-1]))
		return 0						// if first and last match, then you cannot do a repeat
	endif
	return 1			// it must be valid if we got here
End


Function initStackingFaults()			// initialize the package
	ElementDataInitPackage()
	if (!DataFolderExists("root:Packages"))
		NewDataFolder root:Packages
	endif
	if (!DataFolderExists("root:Packages:StackingFaults"))
		NewDataFolder root:Packages:StackingFaults
	endif
	Make/N=3/O/D root:Packages:StackingFaults:kvec
	Make/N=3/O/D root:Packages:StackingFaults:rvec
	Make/N=3/O/D root:Packages:StackingFaults:rN
	setLattices(aoGaN,boGaN,coGaN,alphaGaN,betaGaN, gamGaN)
End
//
//							End of stacking faults package
//	############################################