clear all
% Phase space evolution of the expectation 
% values of a quantum damped harmonic oscillator

tau = linspace(0,100,1000);

% parameters
omega = 1; % natural frequency
gamma = 0.1; % damping strength

q_0 = 0; % initial expectation value of q in dimensionless units
p_0 = 10; % initial expectation value of p in dimensionless units

zeta = sqrt(omega^2 - gamma^2/4);


g_2 = omega/zeta*sin(zeta/omega*tau).*exp(-1/2*gamma/omega*tau);
dg_2 = (cos(zeta/omega*tau)-1/2*gamma/zeta*sin(zeta/omega*tau)).*exp(-1/2*gamma/omega*tau);
ddg_2 = (-gamma/omega*cos(zeta/omega*tau) + (gamma^2/2 - omega^2)/omega^2*omega/zeta*sin(zeta/omega*tau)).*exp(-1/2*gamma/omega*tau);
dg_2n = diff(g_2,1)./diff(tau);
ddg_2n = diff(g_2,2)./(diff(tau(1:end-1)).^2);

figure(1)
subplot(1,2,1)
plot(tau,g_2,'b')
hold on
plot(tau(1:end-1),dg_2n,'g.')
plot(tau,dg_2,'g')
plot(tau(1:end-2),ddg_2n,'r.')
plot(tau,ddg_2,'r')
title('G_2 and derivatives','Fontsize',18)
xlabel('\omega_0 t','Fontsize',18)
axis square
hold off


q = dg_2*q_0 + g_2*p_0;
p = ddg_2*q_0 + dg_2*p_0;


figure(1)
subplot(1,2,2)
plot(q,p,'k')
hold on
plot(q_0,p_0,'r.')
axis square
title('Evolution in phase space','Fontsize',18)
xlabel('Q/Q_0','Fontsize',18)
ylabel('P/P_0','Fontsize',18)

hold off


