clear all
% Time evolution of a closed three level system 


% parameters of the Hamiltonian

En = [1,1,1]; % On-site energy in position
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 = [En(1) b b; b En(2) b; b b En(3)];

[U,Delta] = eig(H);


%initial condition (in position space)
rho_Pos = [1 0 0 ; 0 0 0; 0 0 0]; %left initialization
% rho_Pos = [0 0 0 ; 0 1 0; 0 0 0]; %center initialization
% rho_Pos = [0 0 0; 0 0 0; 0 0 1]; %right initialization
% rho_Pos = rand(3) + i*rand(3); % 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;

%important variables for the analytical solution

rho_En0 = U'*rho_Pos*U;
Bohr = repmat(diag(Delta),1,3) - repmat(diag(Delta)',3,1); 


% time evolution both in position and energy basis
for n = 1:Nt
    %numerical
    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;
    %"analytical"
    rhoEn_an_vec(:,:,n) = rho_En0.*exp(-i*Bohr*t(n));
    rhoPos_an_vec(:,:,n) = U*rhoEn_an_vec(:,:,n)*U';
    
end

% extraction of the relevant quantities

% Numerical

% populations in position
P1 = squeeze(rhoPos_vec(1,1,:));
P2 = squeeze(rhoPos_vec(2,2,:));
P3 = squeeze(rhoPos_vec(3,3,:));

% coherences in position
C12 = squeeze(rhoPos_vec(1,2,:));
C13 = squeeze(rhoPos_vec(1,3,:));
C23 = squeeze(rhoPos_vec(2,3,:));

% populations in energy
Pgr = squeeze(rhoEn_vec(1,1,:));
Pex1 = squeeze(rhoEn_vec(2,2,:));
Pex2 = squeeze(rhoEn_vec(3,3,:));

% coherences in energy
C_gr_ex1 = squeeze(rhoEn_vec(1,2,:));
C_gr_ex2 = squeeze(rhoEn_vec(1,3,:));
C_ex1_ex2 = squeeze(rhoEn_vec(2,3,:));

% "Analytical"

% populations in position
P1_an = squeeze(rhoPos_an_vec(1,1,:));
P2_an = squeeze(rhoPos_an_vec(2,2,:));
P3_an = squeeze(rhoPos_an_vec(3,3,:));

% coherences in position
C12_an = squeeze(rhoPos_an_vec(1,2,:));
C13_an = squeeze(rhoPos_an_vec(1,3,:));
C23_an = squeeze(rhoPos_an_vec(2,3,:));

% populations in energy
Pgr_an = squeeze(rhoEn_an_vec(1,1,:));
Pex1_an = squeeze(rhoEn_an_vec(2,2,:));
Pex2_an = squeeze(rhoEn_an_vec(3,3,:));

% coherences in energy
C_gr_ex1_an = squeeze(rhoEn_an_vec(1,2,:));
C_gr_ex2_an = squeeze(rhoEn_an_vec(1,3,:));
C_ex1_ex2_an = squeeze(rhoEn_an_vec(2,3,:));









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

subplot(2,2,1)
plot(t,[P1,P2,P3],'.-')
hold on
plot(t,[P1_an,P2_an,P3_an],'k')
title('Populations (Pos)')
xlabel('t')
hold off

subplot(2,2,2)
plot(t,[real(C12),imag(C12)],'.-')
hold on
plot(t,[real(C13),imag(C13)],'.-')
plot(t,[real(C23),imag(C23)],'.-')
hold on
plot(t,[real(C12_an),imag(C12_an)],'k')
title('Coherences (Pos)')
xlabel('t')
hold off

subplot(2,2,3)
plot(t,[Pgr,Pex1,Pex2],'.-')
hold on
plot(t,[Pgr_an,Pex1_an,Pex2_an],'k')
title('Populations (En)')
xlabel('t')
hold off

subplot(2,2,4)
plot(t,[real(C_gr_ex1),imag(C_gr_ex1)],'.-')
hold on
plot(t,[real(C_gr_ex2),imag(C_gr_ex2)],'.-')
plot(t,[real(C_ex1_ex2),imag(C_ex1_ex2)],'.-')

plot(t,[real(C_gr_ex1_an),imag(C_gr_ex1_an)],'k')
plot(t,[real(C_gr_ex2_an),imag(C_gr_ex2_an)],'k')
plot(t,[real(C_ex1_ex2_an),imag(C_ex1_ex2_an)],'k')
title('Coherences (En)')
xlabel('t')
hold off


hold off



