%Infinite potential well wavefunction plotter - this program allows the
%user to see the wavefunction solutions to the Schrödinger equation for the
%infinite potential well.  The user enters the value of the quantum number
%"n", the width of the potential well "a", the mass of the particle "m",
%and the time "t".  The program computes the quantum wavefunctions and
%displays the real and imaginary parts as well as the magnitude of that
%wavefunction.

clear all; %Clears any existing data

options.WindowStyle='normal'; %Allows user to manipulate windows while ...
                                %waiting for user input
                                
prompt={'Value of n (integer)','Width of Well (m)','Mass (kg)','Time (seconds)'};
            %Prompts user for input
title1='Infinite Potential Well';%Title of dialog box
dims=[1,60];%Dimensions of dialog box
definput={'1','1e-10','9.11e-31','0'};%Default inputs
params=inputdlg(prompt,title1,dims,definput,options);%Gets user inputs for
    n=str2num(char(params(1))); %quantum number
    a=str2num(char(params(2))); %width of potential well in meters
    m=str2num(char(params(3))); %mass of particle in kg
    t=str2num(char(params(4))); %time in seconds
xinc=0.001*a;%Position increment
x=0:xinc:a-xinc;%Position vector (0 is left side of well, "a" is right side
xnorm=x/a; %Normalized position (unitless)
hbar=1.0546e-34;%Planck's modified constant
kn=n*pi/a; %Quantized wave number
En=hbar^2*kn^2/(2*m); %Quantized energy
psi_x=sqrt(2/a)*sin(kn*x)*exp(-i*En*t/hbar); %Makes position wavefunction

figure('units','normalized','outerposition',[0.1 0.1 0.8 0.8]); %Opens and
                                            %sizes figure window

if max(real(psi_x))>max(imag(psi_x)); %Finds bigger wavefunction in order
    yaxis_lim=1.3*max(real(psi_x)); %to set y-limits of axes
else;
    yaxis_lim=1.3*max(imag(psi_x)); %Factor of 1.3 makes room for text
end;

subplot(1,3,1);%Refers to leftmost subplot of 3 horizontal plots
plot(xnorm,real(psi_x),'k'); %Plots real part of position wavefunction
axis([0 1 -yaxis_lim yaxis_lim],'square'); %Sets axis limits
grid on; %Draws vertical and horizontal grid lines
title('Real part of \Psi'); %Adds title to plot
xlabel('Position (x/a)'); %Labes x-axis
ylabel('Amplitude of \Psi (real)'); %Labels y-axis

%The following three statements make three thick lines to show outline of
%potential well on plot
line([0 0],[-yaxis_lim yaxis_lim],'Linewidth',3,'Color', 'k');
line([1 1],[-yaxis_lim yaxis_lim],'Linewidth',3,'Color', 'k');
line([0 1],[-yaxis_lim -yaxis_lim],'Linewidth',3,'Color', 'k');

text(0.02,0.9*yaxis_lim,['m= ',num2str(m,3),' kg']);%Writes mass on plot
text(0.6,0.9*yaxis_lim,['a= ',num2str(a,3),' m']);%Writes well width on plot
text(0.02,-0.9*yaxis_lim,['n= ',num2str(n)]);%Writes quantum number on plot
text(0.6,-0.9*yaxis_lim,['t= ',num2str(t,3),' s']);%Writes time on plot

subplot(1,3,2); %Refers to middle plot of 3 horizontal plots
plot(xnorm,imag(psi_x),'k');%Plots imaginary part of position wavefunction
axis([0 1 -yaxis_lim yaxis_lim],'square'); %Set axis limits
grid on;%Draws vertical and horizontal grid lines
title('Imaginary part of \Psi'); %Adds title to plot
xlabel('Position (x/a)'); %Labels x-axis
ylabel('Amplitude of \Psi (imaginary)'); %Labels y-axis

%The following three statements make three thick lines to show outline of
%potential well on plot
line([0 0],[-yaxis_lim yaxis_lim],'Linewidth',3,'Color', 'k');
line([1 1],[-yaxis_lim yaxis_lim],'Linewidth',3,'Color', 'k');
line([0 1],[-yaxis_lim -yaxis_lim],'Linewidth',3,'Color', 'k');

text(0.02,0.9*yaxis_lim,['m= ',num2str(m,3),' kg']);%Writes mass on plot
text(0.6,0.9*yaxis_lim,['a= ',num2str(a,3),' m']);%Writes well width on plot
text(0.02,-0.9*yaxis_lim,['n= ',num2str(n)]);%Writes quantum number on plot
text(0.6,-0.9*yaxis_lim,['t= ',num2str(t,3),' s']);%Writes time on plot

subplot(1,3,3);%Refers to rightmost plot of 3 horizontal plots
plot(xnorm,abs(psi_x),'k'); %Plots magnitude of position wavefunction
axis([0 1 0 1.2*max(abs(psi_x))],'square'); %Sets axis limits
grid on;%Draws vertical and horizontal grid lines
title('Magnitude of \Psi'); %Adds title to plot
xlabel('Position (x/a)'); %Labels x-axis
ylabel(['|\Psi|']); %Labels y-axis

%The following three statements make three thick lines to show outline of
%potential well on plot
line([0 0],[0 1.2*max(abs(psi_x))],'Linewidth',3,'Color', 'k');
line([1 1],[0 1.2*max(abs(psi_x))],'Linewidth',3,'Color', 'k');
line([0 1],[0 0],'Linewidth',3,'Color', 'k');

text(0.02,0.9*1.2*max(abs(psi_x)),['E_n= ',num2str(En,3),' J'])
                    %Writes energy on plot
text(0.02,0.05*1.2*max(abs(psi_x)),['n= ',num2str(n)]);%Writes quantum number on plot