reset

# define energy resolution
eps = 0.1 #e-01

# define poisson distribution
p(n, l) = exp(-l) * (l**n) / gamma(n+1)

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

# define dos for different orders of perturbation theory
dos0(n, w, w_qho) = g(n, w, w_qho)
dos1(n, w, w_qho) = 0.5 * (g(n+1, w, w_qho) + g(n-1, w, w_qho))

# we want a sum which we define recursively until n_max, dos is a function evaluating the density of states as function of frequency
sum0(n, l, w, w_qho) = (n<n_max)? p(n,l)*dos0(n, w, w_qho)+sum0(n+1, l, w, w_qho): p(n, l)*dos0(n, w, w_qho)
sum1(n, l, w, w_qho) = (n<n_max)? p(n,l)*dos1(n, w, w_qho)+sum1(n+1, l, w, w_qho): p(n, l)*dos1(n, w, w_qho)

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

# define oscillator frequency
w_qho = 1.0

# define x0 and deltax
x0 = 1.0
deltax = 1.0

# define iteration depth
n_max = 10

q_min = -1.0
q_max = 1.0
w_min = -(n_max/2-0.5)*w_qho
w_max = (n_max+0.5)*w_qho

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

splot [2*q_min:2*q_max][w_min:w_max] sum0(0, (x*x0)**2, y, w_qho) t sprintf("S(q,{/Symbol w}) up to 0th order in (q{/Symbol d}x)^2, {/Symbol e}=%g", eps)
#splot [5*q_min:5*q_max][w_min:w_max] sum0(0, (x*x0)**2, y, w_qho) + (x*deltax)**2*sum1(0, (x*x0)**2, y, w_qho) t sprintf("S(q,{/Symbol w}) up to 1st order in (q{/Symbol d}x)^2, {/Symbol e}=%g", eps)