% @(#) 6.728 Demo code -- Matlab Code
% @(#) This code shows an animated gaussian wave packet in a SHO
% @(#) potential.  The width of the initial gaussian wave packet
% @(#) can be changed to see "breathing" effects.  Also, <x>
% @(#) and the variance of x are plotted vs. time to see
% @(#) the classical result of sinusiodal motion.


% 6.728 -- Simple Hamronic Oscillator
% Mark Schweizer (mschweiz@mit.edu)

clear all;				% Clear workspace
clc;					% Clear the command window
close all;				% Get rid of old figures
wavename=['Gaussian Wave'];		% Used for titles in plots
MOVIE = 1;				% Make this zero to skip the 
					% movie generation


% 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]

% 1 nm is the the first eigenfunction for this SHO.
L = input('What is the gaussian width paramter, L, in nanometers? ')*1E-9;


% Create our x vector
dx=L/8;
x = [-12*1E-9:dx:12*1E-9]';

% Define our initial position
x0 = L*input('What is the initial displacement in units of L? ');

% Create the column vector psi.
psi = (L^2*2*pi)^(-1/4) * exp(-((x-x0) / (2*L)).^2);
 
% Here we set up the q vector.  "dq" is determined by the periodicity
% of the DFT, as is the Max/Min q available.


w0=.5*hbar/Me/(1e-9)^2 ;      % Picked so ground state gaussian has L=1E-9;
g = sqrt(Me*w0/hbar)*x;       % Intermediate value.  Morrison calls this Zeta.

fprintf('Computing hermite polynomials ...\n');
h = herm(25);                  % Computes 1st 25 hermite polynomials

% Compute our basis matrix
fprintf('Computing basis matrix ...\n');
phi = zeros(length(x),size(h,1));
for j =[1:size(h,1)]
  A =  (Me*w0/(hbar*pi))^(1/4)*1/sqrt( prod([1:(j-1)])*2^(j-1));
  phi(:,j) = A*polyval(h(j,:),g).*(exp(-(g.^2)/2));
end;

fprintf('Computing the time dependent wave function...\n');
a   =  dx  * phi' * psi;

% Create our energy and time matrices
tau = 2*pi/w0;      %  Chosen to show one full periods of oscillations
t = tau*[0:.02:2]'; %  10 steps per period
E = w0*hbar*([0:(size(phi,2)-1)]+.5)';
eta = exp(-i* E/hbar*t.');

% Create our time dependent expansion coefficients
a_t = (a*ones(size(t'))).*eta;

% Find our time dependent psi
psi_t = phi *a_t ;
psi_t_pdf = conj(psi_t).*psi_t;

fprintf('Computing Statistical Quantities...\n');
% Compute the statisitical quantitites of interest.
mx = real(diag(psi_t'*diag(x)*psi_t))*dx;
mx2=real(diag(psi_t'*diag(x.^2)*psi_t))*dx;
msdx=sqrt(mx2-mx.^2);


if MOVIE == 1
fprintf('Press any key to generate an animation.\n');
pause;
fprintf('Generating the movie sequence...\n');
% Now we'll generate an animated movie
ymax =1.1*max(max(psi_t_pdf));
pota =ymax / 12E-9^2;
pot  =pota*x.^2;
frames=(length(t)-1)/4;
figure;
M = moviein(frames);
for f=1:frames
  plot(x/1E-9,psi_t_pdf(:,f)); hold on;
  plot(x/1E-9,pot,'g');hold off;
  axis([-12 12 0 ymax]);
  M(:,f) = getframe;
end
delete(gcf);

fprintf('Press any key to view the animation.\n');
pause;
% Create a title axis to play the movie in
figure;
axis;
mh=gca;                                 % Save the axis handle for later
plot(x/1E-9,psi_t_pdf(:,1));		% Setup for proper axis labels
title([wavename,' Function in Time'],'FontSize',18);
axis([-12 12 0 ymax]);
xlabel('X [nm]','FontSize',18);
% 00.10.09 ekd
% sylabel('\18 {\symbol y^*y}  [m^{-1}]');
ylabel('\Psi^* \Psi  [m^{-1}]');
movie(mh,M,-4,6);			% Play the movie 3 times
end

fprintf('Press any key to generate a x-t color plot.\n');
pause;
% Now, we'll generate a two dimensional plot of x vs t, but the
% color will indicate psi^2. The brighter the color, the higher the
% probability!
figure
surf(x/1E-9,t/tau,psi_t_pdf');
view(2);
brighten(.7);
shading flat;
title([wavename,' Function in Time'],'FontSize',18);
xlabel('X [nm]','FontSize',18);
% 00.10.09 ekd
% sylabel('\18 t [{\symbol t}]');
ylabel('t [\tau]');

fprintf('Press any key to plot the statistical qunatitites.\n');
pause;
figure
plot(t/tau,mx);
title('<x> vs. t','FontSize',18);
ylabel('x [nm]','FontSize',18);
% 00.10.09 ekd
% sxlabel('\18 t [{\symbol t}]');
xlabel('t [\tau]');

%label changed by T. Orlando to read x^2, 29-Sept-00


figure
plot(t/tau,msdx.^2);
%stitle('\18 ({\symbol D} x)^2 vs. t');
title('(\Delta x)^2 vs. t');
ylabel('x^2 [nm^2]','FontSize',18);
% 00.10.09 ekd
% sxlabel('\18 t [{\symbol t}]');
xlabel('t [\tau]');