%Mason Averill
%PHYS-442: Quantum Mechanics
%Fall 2021
%Computational Project

clear;
clc;

%Input Desired Constants
alpha=100;%changes initial wave function
k_0=15; %changes initial wave function
a=1; %dictates the size of the square well
m=10^-30;%mass of particle
n_max=30; %Number of terms to consider for solution
x_step_size=0.0001; %Specify resolution in x steps
t_step_size=2; %Specify resolution in time steps (seconds)
max_time=1000; %Maximum time for animation (seconds)

%Other Physical Constants
h_bar=6.62607015*10^-34/(2*pi);

%Initialize x and t vector 
x=0:x_step_size:a;
t=0:t_step_size:max_time;

%Solve for c_n up to specified number of terms
c_n=[];
for n=1:n_max
%Compute c_n

to_int=sin((n*pi/a).*x).*exp(-alpha.*(x-a/2).^2+1i*k_0.*(x-a/2));
int_value=trapz(x,to_int);

c_n(1,n)=sqrt(2/a)*(2*alpha/pi)^(1/4)*int_value;

end

%Compute the linear combination of the wave function

Si_t=[];
counter=1;
[row,column]=size(t);
while(counter<=column)
n=1;
Si_n=0;
for n=1:n_max    
Si_n=Si_n+c_n(1,n)*sqrt(2/a)*sin(((n*pi)/a).*x).*exp(-1i*((pi^2*h_bar*n^2)/(2*m*a^2))*t(1,counter));
end
Si_t=[Si_t;Si_n];
counter=counter+1;
end

%Si_n is si(x) for a specific value in time
%Si_t stores all Si_n, so each row in Si_t is for a new value in time 
Si_t_star=conj(Si_t);
Si_star_Si=Si_t_star.*Si_t;
shift_x=x-a/2;
check_normalization=trapz(shift_x,Si_star_Si(1,:));

%Begin Animation
p=plot(shift_x,Si_star_Si(1,:));
xlabel('Position in the Square Well')
ylabel('\psi^\ast\psi','FontSize',24)
xlim([-a/2 a/2])
ylim([0 round(max(Si_star_Si,[],'all'),12)])
text(a/4,.75*round(max(Si_star_Si,[],'all'),12), sprintf('alpha=%.1f \nk_0=%.1f \nn=%i \nm=%.2e \na=%.1f \nmax time=%.1f',[alpha, k_0,n,m,a,max_time]))
counter=1;

%normalization_store=[];
while(counter<=column)
title(sprintf('Probability Density Function of Particle vs Position in Well at Time=%.1f seconds',t(1,counter)))
%declare variable data source
p.YDataSource='Si_star_Si(counter,:)';
%check_normalization=trapz(shift_x,Si_star_Si(counter,:));
%normalization_store(1,counter)=check_normalization;
refreshdata
drawnow

counter=counter+1;
end

%max_area=max(normalization_store)
%min_area=min(normalization_store)
%average_area=mean(normalization_store)