reset

# define energy resolution
eps = 0.01 #e-01

# define Debye-Waller factor (and set it to 1 for simplicity
G(x) = 1.0

# phonon dispersion
# v_s = 1.0 # speed of sound waves
Omega(k) = 2.0 * abs(sin(k*pi/2)) #v_s*abs(k)

# define gaussian to model delta function with finite energy resolution eps
g(x) = exp(-0.5 * x**2 / (eps**2)) / sqrt(2*pi*eps**2)

# define Bose-Einstein distribution function
n(T,k) = 1.0/(exp(Omega(k)/T)-1)

# define dos 1-Phonon contribution 
dos1(T,w,k) = n(T,k) * g(w + Omega(k)) + (1.0+n(T,k))*g(w - Omega(k))

set terminal qt enhanced font "Verdana, 18"
set pm3d
set isosamples 100,100
set hidden3d
set xyplane at 0
set view 0,0

# define temperature (in units of hbar/k_B)
T=0.5

q_min = -0.95
q_max = 0.95
w_min = -1.0*1.2*Omega(q_min)
w_max = 1.2*Omega(q_max)

set xlabel "wave vector q [{/Symbol p}/a]"
set ylabel "frequency {/Symbol w}" rotate parallel offset -5.,0.
set cblabel "S(q,{/Symbol w})"

splot [q_min:q_max][w_min:w_max] (x**2*G(x)/Omega(x))*pi**3*dos1(T,y,x) t sprintf("1-phonon DSF S_1(q,{/Symbol w}), {/Symbol e}=%g", eps)