#pragma rtGlobals=1		// Use modern global access method.
#pragma version = 1.5
#include "FWHM"

// consider also the following:
// #include "JelliumResponse"


Menu "inelastic"
	"Q of scan", Q_of_scan(NaN)
	"Q->kF",  pQau2kF(-1)
	"QkF->au", pQkF2au(-1)
	"Q->au", pQang2au(-1)
	"Qau->", pQau2ang(-1)
	"hkl->au", phkl2au(-1)
	"hkl->", phkl2ang(-1)
	"hkl->kF", phkl2kF(-1)
	"f-sum rule",f_sumRule(NaN)
	"List all allowed reflections",allowedReflections(75,NaN)
	"List ALL reflections",AllReflections(0,NaN)
	"info about 1 reflection",OneReflections("",NaN)
End


Proc addMoreAnnotation2specPlot(scanNum)
	Variable scanNum
	ModifyGraph mode=4,marker=19,msize=2,lstyle=1,rgb=(1,4,52428)
	Variable Q = Q_of_scan(scanNum)
	AppendText "Q = "+num2str(Q)+" (1/nm)"
	AppendText "Q = "+num2str(Q/10)+" (1/)"
EndMacro



Function f_sumRule(Q)					// f-sum rule = n * đ * (e^2/aB) * (Q*aB)^2s   (eV/nm^3)
	Variable Q							// Q in 1/
	Variable printIt=0
	if (numtype(Q) || Q<0)
		Q = 1.75
		Prompt Q, "Q (1/)"
		DoPrompt "Give the Q (1/)", Q
		if (V_flag)
			return NaN
		endif
		printIt=1
	endif
	if (numtype(Q) || Q<0)
		return NaN
	endif
	Variable f_sum
	NVAR aB=root:Packages:Constants:aB
	NVAR Rydberg=root:Packages:Constants:Rydberg		// 2*Rydberg = Hartree = e^2 / aB
										// Rydberg = (alpha^2) * me_kg * c / (2*h_Js) * h_Js*c / e  [eV]
	Variable ao = NumVarOrDefault("root:Packages:inelasticUtility:ao_A", 4.048 )	// current lattice const. ()
	Variable Ve = NumVarOrDefault("root:Ve", ao^3/12 )// volume of one electron
	Variable Ve_Al = 4.048^3/12						// Ve for Al
	Variable not_aluminum = (abs(Ve-Ve_Al)>.01)

	f_sum = (1/Ve)*PI*(2*Rydberg)*(aB*Q)^2		// n = 1/Ve, electron density (electrons/^3)
	if (printIt || not_aluminum || ItemsInList(GetRTStackInfo(0))<2)
		printf "f-sum(Q=%g (1/), %g (^3/electron)) = %g (eV/^3)", Q, Ve, f_sum
		printf "\t\tf-sum(Q=%g (1/), %g (^3/electron)) = %g (eV/nm^3)", Q, Ve, f_sum*1000
		if (not_aluminum)
			printf "        this is NOT Aluminum"
		endif
		print ""
	endif
	return f_sum*1000								// return answer in nm, not 
End



Function QkF_of_scan(scanNum)
	Variable scanNum
	String funcName="specInfo"
	if (exists("specInfo")!=6)
		funcName="specInfoProtoTypeFunc"
		Abort "need to load the 'spec files' Package if you want to work with spec scans"
	endif
	FUNCREF specInfoProtoTypeFunc refspecInfo = $funcName

	if (scanNum<1)
		scanNum=NumVarOrDefault("root:lastScan", 0)
		Prompt scanNum,"spec scan number"
		DoPrompt "scan number", scanNum
	endif
	if (V_Flag)
		Abort
	endif
	if (scanNum<1)
		Abort "scan number must be >0"
	endif
	Variable theta = refspecInfo(scanNum,"delta")/2
	if (numtype(theta))
		theta = refspecInfo(scanNum,"2theta")/2					// try again if kappa not psic
	endif
	NVAR hc = root:Packages:Constants:hc
	Variable QkF=theta2QkF(theta,hc/refspecInfo(scanNum,"Eo"))	// was DCM_energy, instead of Eo
	return QkF
End
Function specInfoProtoTypeFunc(scanNum,key)	
	Variable scanNum
	String key
	print "bad call, scanNum=",scanNum,"     key=",key
