% Animation of time evolution of a Gaussian wavepacket in an infinite
% potential well bounded between 0 and L.  

% Well parameters
L = 100;

% Particle parameters
m = 1;
hbar = 1;
x0 = 20; % initial position
sigmax = 5; % initial spread
p0 = 4*hbar/sigmax; % initial momentum
nmax = 100; % largest n component of the wavefunction

x = linspace(0,100,1000);
% initial wavefunction
psi0 = sqrt(sqrt(2*pi)*sigmax).*exp(-(x-x0).^2/(4*sigmax^2)).*exp(i*p0*x/hbar);
% basis wavefunctions, expansion coefficients, energies
psi = @(n) sqrt(2/L)*sin(n*pi*x/L);
c = @(n) sum(sqrt(2/L)*sin(n*pi*x/L).*psi0);
E = @(n) n^2*pi^2*hbar^2/(2*m*L^2); 

% probability density animation
clear M;
tlim = 500;
clf;
for t = 1:tlim
    psi_t = 0;
    for n = 1:nmax
        psi_t = psi_t + exp(-i*E(n)*t/hbar)*c(n)*psi(n);
    end
    plot(x,abs(psi_t).^2);
    axis([0,100,0,2000])
    ylabel('Probability density')
    xlabel('Position')
    title('Particle in a box')
    M(t) = getframe(gcf);
end

%v = VideoWriter('ParticleInBox.avi');
%open(v);
%writeVideo(v,M);
%close(v);

%figure
%movie(M)

% E. Boyda / SMC Phys125 / Fall '17. 