#pragma rtGlobals=1		// Use modern global access method.
#include "JelliumResponse"
  
// Variable/C Gconst = cmplx(1.,0.)
// Variable/C :root:one = cmplx(1.,0.)

Function Sqw_const(hkl,qkF,Gc,eV)
	Variable hkl	// hkl = 1, 11, 13
	Variable qkF
	Variable/C Gc
	Variable eV

	if (exists("root:kernels:kernel") == 1)
		WAVE kernel=root:kernels:kernel
	else
		if (exists("kernel") != 1)
			print "Making a local copy of the kernel"
			Make kernel = {0.0339806, 0.0776699, 0.194175, 0.38835, 0.194175, 0.0776699, 0.0339806}
			SetScale/I x -1,1,"eV", kernel
		endif
		WAVE kernel=kernel
	endif

	Variable yy
	Variable i
	Variable Nkernel=numpnts(kernel)
	NVAR kF=root:Packages:inelasticUtility:kF_invA

	i = 0
	yy = 0
	do
		yy += S_of_G(hkl,-eV-pnt2x(kernel,i),qkF*kF,Gc) * kernel[i]
		i += 1
	while (i<Nkernel)
	return yy
End


Function Sqw_jel_const(qkF,Gc,eV)
	Variable qkF
	Variable/C Gc
	Variable eV

	if (exists("root:kernels:kernel") == 1)
		WAVE kernel=root:kernels:kernel
	else
		if (exists("kernel") != 1)
			print "Making a local copy of the kernel"
			Make kernel = {0.0339806, 0.0776699, 0.194175, 0.38835, 0.194175, 0.0776699, 0.0339806}
			SetScale/I x -1,1,"eV", kernel
		endif
		WAVE kernel=kernel
	endif

	Variable yy
	Variable i
	Variable Nkernel=numpnts(kernel)
	NVAR kF=root:Packages:inelasticUtility:kF_invA

	i = 0
	yy = 0
	do
		yy += SjelliumG(-eV-pnt2x(kernel,i),qkF*kF,Gc) * kernel[i]
		i += 1
	while (i<Nkernel)
	return yy
End




Function S_of_G(hkl,eV,q,G)		// this gives X as a function of Xo, v, and G with q in 1/�
	Variable hkl		// hkl of direction, used to pick X0_013 or X0_001
	Variable q			// magnitude of q  (1/�)
	Variable eV			// energy loss, is normally negative
	Variable/C G		// complex value of G to use

	Variable/C xx
	return -imag(X_of_G(hkl,eV,q,G))
End


Function/C X_of_G(hkl,eV,q,G)		// this gives X as a function of Xo, v, and G with q in 1/�
	Variable hkl		// hkl of direction, used to pick X0_013 or X0_001
	Variable q			// magnitude of q  (1/�)
	Variable eV			// energy loss, is normally negative
	Variable/C G		// complex value of G to use
	Variable/C one = cmplx(1.,0.)
	WAVE X0_013=root:AL_X0:X0_013
	WAVE X0_001=root:AL_X0:X0_001
	Variable/C xx
	Variable iq,iqlo,iqhi
	ev = -abs(eV)
	if (hkl==013)
		iq = (q-DimOffset(X0_013, 1)) / DimDelta(X0_013, 1)
		iqlo = floor(iq)
		iqhi = ceil(iq)
		if (iqhi==iqlo)
			xx = X0_013(eV)[iq]
		else
			 xx = X0_013(eV)[iqlo]*(iqhi-iq)  +   X0_013(eV)[iqhi]*(iq-iqlo)
		endif
		return (xx / ( one - vq_eV(q)*( one-G)*xx))
	endif
	if (hkl==001)
		iq = (q-DimOffset(X0_001, 1)) / DimDelta(X0_001, 1)
		iqlo = floor(iq)
		iqhi = ceil(iq)
		if (iqhi==iqlo)
			xx = X0_001(eV)[iq]
		else
			 xx = X0_001(eV)[iqlo]*(iqhi-iq)  +   X0_001(eV)[iqhi]*(iq-iqlo)
		endif
		return (xx / ( one - vq_eV(q)*( one-G)*xx))
	endif
	if (hkl==011)
		NVAR kF=root:Packages:inelasticUtility:kF_invA
		if ((abs(q-(1.45*kF)) < 0.47) %& (abs(q-(1.33*kF))>.001) %& (abs(q-(1.5*kF))>.001))
//		if ((abs(q-(1.45*kF)) < 0.47) %& (abs(q-(1.33*kF))>.001))
			Abort "uncomment this if you want this feature"
			Variable fraclo = (2.62 - q)/(2.62-2.33)
			Wave X233=root:Al_X0:X0_011_233
			Wave X262=root:Al_X0:X0_011_262
			xx = X233(eV) * fraclo + X262(eV) * (1-fraclo)
			return (xx / ( one - vq_eV(q)*( one-G)*xx))
		endif
		Variable ii = round(q*100-.1)
		if (ii==124)
			ii = 125		
		endif
		Wave X0 = $("root:Al_X0:X0_011_"+num2istr(ii))
		xx = X0(eV)
		return (xx / ( one - vq_eV(q)*( one-G)*xx))
	endif
	Abort ("illegal value of hkl entered into X_of_G, hkl must be 001 or 013")
End



Function vq_eV(q)					// this gives v in 1/eV, using q in 1/�
	Variable q
	//  8.1711*4*PI = 102.681
	return (102.681/(q*q))
EndMacro



//Macro MakeS0surface(hkl)
//	Variable hkl
//	Silent 1
//	if (hkl==013)
//		Make/N=(501,22)/O S0_013
//		S0_013 = - imag(X0_013)
//		SetScale/I x 0,-100,"eV", S0_013
//		SetScale/I y 2.5,4.6,"1/�", S0_013
//		Graph_S0_013()
//		return
//	endif
//	if (hkl==001)
//		Make/N=(501,28)/O S0_001
//		SetScale/I x 0,-100,"eV", S0_001
//		SetScale/I y 1.2,3.9,"1/�", S0_001
//		S0_001 = - imag(root:AL_X0:X0_001)
//		Graph_S0_001()
//		return
//	endif
//
//	Abort ("illegal value of hkl entered into X_of_G, hkl must be 001 or 013")
//EndMacro
//
//
//Window Graph_S0_013()
//	PauseUpdate; Silent 1		// building window...
//	Display /W=(5,42,600,564)
//	AppendMatrixContour S0_013
//	ModifyGraph gfSize=18
//	ModifyGraph tick=2
//	ModifyGraph mirror=1
//	ModifyGraph minor=1
//	ModifyGraph lblMargin(left)=9,lblMargin(bottom)=7
//	ModifyGraph standoff=0
//	Label left "Q  (\\U)"
//	Label bottom "Energy Loss  (\\U)"
//	SetAxis bottom 0,-100
//EndMacro
//
//Window Graph_S0_001()
//	PauseUpdate; Silent 1		// building window...
//	Display /W=(351,78,932,574)
//	AppendMatrixContour S0_001
//	ModifyGraph gfSize=18
//	ModifyGraph tick=2
//	ModifyGraph mirror=1
//	ModifyGraph minor=1
//	ModifyGraph lblMargin(left)=5
//	ModifyGraph standoff=0
//	ModifyGraph axOffset(left)=-1.1,axOffset(bottom)=-0.388889
//	Label left "Q  (\\U)"
//	Label bottom "Energy Loss  (\\U)"
//	SetAxis bottom 0,-100
//EndMacro