End


Function Q_of_scan(scanNum)
	Variable scanNum
	String funcName="specInfo"
	if (exists("specInfo")!=6)
		funcName="specInfoProtoTypeFunc"
		Abort "need to load the 'spec files' Package if you want to work with spec scans"
	endif
	FUNCREF specInfoProtoTypeFunc refspecInfo = $funcName

	Variable doPrint=0
	if (scanNum<1 || numtype(scanNum))
		scanNum=NumVarOrDefault("root:Packages:spec:lastScan", 0)
		Prompt scanNum,"spec scan number"
		DoPrompt "scan number", scanNum
		doPrint=1
	endif
	if (V_Flag)
		Abort
	endif
	if (scanNum<1)
		Abort "scan number must be >0"
	endif
	Variable theta = refspecInfo(scanNum,"delta")/2	// usual for psic
	if (numtype(theta))
		theta = refspecInfo(scanNum,"Two-theta")/2
	endif
	if (numtype(theta))
		theta = refspecInfo(scanNum,"2theta")/2		// usual for kappa
	endif
	NVAR hc = root:Packages:Constants:hc
	Variable wavelength=hc/refspecInfo(scanNum,"Eo")
	wavelength = numtype(wavelength) ? specInfo(scanNum,"wavelength") : wavelength
	Variable Q=theta2Q(theta,wavelength)*10			// was DCM_energy, instead of Eo
	if (numtype(Q))
		Q = specInfo(scanNum,"Q_nm")
		if (numtype(Q))
			Q = specInfo(scanNum,"Q_angstrom")*10
		endif
	endif
	if (doPrint || ItemsInList(GetRTStackInfo(0))<2)
		printf "Q(scan %d) = %g(1/) = %g(1/nm)\r", scanNum, Q/10,Q
	endif
	return Q											// return Q (1/nm)
End




Function theta2QkF(theta,lambda)
	Variable theta								// Bragg angle (degrees)
	Variable lambda								// wavelength ()
	NVAR kF_invA=root:Packages:inelasticUtility:kF_invA	// kFermi (1/)
	// Q = 4đ*sin(theta)/lambda
	return 4*PI*sin(theta*PI/180)/lambda / kF_invA
End


Function theta2Q(theta,lambda)
	Variable theta								// Bragg angle (degrees)
	Variable lambda								// wavelength ()
	// Q = 4đ*sin(theta)/lambda
	return 4*PI*sin(theta*PI/180)/lambda
End


Function hkl2kF(hkl)
	Variable hkl								// magnitude of G in hkl units
	NVAR ao=root:Packages:inelasticUtility:ao_A	// lattice const ()
	NVAR kF_invA=root:Packages:inelasticUtility:kF_invA	// kFermi (1/)
	return  (2*PI * hkl/ao /kF_invA)
End

Function Qau2kF(Qau)
	Variable Qau								// q in atomic units
	NVAR aB=root:Packages:Constants:aB
	NVAR kF_invA=root:Packages:inelasticUtility:kF_invA	// kFermi (1/)
	return Qau/(kF_invA*aB)					// kF_invA is 1/->kF, aB is -> au
End

Function QkF2au(QkF)
	Variable QkF								// q in kFermi
	NVAR aB=root:Packages:Constants:aB
	NVAR kF_invA=root:Packages:inelasticUtility:kF_invA	// kFermi (1/)
	return QkF * (kF_invA*aB)					// kF_invA is 1/->kF, aB is -> au
End

Function Qang2au(qAngstroms)
	Variable qAngstroms						// q in 1/
	NVAR aB=root:Packages:Constants:aB
	return qAngstroms * aB						// convert to q in atomic units
End

Function Qau2ang(Qau)
	Variable Qau								// q in atomic units
	NVAR aB=root:Packages:Constants:aB
	return Qau / aB								// convert to q in 1/
End

Function hkl2au(hkl)
	Variable hkl								// magnitude of G in hkl units
	NVAR aB=root:Packages:Constants:aB
	NVAR ao=root:Packages:inelasticUtility:ao_A	// lattice const ()
	return  (2*PI * hkl/ao *aB)
End

