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

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


Menu "inelastic"
	"Q of scan", Q_of_scan(0)
	"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)
End


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



Function f_sumRule(Q)					// f-sum rule = n * � * (e^2/aB) * (Q*aB)^2s   (eV/�^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)
		printf "f-sum(Q=%g (1/�), %g (�^3/electron)) = %g (eV/�^3)", Q, Ve, f_sum
		if (not_aluminum)
			printf "        this is NOT Aluminum"
		endif
		print ""
	endif
	return f_sum
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
	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)
		scanNum=NumVarOrDefault("root: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
	if (numtype(theta))
		theta = refspecInfo(scanNum,"Two-theta")/2
	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)		// was DCM_energy, instead of Eo
	if (numtype(Q))
		Q = refspecInfo(scanNum,"Q_angstrom")
	endif
	if (doPrint)
		printf "Q(scan %d) = %g (1/�)\r", scanNum, Q
	endif
	return Q
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 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