% @(#) 6.728 Homework Problem 3.3 -- Matlab Code
% ps3p3_movie.m file
% @(#) This program computes decomposes an the initial state of a particle
% in a box state into the M lowest energy eigenstates (i.e., cosines and
% sines). It then makes a movie of the time evolution.
% @(#) See Also: Eigenfunction Expansion Tutorial

% preamble

clear all; % Clear workspace
clc; % Clear the command window
close all; % Get rid of old figures
format compact; % suppresses extra lines

% Defines some standard constants used in quantum mechanics and
% statistical mechanics.
h = 6.62607E-34; % Planck constant [J*s]
hbar = h/(2*pi); % Rationalized Planck constant [J*s]
Me = 9.1094E-31; % Electron rest mass [kg]
Mp = 1.672623E-27; % Proton rest mass [kg]
Mn = 1.674929E-27; % Neutron rest mass [kg]
c = 2.99792458E8; % Speed of light in vacuum [m/s]

L = 1; % 1D Box [-1/2, 1/2] is of unit length
n_max = 20;

% Create our x vector
dx= 0.5E-3;
x = [-0.5:dx:0.5]';

% Create the column vectors psi.
psi = zeros(2001,1);
psi(876:1126) = sqrt(8);
% Initial state: psi(x,0) = \sqrt{\frac{8}{L}} if |x| \leq L/16
% else 0. x(876) = -0.5+875*0.0005 = -0.0625 = -L/16 as L = 1
% x(1126) = -0.5+1125*0.0005 = 0.0625 = L/16 as L = 1
psi_pdf = conj(psi).*psi;

% Create matrix of basis functions phi where energy level n is incremented
% column by column and location x is incremented row by row.
%
% phi_1(x(1)) phi_2(x(1)) ... phi_{n_max}(x(1))
% phi = phi_1(x(2)) phi_2(x(2)) ... phi_{n_max}(x(2))
% | | \ |
% phi_1(x(2001)) phi_2(x(2001)) ... phi_{n_max}(x(2001))
%
% in our case
% phi_{n odd} = sqrt{2/L) * cos(n*pi*x/L)
% phi_{n even} = sqrt(2/L) * sin(n*pi*x/L)
%
phi = zeros(2001,n_max); % allocate space for phi matrix
N = ones(2001,1)*[1:n_max];
X = x*ones(1,n_max);

if rem(n_max,2) == 0 % that is, if M is even
phi(:,1:2:n_max-1) = sqrt(2/L)*cos(N(:,1:2:n_max-1).*pi.*X(:,1:2:n_max-1)/L); % odd columns are cosines
phi(:,2:2:n_max) = sqrt(2/L)*sin(N(:,2:2:n_max).*pi.*X(:,2:2:n_max)/L); % even columns are sines
else % that is, if M is odd
phi(:,1:2:n_max) = sqrt(2/L)*cos(N(:,1:2:n_max).*pi.*X(:,1:2:n_max)/L); % Again, odd columns are cos
phi(:,2:2:n_max-1) = sqrt(2/L)*sin(N(:,2:2:n_max-1).*pi.*X(:,2:2:n_max-1).L); % even columns are sines
end

a = dx * phi' * psi; % Decomposition of initial state psi into the
% M lowest energy levels (cosines and sines)
a_pdf = abs(a).^2;

% Create our energy and time matrices

m = Me; % we'll say the particle is an electron
Eg = (hbar*pi)^2 / (2*m*L^2); % Ground state energy (n = 1)
varrho = m*L^2/h; % We'll measure time as fractions of the
t = varrho*[0:0.01:1]'; % reconstruction time mL^2 / h

Ecol = [(hbar*pi)^2/(2*m*L^2)]*[1:n_max].^2' ;

E = Ecol*ones(1,length(t));
T = ones(n_max,1)*t';

eta = exp(-i*E.*T/hbar);

a_t = (a*ones(1,length(t))).*eta;

psi_t = phi*a_t;
psi_t_pdf = conj(psi_t).*psi_t;

% Make movie
for j = 1:length(t);
plot(x,psi_t_pdf(:,j))
text(0.2,8,['Time t = ',num2str(t(j)/varrho),'mL^2/h'])
xlabel('Position x (in units of box length L)')
ylabel('|\Psi(x,t)|^2 (probability per unit length)')
axis([-0.5 0.5 0 max(max(psi_t_pdf))]);
M(j) = getframe;
end