Function hkl2ang(hkl)
	Variable hkl								// magnitude of G in hkl units
	NVAR ao=root:Packages:inelasticUtility:ao_A	// lattice const ()
	return  (2*PI * hkl/ao)
End




Function pQau2kF(Qau)
	Variable Qau								// q in atomic units
	if (Qau<0)
		Qau = 0
		Prompt Qau, "give Q in atomic units"
		DoPrompt "au -> kF", Qau
	endif
	printf "convert Q,    %g (au) == %g (kF)\r",Qau, Qau2kF(Qau)
End
Function pQkF2au(QkF)
	Variable QkF								// q in kFermi
	if (QkF<0)
		QkF = 0
		Prompt QkF, "give Q in units of k-Fermi"
		DoPrompt "kF -> au", QkF
	endif
	printf "convert Q,   %g (kF) == %g (au)\r",QkF, QkF2au(QkF)
End
Function pQang2au(qAngstroms)
	Variable qAngstroms						// q in 1/
	if (qAngstroms<0)
		qAngstroms = 0
		Prompt qAngstroms, "give Q in units of (1/)"
		DoPrompt "(1/) -> au", qAngstroms
	endif
	printf "convert Q,   %g (1/) == %g (au)\r",qAngstroms, Qang2au(qAngstroms)
End
Function pQau2ang(Qau)
	Variable Qau								// q in atomic units
	Variable qAngstroms						// q in 1/
	if (Qau<0)
		Qau = 0
		Prompt Qau, "give Q in units of atomic"
		DoPrompt "au -> (1/)", Qau
	endif
	printf "convert Q,   %g (au) == %g (1/)\r",Qau, Qau2ang(Qau)
End
Function phkl2au(hkl)
	Variable hkl								// magnitude of G in hkl units
	Prompt hkl, "give |hkl| in rlu"
	if (hkl<0)
		hkl = 0
		Prompt hkl, "give |hkl| in rlu"
		DoPrompt "|hkl| -> au", hkl
	endif
	printf "convert Q,   |hkl| = %g (rlu) == %g (au)\r",hkl, hkl2au(hkl)
End
Function phkl2ang(hkl)
	Variable hkl								// magnitude of G in hkl units
	Prompt hkl, "give |hkl| in rlu"
	if (hkl<0)
		hkl = 0
		Prompt hkl, "give |hkl| in rlu"
		DoPrompt "|hkl| -> ", hkl
	endif
	printf "convert Q,   |hkl| = %g (rlu) == %g (1/)\r",hkl, hkl2ang(hkl)
End
Function phkl2kF(hkl)
	Variable hkl								// magnitude of G in hkl units
	if (hkl<0)
		hkl = 0
		Prompt hkl, "give |hkl| in rlu"
		DoPrompt "hkl -> kF", hkl
	endif
	printf "convert Q,   |hkl| = %g (rlu) == %g (kF)\r",hkl, hkl2kF(hkl)
End




Function allowedReflections(Qmax,ao)
	Variable Qmax			// max Q (1/nm)
	Variable ao				// =4.177				// this needs to be changed to match the material
	String aoStr=""
	if (numtype(Qmax+ao) || Qmax<1 || ao<=0)
