% free_particle.m   :   Wave packet for free particle
clear all; 
close all; 
format long;
% table with x-values
x0 = -6; XL = 5; x2 = 4*XL; x1=-x2;  NX = 10000;
x = x1 : (x2-x1)/NX : x2; 
% central wave number k0
k0 = 6;
% formula for f at t=0   %  (x0 is the initial position, minimal uncertainty at t=0)
commonfactor = sqrt(pi)*exp(-(x-x0).^2/4)/5;
rf0 = cos(k0*(x-x0)).*commonfactor;                     % real part at t = 0
if0 = sin(k0*(x-x0)).*commonfactor;                     % imaginary part at t = 0
fullscreen = get(0,'ScreenSize');
% fullscreen is a matrix. fullscreen(1,2) are the pixel coordinates of the lower left corner, 
% fullscreen(3) is the width, fullscreen(4) is the height, all in pixels 

figure('Position',[0 0 0.8*fullscreen(3) 0.8*fullscreen(4)])  % determines the position and size on the screen
  p = plot3(x,rf0,if0,'-','EraseMode','xor','linewidth',2);
grid;
axis([-10 20 -0.3 0.3 -0.3 0.3]);
view(-20,30);   % "angles of observation"
daspect([10 1 1])
hold on;
%figure
%plot(x,f0,x,fh,x,-fh);
xlabel('x');
ylabel('real part');
zlabel('imaginary part');
title('free particle wave packet with dispersion; initial state included for comparison (red curve)');
%
t0 =0; T =10/k0; NT = 400;	   
t = t0 : T/NT : T;
prefactor = sqrt(pi./(1+i*t)).*exp(-i*t*k0^2)/5;
for nt = 1 : length(t)                                 % loop in time
	  f = prefactor(nt)*exp(i*k0*(x-x0)).*exp(-(x-x0-2*k0*t(nt)).^2/(4*(1+i*t(nt))));	     
set(p,'XData',x,'YData',real(f),'ZData',imag(f))       % erases old curve, displays new one
drawnow
	   pause(0.05);  % increase if you want a "slower" animation
end % t-loop
line(x,rf0,if0,'color','r');                           % initial state 
hold off;
