% @(#) 6.728 Homework Problem 3.1 -- Matlab Code
% ps3p1mod.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, part a
% .m file originally written by Mark Schweizer
% modified by V. Anant
% Further modified 9/26/2006 by W.M. Kaminsky
% Corrected psi to agree with textbook/prob set convention (line 38 ff.)
% Made new code for analytical solution as old code seemed to disagree
% with numerical solution for all times t > 0. (lines 90-107)

% 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 = (pi*L^2)^(-1/4) * exp(-[x.^2]./[2*L.^2]);
%%%% WMK 9/26/2006 %%%%
% Code previously had psi = (L^2*2*pi)^(-1/4) * exp(-(x / (2*L)).^2);
% which implies abs(psi).^2 = (L^2*2*pi)^(-1/2) * exp(-[x.^2]./[2*L.^2]);
% This is a correctly normalized Gaussian with mean 0 and stand dev L
% However, PS3 lists psi as
% psi = (pi*L^2)^(-1/4) * exp(-[x.^2]./[2*L.^2]);
% --> abs(psi).^2 = (pi*L^2)^(-1/2) * exp(-[x.^2]./[L.^2]);
% which is a correctly normalized Gaussian with mean 0 and stand dev
% L/sqrt(2). Personally, I think the former convention makes more sense
% than the latter since abs(psi).^2 is the actual probability
% distribution. Nevertheless, I corrected m-file to agree with textbook
% and problem set convention.
%%%%%%%%%%%%%%%%%%%%%%%
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;

%%%%% WMK 9/26/2006 %%%%%%%%
%% Make my own code for plotting analytical solution since it seems to me
%% that the original code, namely
%
%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;
%
%% has some sort of error as it only agrees well with numerical solution at
%% time = 0.

X = x*ones(1,length(t)); % Make matrices (X,T) with x incremented row by
T = ones(length(x),1)*t'; % row and t incremented column by column

var_t = (0.5*L^2)*[1 + (hbar/(m*L^2))^2 * T.^2]; %% Textbook Equ 4.20
psia_pdf = [(2*pi*var_t).^(-1/2)].*exp(-[X.^2]./(2*var_t)); %% Equ. 4.21

%%%% END WMK 9/26/2006 %%%%%%%%


% 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,'-',t/1e-15,1E9*sqrt(var_t(1,:)),'x');
title('P3.1: Standard Deviation in x-space');
xlabel('Time [fs]');
ylabel('\Delta x [nm]');
legend('Numerical solution','Analytical solution')

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;