% @(#) 6.728 Homework Problem 3.1  -- Matlab Code
% ps3p1.m file
% @(#) This program computes a DFT by hand. It then performs a time evolution.
% @(#) See Also:  Eigenfunction Expansion Tutorial
% 6.728 -- Problem 3.1
% .m file written by Mark Schweizer
% modified by V. Anant

% preamble
diary ps3p1.out				% keep a log of our output
diary on;

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 = .5 *1E-9;				% Gaussian Width Paramters

% Create our x vector
dx=.2E-9;
x = [-20E-9:dx:20E-9]';

% Create the column vectors psi.
psi = (L^2*2*pi)^(-1/4) * exp(-(x / (2*L)).^2);
psi_pdf = conj(psi).*psi;

% Here we set up the q vector.  "dq" is determined by the periodicity
% of the DFT, as is the Max/Min q available.
N=size(x,1); % 				N = total number of points in x
dq=2*pi/(N*dx);
q=dq*[-(N-1)/2:(N-1)/2]';

% Compute our basis matrix and determine the expansion coefficients
phi = exp(i*(x*q'));
a   =  dx  * phi' * psi;
a_pdf = abs(a).^2;

% Create our energy and time matrices
% Note, that in all that follows, I've again combined a series of
% vectors into one matrix.  For example, the energy matrix defined
% below consists of columns of energy vectors as defined in the 
% tutorial.  Each column corresponds to a different time instance.
tau = Me*L^2/hbar
t = tau*[0 2 5 10 15]';
E = hbar^2/(2*Me) * q.^2;
eta = exp(-i* E/hbar*t.');

a_t = (a*ones(size(t'))).*eta;

psi_t = phi *a_t * dq /(2*pi);
psi_t_pdf = conj(psi_t).*psi_t;

% Plot |psi|^2 for each t.  
figure;
colormap('white');
plot(x/1E-9,psi_t_pdf);
% make sure cyclelines.m is in the same directory as this file
cyclelines;
axis([-10 10 0 max(max(psi_t_pdf))]);

% Create a vector with the analytical solution
m = Me;
psia = ones(size(x))*(((2*pi*L^2*(1+i*t*hbar/(L^2*2*m)).^2).^(-1/4)).') ...
  .*exp(-(1/(4*L^2))* (x.^2)*(1./(1+i*t*hbar/(2*m*L^2))).');
psia_pdf = conj(psia).*psia;

% Plot the analytical solution on the same graph
hold on;
plot(x/1E-9,psia_pdf,'x');

legend('t=0','t=2 tau','t=5 tau','t=10 tau','t=15 tau','Analytic solutions');
title('P3.1: Gaussian Wave Functions in Time');
xlabel('X [nm]');
ylabel('\Psi* \Psi  [m^{-1}]');

prob_sums_x= diag(psi_t'*psi_t)*dx
mean_x = (diag(psi_t'*diag(x,0)*psi_t)*dx)
mean_sq_x = (diag(psi_t'*diag(x.^2,0)*psi_t)*dx)
std_dev_x=sqrt(mean_sq_x - mean_x.^2)

prob_sums_q= diag(a_t'*a_t)*dq/(2*pi)
mean_q = (diag(a_t' * diag(q,0) * a_t ) * dq/(2*pi))
mean_sq_q = (diag(a_t' * diag(q.^2,0) * a_t )* dq/(2*pi))
std_dev_q=sqrt(mean_sq_q - mean_q.^2)

hup = std_dev_x.*std_dev_q

figure;
% we take the real part because the calculations above
% generates a small imaginary part due to machine precision
% errors.
plot(t/1e-15,real(std_dev_x)/1E-9);
title('P3.1: Standard Deviation in x-space');
xlabel('Time [fs]');
ylabel('\Delta x [nm]');

figure
plot(t/1e-15,real(std_dev_q)*1E-9);
title('P3.1: Standard Deviation in q-space');
xlabel('Time [fs]');
ylabel('\Delta q [nm^{-1}]');

diary off;