//		String latticeConstants="4.177  NiO;4.14  CoO;5.462  CaF2;4.05  Al;3.52  Ni;2.88  Cr;5.43102088  Si;5.657764  Ge;TiO2  a=4.59 b=4.59 c=2.96 tetragonal;Sc  a=3.3088  c=5.268, hex"
//		String latticeConstants="4.177  NiO;4.14  CoO;5.462  CaF2;4.05  Al;3.52  Ni;2.88  Cr;5.43102088  Si;5.657764  Ge;TiO2  a=4.59 b=4.59 c=2.96 tetragonal;Sc  a=3.3088  c=5.268, hex;4.02767  LiF"
//		String latticeConstants="4.177 NiO fcc;4.14 CoO fcc;5.462 CaF2 fcc;4.05 Al fcc;3.52 Ni fcc;2.88 Cr bcc;5.43102088 Si fcc;5.657764 Ge fcc;TiO2 a=4.59 b=4.59 c=2.96 tetragonal;Sc a=3.3088 c=5.268, hex;4.02767 LiF fcc"
//		String latticeConstants="4.177 NiO fcc;4.14 CoO fcc;5.462 CaF2 fcc;4.05 Al fcc;3.52 Ni fcc;2.88 Cr bcc;5.43102088 Si fcc;5.657764 Ge fcc;TiO2 a=4.59 b=4.59 c=2.96 tetragonal;Sc a=3.3088 c=5.268, hex;4.02767 LiF fcc;4.307 FeO fcc;4.4457 MnO fcc"
		String latticeConstants="4.177 NiO fcc;4.14 CoO fcc;5.462 CaF2 fcc;4.05 Al fcc;3.52 Ni fcc;2.88 Cr bcc;5.43102088 Si fcc;5.657764 Ge fcc;"
		latticeConstants += "TiO2 a=4.59 b=4.59 c=2.96 tetragonal;Sc a=3.3088 c=5.268, hex;4.02767 LiF fcc;4.307 FeO fcc;4.4457 MnO fcc;3.905 SrTiO3"
		Prompt Qmax, "max Q (1/nm)"
		Prompt aoStr, "lattice constant",popup,latticeConstants
		DoPrompt "max Q",Qmax,aoStr
		if (V_flag)
			return NaN
		endif
		ao = str2num(aoStr)
	endif
	if (numtype(Qmax+ao) || Qmax<1 || ao<=0)
		printf "Unable to proceed with Qmax = %g  and ao=%g\r",Qmax,ao
		return NaN
	endif
	if (strlen(aoStr)<1)
		sprintf aoStr,"ao=%g",ao
	endif

	String lattice = StringFromList(2,aoStr," ")
	if (stringmatch(lattice,"FCC"))
		FUNCREF allowAll f=allowFCC
	elseif (stringmatch(lattice,"BCC"))
		FUNCREF allowAll f=allowBCC
	else
		FUNCREF allowAll f=allowAll
	endif

	Variable maxhkl = ceil(ao*Qmax/10/2/PI)	// note, here Qmax(1/nm), ao()
	Make/N=3/O hkl
	Make/N=1000/O hh,kk,ll,QQ
	Variable i,h,k,l,Q,dspace
	i = 0
	for (h=0;h<=maxhkl;h+=1)
		hkl[0] = h
		for (k=h;k<=maxhkl;k+=1)
			hkl[1]=k
			for (l=k;l<=maxhkl;l+=1)
				hkl[2] = l
					dspace = ao/norm(hkl)
					Q = 2*PI/dspace*10		// Q (1/nm)
					if (f(hkl) && Q<Qmax)
						// printf "(%d %d %d),   Q=%.3f\r",h,k,l,Q
						hh[i] = h
						kk[i] = k
						ll[i] = l
						QQ[i] =Q
						i += 1
					endif
			endfor
		endfor
	endfor
	Variable N=i
	Redimension/N=(i) hh,kk,ll,QQ
	Sort QQ, hh,kk,ll,QQ
	printf "\rfor   %s\r\r",aoStr
//	printf "(h k l),    (1/)\r"
	printf "(h k l)\r"
	for (i=0;i<N;i+=1)
		printf "(%d %d %d),    Q = %.3f(1/)  =  Q=%.2f(1/nm)\r",hh[i],kk[i],ll[i],QQ[i]/10,QQ[i]
	endfor
	printf "to a maximum Q of %g (1/nm)\r",Qmax
	KillWaves/Z hkl,hh,kk,ll,QQ
	return N
End
Function allowAll(H)		// everybody is allowed
	Wave H
	return 1
End
Static 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
Static Function allowBCC(H)	// flags whether H=(hkl) is an allowed BCC reflection (sum is even)
	Wave H
	if ( (abs(H[0])+abs(H[1])+abs(H[2]))==0 )	// the (000) is not allowed
		return 0
	endif
	return !mod(sum(H,-inf,inf),2)
End


