% A simple MATLAB code to plot the momentum
% wave function for a particle in an escape-proof "box"

clear all

hbar = 1.054e-34;
L = 101e-10;

% Set up k' and p' axes

k_prime = linspace(-2*pi/L,2*pi/L,100);
dk_prime = k_prime(2) - k_prime(1);
p_prime = hbar * k_prime;
dp_prime = p_prime(2) - p_prime(1);

% k1 and p1 

k1 = (pi/L); p1 = hbar * k1;

% Fourier transform A(k')

A1 = 0.5*sqrt(L/pi)*( sinc((k_prime+k1)*L/(2*pi)) ); 
A2 = 0.5*sqrt(L/pi)*( sinc((k_prime-k1)*L/(2*pi)) ); 
A = A1 + A2;
A_sum = sum(abs(A).^2)*dk_prime;

% Momentum wave function Phi(p')

Phi1 = 0.5*sqrt(L/(pi*hbar))*( sinc((p_prime+p1)*L/(2*pi*hbar)) );
Phi2 = 0.5*sqrt(L/(pi*hbar))*( sinc((p_prime-p1)*L/(2*pi*hbar)) );
Phi = Phi1 + Phi2;
Phi_sum = sum(abs(Phi).^2)*dp_prime;

% Normalized p' axis for plotting purposes; the points of classical
% momenta occur when the variable pp_N is plus or minus unity 

pp_N = p_prime / p1;

figure(1); clf; 
h = plot(pp_N,Phi1,'kx',pp_N,Phi2,'ko',pp_N,Phi,'k--');
grid on; set(h,'linewidth',[2.0]); set(gca,'Fontsize',[18]);
xlabel('NORMALIZED MOMENTUM  [p^{\prime}/p_1]'); 
ylabel('WAVE FUNCTION  [sqrt(s / kg m)]');
legend('Phi1','Phi2','Phi');

figure(2); clf; h = plot(pp_N,abs(Phi).^2,'k-');
grid on; set(h,'linewidth',[2.0]); set(gca,'Fontsize',[18]);
xlabel('NORMALIZED MOMENTUM  [p^{\prime}/p_1]'); 
ylabel('PROBABILITY DENSITY  [s / kg m]');