clear all
% Time evolution of a closed N level system  

N = 10;
% parameters of the Hamiltonian

En = ones(1,N); % On-site energy in position
b = -1 % hopping between the levels

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



% The Hamiltonian
% The Hamiltonian
H = diag(En) +...
    diag(b*ones(1,N-1),1) + diag(b*ones(1,N-1),-1) +...
    diag(b*ones(1,1),N-1) + diag(b*ones(1,1),-N+1);

[U,Delta] = eig(H);

%initial condition (vector -> pure state)
l_1 = 0;
l_2 = 3;

psi_Pos = exp(i*2*pi/N*l_1*[1:N])+exp(i*2*pi/N*l_2*[1:N]);
rho_Pos = psi_Pos'*psi_Pos;

% definition of the time vector for only one step rotation
% T_1 = 2*pi/N*(l_1-l_2)/(2*b*(cos(2*pi/N*l_1)-cos(2*pi/N*l_2)));
% t = linspace(0,T_1,Nt);

% initial condition (energy statistical mixture)
p = rand(1,N);
% I'm  taking advantage that U is real since it diagonalizes an Hermitian
% matrix -> ' = transpose
% rho_Pos = zeros(N);
% for n = 1:N
% rho_Pos = rho_Pos + p(n)*(U(:,n)*U(:,n)');
% end

%initial condition (density in position space)
% rho_Pos = diag([1,zeros(1,N-1)]); %|1><1| initialization
% rho_Pos = rand(N) + i*rand(N); % 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));

%important variables for the analytical solution
rho_En0 = U'*rho_Pos*U; %initial condition 
Bohr = repmat(diag(Delta),1,N) - repmat(diag(Delta)',N,1); % Bohr frequencies 


% generation of the sphere coordinates necessary for the movie 
[X,Y,Z] = sphere(100);

% coordinates of the position of the atoms
R = [cos([1:N]'*2*pi/N),sin([1:N]'*2*pi/N),zeros(N,1)];

% time evolution both in the position and energy basis
for n = 1:Nt 
    rhoEn = rho_En0.*exp(-i*Bohr*t(n));
    rhoPos = U*rhoEn*U';
    % vector of probabilities
    P = real(diag(rhoPos));
    
    % construction of the movie
    figure(1)
    for n = 1:N
    surf(P(n)*X+R(n,1),P(n)*Y+R(n,2),P(n)*Z+R(n,3),ones(size(X)))
    hold on
    plot3(R(n,1),R(n,2),R(n,3)+2,'k.')
    hold on
    view([30 90])
    axis equal
    axis off
    xlim([-3 3])
    ylim([-4 4])
    shading interp
    light
    lighting phong
    
    hold on
    
    end
    hold off
    pause(0.01);
end






