# Solving Schrodinger's equation: Free particles and uncertainty

# We'll be exploring the evolution of a Gaussian packet, corresponding 
# to the ground state of Hydrogen atom in a harmonic oscillator with 
# frequency one Terahertz. We start by defining some constants.


hbar = 1.054571628251774e-27;
omega = 1.e12;
protonMass = 1.672621637e-24;
m = protonMass;
a0 = sqrt(...); # RMS width of Gaussian

# We define psi0 on a lattice of Np points x, from -L/2 to L/2

L = ...;
Np = ...;
dx = ...;
x = linspace(-L/2,L/2-dx,Np);

# Now we define the initial wavefunction psi0 on our grid.
# Remember to use x.^2 to do the square pointwise (rather than as
# a matrix multiplication)

psi0 = (...)^(1./4.) * exp(-...);

plot(x, abs(psi0).^2)
xlim([-2*a0,2*a0])

# If psi(k,t) is expressed in Fourier space, applying the kinetic energy is 
# pointwise multiplication. p=−i hbar k, so p^2/2m=−ℏ^2 k^2/2m. We convert 
# from real space into Fourier space using a FFT (Fast Fourier Transform), 
# which returns psi(k) at Np points separated by dk=2 pi/L:
#    k=[0,dk,2dk,…,(Np/2)dk,−(Np/2−1)dk,…,−2dk,−dk]
# Notice that k grows until halfway down the list, and then jumps to 
# negative values and shrinks. We define k and k2=k2 as arrays
# Note that 'zeros' creates a matrix, so for a vector we initialize 
#    k=zeros(Np,1)
# Note that Matlab/Octave arrays start at one, so k(1)=0, k(2)=dk, ...

dk = ...;
k = zeros(1,Np);
for j  = 0:(...)
    k(j+1) = ... * dk;
end
for j = (...):(Np-1)
    k(j+1) = (j-Np) * ...;
end
k2 = k.^2;

# and we define the psi(t) using Fourier transforms.
# I still haven't figured out how to define functions
# like psi(t) without passing lots of variables like psi0, hbar, k2, ...

P = ...;

psiPby4 = ifft(exp(-...*(P/4.)/(...)).*fft(psi0));
plot(x,real(psiPby2), x, imag(psiPby2))
figure()
psiP = ...
plot(...)
figure()
...

# Calculate the width of our wavepackets at various times.
# The width is given by the sum of x.^2 .* abs(psi).^2 * dx

# Our initial wavepacket has RMS with DeltaX = a0. 
sqrt(sum(... .^2 .* abs(...).^2 * ...))
a0

# We calculate the widths for a series of times

times = linspace(0,2*P,101);
widths = 0.*times;
for j = 1:length(times)
    psi = ...
    widths(j) = ... abs(psi).^2 ...
end

# We know that a minimum uncertainty wavepacket like ours has
# Delta x Delta p=hbar/2, so we expect the packet width to grow like 
# vt with v given by the momentum uncertainty. We plot the comparison...

deltaP = ...
v = ...

plot(times, widths, times, ...)