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

Proc doNotUse()
	x = thetaM / T
	Phi(x) = 1/x * integral(0,x,t/(exp(t)-1)  dt)
	2M = (3 * Q^2 * hbar^2 * T) / (mass * kB * thetaM^2) *  { Phi(x) + x/4 }

	Note that mass is in mass, not energy.  E = m c^2

	//	make/N=501 phi
	//	SetScale/I x 0,20,"", phi
	//	phi = compute_one_phi(x)
	//	KillWaves  phi
EndMacro

Function Phi(xx)		// the Phi function
	Variable xx
	if (!exists("phi_temp_"))
		Make/N=501/O phi_temp_
	else
		WAVE phi_temp_=phi_temp_
	endif
	SetScale/I x 0,(xx),"", phi_temp_

	if (xx<=0)
		return 1
	endif
	phi_temp_ = x/(exp(x)-1)
	phi_temp_[0]=1
	return (area(phi_temp_, 0, xx)/xx)
End

Function DW_factor_M(T,thetaM,Q,amu)
	Variable T			// Temperature (K)
	Variable thetaM		// Debye Temperature (K)
	Variable Q			// length of q vector (1/�)
	Variable amu		// mass of atom (amu)

	Variable amu_eV = 931.49432e6	// energy of one amu (eV)
	Variable hbar = 6.582122e-16  	// eV - sec
	Variable kB = 8.617385e-5		// eV / K
	Variable c = 2.99792458e8		// speed of light (m/sec)

	if (T<=0)
		return 0
	endif

	Variable xx = thetaM / T
	Variable M,B

	B = (3*hbar^2 * T * c^2)/(2*(amu*amu_eV)*kB*thetaM^2)* (Phi(xx) + xx/4)
	M = B *(Q*1.e10)^2		// Q is in 1/�, but we need it in 1/m
	return M
End


Macro Get_Debye_Waller_Panel()
	Silent 1
	String name = WinList("Debye_Waller_Panel", "", "WIN:64")

	if (exists("temperature")!=2)
		Variable/G temperature = 295
	endif
	if (exists("Qang")!=2)
		Variable/G Qang = 2.5
	endif
	if (exists("DWM")!=2)
		Variable/G DWM = 0
	endif
	if ( exists("elementNo")!=2)
		Variable/G elementNo = 1
	endif
	if ( NotWave("root:thetaMs") || NotWave("root:amus") || NotWave("root:eNames"))
		Make/O root:thetaMs={394,315,625,360,275}
		Make/O root:amus={26.98,63.54,28.086,72.59,106.4}
		Make/O/T root:eNames={"Al","Cu","Si","Ge","Pd"}
	endif
	if (strlen(name)>1)
		DoWindow /F $name
	else
		Debye_Waller_Panel()		// create it since it does not exist
	endif
EndMacro
Function NotWave(name)
	String name
	Variable i = (exists(name)!=1)
	return i
End


Macro Debye_Waller(element,T,Q)
	Variable element=0
	Prompt element, "element", popup, "Al;Cu;Si;Ge;Pd"
	Variable T=300				// Temperature (K)
	Prompt T, "Temperature (K)"
	Variable Q=2.975			// length of Q vector (1/�)
	Prompt Q, "Length of Q vector (1/�)"

	Silent 1
	element -= 1
	Make/O root:thetaMs={394,315,625,360,275}
	Make/O root:amus={26.98,63.54,28.086,72.59,106.4}
	Make/O/T root:eNames={"Al","Cu","Si","Ge","Pd"}

	Variable thetaM=root:thetaMs[element]	// Debye Temperature (K)
	Variable amu=root:amus[element]		// mass of atom (amu)

	print "element =",root:eNames[element],"      Debye Temperature =",thetaM, "     ",amu,"amu"

	Variable xx = thetaM / T
	Variable M = DW_factor_M(T,thetaM,Q,amu)

	print "ThetaM / T = ",thetaM / T,"           M =",M,"         exp(-2M) =",exp(-2*M)
	KillWaves/Z phi_temp_
End


Function PopMenuProc(ctrlName,popNum,popStr) : PopupMenuControl
	String ctrlName
	Variable popNum
	String popStr
	NVAR elementNo = root:elementNo
	elementNo=popNum
	ButtonProc(ctrlName)
End

Function ButtonProc(ctrlName) : ButtonControl
	// set exp(-M) and M based on other things Qang, temperature, elementNo
	String ctrlName

	NVAR M = root:DWM
	NVAR elementNo=root:elementNo
	NVAR Q=root:Qang
	NVAR T = root:temperature

	WAVE thetaMs=root:thetaMs
	WAVE amus = root:amus
	WAVE/T eNames = root:eNames

	Variable thetaM=thetaMs[elementNo-1]		// Debye Temperature (K)
	Variable amu=amus[elementNo-1]			// mass of atom (amu)
	Variable xx = thetaM / T

	M = DW_factor_M(T,thetaM,Q,amu)
	KillWaves/Z phi_temp_
End

Window Debye_Waller_Panel() : Panel
	PauseUpdate; Silent 1		// building window...
	NewPanel /W=(7,68,307,237)
	PopupMenu elementName,pos={0,2},size={97,19},proc=PopMenuProc,title="element"
	PopupMenu elementName,mode=1,popvalue="Al",value= #"\"Al;Cu;Si;Ge;Pd\""
	SetVariable temperature,pos={8,30},size={230,20},title="Temperature (K�)"
	SetVariable temperature,fSize=18,limits={0,INF,1},value= temperature
	SetVariable Q,pos={130,4},size={107,20},title="Q",fSize=18
	SetVariable Q,limits={0,INF,0.1},value= QAng
	Button button0,pos={115,63},size={50,30},proc=ButtonProc,title="Go"
	ValDisplay val_DebyeTemp,pos={11,102},size={253,20},title="Debye Temperature ="
	ValDisplay val_DebyeTemp,fSize=18,format="%g� (K)",frame=0
	ValDisplay val_DebyeTemp,limits={0,0,0},barmisc={0,1000}
	ValDisplay val_DebyeTemp,value= #"thetaMs[elementNo-1]"
	ValDisplay valdisp_M,pos={17,134},size={95,20},title="M =",fSize=18
	ValDisplay valdisp_M,format="%.5f",frame=0,limits={0,0,0},barmisc={0,1000}
	ValDisplay valdisp_M,value= #"root:DWM"
	ValDisplay valdisp_exp_M,pos={143,134},size={147,20},title="exp(-M) =",fSize=18
	ValDisplay valdisp_exp_M,format="%.4f",frame=0,limits={0,0,0},barmisc={0,1000}
	ValDisplay valdisp_exp_M,value= #"exp(-DWM)"
EndMacro