Function/S OneReflections(hklStr,E)
	String hklStr
	Variable E					// Energy (keV)

	Variable ao=4.6837, bo=3.4226, co=5.1288
	Variable alpha=90, beta=99.54, gam=90

	if (numtype(E) || E < 2)
		E = 	(numtype(E) || E < 2) ? 7.59416 : E
		Prompt hklStr,  "string with the hkl (comma separated)"
		Prompt E, "energy (keV)"
		DoPrompt "hkl and E",hklStr,E
		if (V_flag)
			return ""
		endif
	endif
	hklStr = ReplaceString(",", hklStr,";" )
	hklStr = ReplaceString(" ", hklStr,";" )
	Make/N=3/O/D hkl
	hkl[0] = str2num(StringFromList(0,hklStr))
	hkl[1] = str2num(StringFromList(1,hklStr))
	hkl[2] = str2num(StringFromList(2,hklStr))
	if (numtype(E+sum(hkl)) || E < 2)
		printf "Unable to proceed with hkl = %s,  E = %g\r",hklStr, E
		return ""
	endif

	// Make the real lattice vectors
	Variable ca=cos(alpha*PI/180), cb = cos(beta*PI/180), cg = cos(gam*PI/180)
	Variable sg=sin(gam*PI/180)

	Make/N=3/O/D a1,a2,a3, a3hat
	a1 = {ao,0,0}
	a2 = { bo*cg, bo*sg, 0}

	a3hat[0] = cb
	a3hat[1] = ( ca - cg*cb ) / sg
	Variable z2 = 1-a3hat[0]^2-a3hat[1]^2
	a3hat[2] = (z2<0) ? 0 : sqrt(z2)
	a3 = a3hat*co

	// make the reciprocal lattice vectors
	Make/N=3/O/D r1,r2,r3
	Make/N=(3,3)/O/D recip				// make the reciprocal lattice
	Variable Vc
	Cross a2,a3
	Wave W_Cross = W_Cross
	recip[][0] = W_Cross[p]	// r1 = (a2 x a3) / Vc
	Vc = MatrixDot(a1,W_Cross)
	Cross a3,a1
	recip[][1] = W_Cross[p]	// r2 = (a3 x a1) / Vc
	Cross a1,a2
	recip[][2] = W_Cross[p]	// r3 = (a1 x a2) / Vc
	recip *= 2*PI/Vc

	MatrixMultiply recip,hkl
	Wave M_product=M_product
	Variable Q = norm(M_product)

	NVAR hc = root:Packages:Constants:hc
	Variable tth, lambda = hc/E
	tth = 2*asin(Q*lambda/(4*PI)) * 180/PI

	String out,str
	sprintf out, "(% d % d % d)@%.4f keV; \tQ=%.3f (1/),\td=%.4f ,",hkl[0],hkl[1],hkl[2],E,Q,2*PI/Q
	if (!numtype(tth) && tth > 0.1)
		sprintf str,"\t\t2th=%.3fĄ,",tth
		out += str
	endif
	sprintf str,"\t\t a=%g , b=%g , c=%b ,   alpha=%gĄ, beta=%gĄ, gamma=%gĄ",ao,bo,co,alpha,beta,gam
	out += str
	if (ItemsInList(GetRTStackInfo(0))<2)
		printf "%s\r", out
	endif
	KillWaves/Z hkl, a1,a2,a3, a3hat, W_Cross, r1,r2,r3, recip, M_product
	return out
End

