clear all
% Time evolution of a closed two level system 


% parameters of the Hamiltonian

D = 0 % splitting of the levels
b = -1 % hopping between the levels

%discretization parameters
dt = 0.01; % step
Nt = 500; % number of steps
t = [0:Nt-1]*dt; % time vector


% The Hamiltonian
H = [D/2 b; b -D/2];

[U,Delta] = eig(H);


%initial condition (in position space)
% rho_Pos = [1 0; 0 0]; %left initialization
% rho_Pos = [0 0; 0 1]; %right initialization
% rho_Pos = [1 1; 1 1]; %even initialization
% rho_Pos = [1 -1; -1 1]; %odd initialization
% rho_Pos = rand(2) + i*rand(2); %random initialization

% this procedure takes care that the 
% input for the initial condition 
% is Hermitian and normalized
rho_Pos = (rho_Pos + rho_Pos')/(2*trace(rho_Pos));

%initial condition (in energy space)
rho_En = U'*rho_Pos*U;


% time evolution both in position and energy basis
for n = 1:Nt
    rhoPos_vec(:,:,n) = rho_Pos;
    rhoEn_vec(:,:,n) = rho_En;
    rho_Pos = rho_Pos -i*(H*rho_Pos - rho_Pos*H)*dt;
    rho_En = U'*rho_Pos*U;
end

% extraction of the relevant quantities
P1 = squeeze(rhoPos_vec(1,1,:));
P2 = squeeze(rhoPos_vec(2,2,:));
C12 = squeeze(rhoPos_vec(1,2,:));

Pe = squeeze(rhoEn_vec(1,1,:));
Po = squeeze(rhoEn_vec(2,2,:));
Ceo = squeeze(rhoEn_vec(1,2,:));




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                   %
%     GRAPHICAL REPRESENTATION      %
%                                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subplot(2,2,1)
plot(t,P1,'.-')
hold on
plot(t,P2,'ro-')
title('Populations (Pos)')
xlabel('t')
hold off

subplot(2,2,2)
plot(t,real(C12),'.-')
hold on
plot(t,imag(C12),'ro-')
title('Coherences (Pos)')
xlabel('t')
hold off

subplot(2,2,3)
plot(t,Pe(1:end),'.-')
hold on
plot(t,Po(1:end),'ro-')
title('Populations (En)')
xlabel('t')
hold off

subplot(2,2,4)
plot(t,real(Ceo),'.-')
hold on
plot(t,imag(Ceo),'ro-')
title('Coherences (En)')
xlabel('t')
hold off


hold off



