%****************************************************************
% Program 6: Wavepacket propagation using exponential of H
% Original at https://arxiv.org/abs/0704.1622
% Modified: R. Herman 10/15/2016
%****************************************************************
clear
close all

% Parameters for solving the problem in the interval 0 < x < L
L = 100;                  % Interval Length
N = 400;                  % No of points
x = linspace(0,L,N)';     % Coordinate vector
dx = x(2) - x(1);         % Coordinate step

hbar = 1; m = 1;          % constants for Hamiltonian

% Parameters for making intial momentum space wavefunction phi(k)
ko = 2;                   % Peak momentum
a = 20;                   % Momentum width parameter
dk = 2*pi/L;              % Momentum step
km=N*dk;                  % Momentum limit
k=linspace(0,+km,N)';     % Momentum vector

% Obtain psi(x,0) from Gaussian wavefunction phi(k) using FFT
phi = exp(-a*(k-ko).^2).*exp(-i*6*k.^2); % unnormalized phi(k)
psi = ifft(phi);                   % integrate phi*expikx wrt x
psi = psi/sqrt(psi'*psi*dx);       % normalize the psi(x,0)

% Expectation value of energy; e.g. for the parameters
% chosen above <E> = 2.062.
avgE = phi'*0.5*diag(k.^2,0)*phi*dk/(phi'*phi*dk);

% CHOOSE POTENTIAL U(X): 
H = @(t) zeros(size(t)) + (t >= 0).*(1);
N0=N/2+1;
% U = step potential of height avgE that is located at x=L/2
% U = 0*H(x-(L/2));         % free particle wave packet evolution
% U = avgE*H(x-(L/2));      % scattering off step potential
                            % scattering through barrier width w
 x0=L/2; w=10; U = avgE*(H(x-x0)-H(x-(x0+w))); N0=ceil((x0+w)/dx);



% Finite-difference representation of Laplacian and Hamiltonian
e = ones(N,1); 
Lap = spdiags([e -2*e e],[-1 0 1],N,N)/dx^2;
H = -1/2*(hbar^2/m)*Lap + spdiags(U,0,N,N);

% Parameters for computing psi(x,T) at different times 0 < T < TF
NT = 200;                           % No. of time steps
TF = 29; T = linspace(0,TF,NT);     % Time vector
dT = T(2)-T(1);                     % Time step
P = max(abs(psi));                  
Usc = U/max(U)*P/2;                 % Rescale potential for plot

% Time displacement operator E=exp(-iHdT/hbar)
E = expm(-i*full(H)*dT/hbar);       % time evolution operator

%***************************************************************
% Simulate and plot probability density rho(x,T)
%***************************************************************
for t = 1:NT;                        
    psi = E*psi;                     % calculate new psi
    rho = conj(psi).*psi;            % rho(x,T) =psi*psi
    plot(x,rho,'b',x,Usc,'b--');     % plot rho(x,T) vs. x
    axis([0 L 0 0.15]);              % set x,y axis scale
    xlabel('x [m]', 'FontSize', 12);
    ylabel('probability density [1/m]','FontSize', 12);
    title('|\psi|^2 vs x')
    pause(0.05); 
end

% Calculate Reflection probability
R=0;
for a=1:N/2;
R=R+rho(a);
end

T=0;
for a=N0:N;
T=T+rho(a);
end
R=R*dx
T=T*dx