Function AllReflections(Qmin,Qmax)
	Variable Qmin,Qmax

	Variable ao=4.6837, bo=3.4226, co=5.1288
	Variable alpha=90, beta=99.54, gam=90

	if (numtype(Qmin+Qmax) || Qmax<1 || Qmin<0)
		Prompt Qmin, "min Q (1/)"
		Prompt Qmax, "max Q (1/)"
		DoPrompt "Q range",Qmin,Qmax
		if (V_flag)
			return NaN
		endif
	endif
	if (numtype(Qmin+Qmax) || Qmax<1 || Qmin<0)
		printf "Unable to proceed with Qmax = %g,  Qmax = %g\r",Qmin,Qmax
		return NaN
	endif
	Qmin = max(Qmin,1e-3)

	// Make the real lattice vectors
	Variable ca=cos(alpha*PI/180), cb = cos(beta*PI/180), cg = cos(gam*PI/180)
	Variable sg=sin(gam*PI/180)

	Make/N=3/O/D a1,a2,a3, a3hat
	a1 = {ao,0,0}
	a2 = { bo*cg, bo*sg, 0}

	a3hat[0] = cb
	a3hat[1] = ( ca - cg*cb ) / sg
	Variable z2 = 1-a3hat[0]^2-a3hat[1]^2
	a3hat[2] = (z2<0) ? 0 : sqrt(z2)
	a3 = a3hat*co

	// make the reciprocal lattice vectors
	Make/N=3/O/D r1,r2,r3
	Make/N=(3,3)/O/D recip				// make the reciprocal lattice
	Variable Vc
	Cross a2,a3
	Wave W_Cross = W_Cross
	recip[][0] = W_Cross[p]	// r1 = (a2 x a3) / Vc
	Vc = MatrixDot(a1,W_Cross)
	Cross a3,a1
	recip[][1] = W_Cross[p]	// r2 = (a3 x a1) / Vc
	Cross a1,a2
	recip[][2] = W_Cross[p]	// r3 = (a1 x a2) / Vc
	recip *= 2*PI/Vc

	Make/N=3/O/D hkl
	Variable amax = max(max(ao,bo),co)
	Variable maxhkl = ceil(amax*Qmax/2/PI)
	Variable Nmax = (2*maxhkl+1)^3

	printf "maxhkl = %d,   Nmax = %d\r",maxhkl,  Nmax
	Make/N=(Nmax)/O hh,kk,ll,QQ
	Variable i,h,k,l,Q
	i = 0
	for (l=maxhkl;l>=-maxhkl;l-=1)
		hkl[2] = l
		for (k=maxhkl;k>=-maxhkl;k-=1)
			hkl[1]=k
			for (h=maxhkl;h>=-maxhkl;h-=1)
				hkl[0] = h
				MatrixMultiply recip,hkl
				Wave M_product=M_product
				Q = norm(M_product)
				if (Q<=Qmax && Q>=Qmin)
					// printf "(%d %d %d),   Q=%.3f\r",h,k,l,Q
					hh[i] = h
					kk[i] =  k
					ll[i] = l
					QQ[i] =Q
					i += 1
				endif
			endfor
		endfor
	endfor

	Variable N=i
	Redimension/N=(i) hh,kk,ll,QQ
	Sort QQ, hh,kk,ll,QQ

	printf "\rfor a=%g , b=%g , c=%b ,   alpha=%gĄ, beta=%gĄ, gamma=%gĄ\r\r",ao,bo,co,alpha,beta,gam
	printf "(h k l),    (1/)\r"
	Variable lastQ=-1, tth
	NVAR hc = root:Packages:Constants:hc
	Variable lambda = hc/7.59416
	for (i=0;i<N;i+=1)
		if (QQ[i]-lastQ > .0001)
			tth = 2*asin(QQ[i]*lambda/(4*PI)) * 180/PI
			if (numtype(tth) || tth < 0.1)
				printf "(% d % d % d), \t\tQ=%.3f,\td=%.4f\r",hh[i],kk[i],ll[i],QQ[i],2*PI/QQ[i]
			else
				printf "(% d % d % d), \t\tQ=%.3f,\td=%.4f\t2th=%.2fĄ\r",hh[i],kk[i],ll[i],QQ[i],2*PI/QQ[i],tth
			endif
			lastQ = QQ[i]
		endif
	endfor
	printf "to a maximum Q of %g (1/)\r",Qmax
	KillWaves/Z hkl,hh,kk,ll,QQ
	KillWaves/Z a1,a2,a3, a3hat, W_Cross, r1,r2,r3, recip, M_product
	return N
End






Function inelasticInitPackage()
	Execute/P "INSERTINCLUDE  \"PhysicalConstants\""		// this package needs the PhysicalConstants package
	Execute/P "COMPILEPROCEDURES "
	Execute/P "PhysicalConstantsInitPackage()"
	Execute/P "DELETEINCLUDE  \"PhysicalConstants\""
	Execute/P "COMPILEPROCEDURES "

	NewDataFolder /O root:Packages
	NewDataFolder /O root:Packages:inelasticUtility
	if (exists("root:Packages:inelasticUtility:ao_A")!=2)
		// default lattice constant
		Variable /G root:Packages:inelasticUtility:ao_A=4.048		// lattice const. of Al ()
	endif
	if (exists("root:Packages:inelasticUtility:kF_invA")!=2)
		// length of Fermi wave-vector in 1/
		Variable /G root:Packages:inelasticUtility:kF_invA=1.75		// kFermi of Al (1/)
	endif
End