%Infinite potential well probability density plotter.  This program allows
%the user to add two or more infinite-well position wavefunctions together
%and to see how the combined wavefunction changes over time.

clear all;%Clears any existing data

options.WindowStyle='normal'; %Allows user to manipulate windows while ...
                                %waiting for user input

prompt={'Number of wavefunctions to add','Width of Well (meters)','Mass (kg)'};
            %Prompts user for input
title1='Infinite Potential Well';%Title of dialog box
dims=[1,60];%Dimensions of dialog box
definput={'2','1e-10','9.11e-31'};%Default inputs
params=inputdlg(prompt,title1,dims,definput,options);%Gets user inputs for
    numfuncs=str2num(char(params(1))); %number of wavefunctions to add
    a=str2num(char(params(2))); %width of potential well in meters
    m=str2num(char(params(3))); %mass of particle in kg

for count=1:numfuncs;%Loop to get quantum number and amount of each function
    prompt={['Quantum Number (n) of wavefunction #',num2str(count)],...
        ['Relative amount of wavefunction #',num2str(count),' (0 to 1)']};
            %Prompts user for input
title1='Wavefunction Numbers and Amounts';%Title of dialog box
dims=[1,60];%Dimensions of dialog box
definput={num2str(count),'0.5'};%Default inputs
params=inputdlg(prompt,title1,dims,definput,options);%Gets user inputs for
    n(count)=str2num(char(params(1))); %quantum number of wavefunction
    amt(count)=str2num(char(params(2))); %amount of wavefunction
end;  

total_amt=sum(amt);%Finds total of amounts for all wavefunctions
frac=sqrt(amt/total_amt);%Finds fraction of total for each wavefunction
hbar=1.05457e-34;%Planck's modified constant
xinc=0.01*a;%Position (x) increment
xvec=0:xinc:a+xinc;%Makes vector of position values
xvec_norm=xvec/a;%Normalizes position vector by dividing by width of well
kn=n*pi/a;%Finds wavenumber for each wavefunction
En=hbar^2*kn.^2/(2*m);%Finds energy for each wavefunction

%The following statements find the period of the composite wavefunction.
%If only two wavefunctions are added, the period is related to the
%energy difference between those two wavefunctions.  If three or more
%wavefunctions are added, the period depends on the least common multiplier
%(lcm) of the inverse square of the quantum numbers (n^-2).
if numfuncs==2;%Handles the case of only two wavefunctions being added
    omega=abs(En(2)-En(1))/hbar;%Finds angular frequency (omega) based on the
                                %energy difference between the wavefunctions 
    T=2*pi/omega;%Finds the period based on the angular frequency (omega)
    Sampletime=T/200;%Sets the time between samples as 1/200th of the period
else;%Handles the case of more than two wavefunctions being added
    factor_lcm=double(lcm(sym(n.^-2)));%This statement finds a factor needed
                    %for the period when 3 or more wavefunctions are added.
                    %To find the least common multiple of the inverse
                    %square of the quantum numbers, the MATLAB lcm function
                    %requires n^-2 to be converted into a symbolic expression.
                    %Once the lcm function operates has found the least
                    %common multiple of n^-2, that values is converted into
                    %a double-precision numerical value and stored as the
                    %variable factor_lcm.
    omega=min(En)/hbar;%Finds the smallest angular frequency based on the
                        %lowest energy
    T=2*pi*factor_lcm*min(n.^2)/omega;%Uses the least common multiple factor
                                        %and the smallest angular frequency
                                        %to find the period
    omega_max=max(En)/hbar;%Finds the biggest angular frequency of the
                            %wavefunctions
    Min_period=2*pi/omega_max;%Find the shortest period of the wavefunctions
    Sampletime=Min_period/20;%Sets the time between samples as 1/20th of the
                                %shortest period
end;

figure('units','normalized','outerposition',[0.25 0.42 0.5 0.5]);%Opens and
                                                %sizes figure window
          
for t=(0:Sampletime:T);%Makes time variable (t)
    for count=1:numfuncs;%Counts through the number of wavefunctions
                            %selected by the user
          psi_x(count,:)=frac(count)*sqrt(2/(a))*sin(kn(count)*xvec)...
              *exp(-i*En(count)*t/hbar);%Makes the position wavefunctions
                                %as 2-dimensional matrices.  The first
                                %index specifies the wavefunction number
                                %and the second index specifies the time,
                                %from t=0 to t=T with step size equal to
                                %Sampletime.
    end;
          psi_xtot=sum(psi_x);%Adds the wavefunctions together.  This sum
                                %operates on the first index, so it adds
                                %wavefunction 1 to wavefunction 2 (and
                                %wavefunctions 3, 4, etc. if the user has
                                %specified more than 2 wavefunctions).
                                
          x_probden_norm=a*(conj(psi_xtot).*psi_xtot);%The combined position
                                %wavefunction (psi_xtot) is complex, so the
                                %position probability density is formed by
                                %multiplying the wavefunction by its complex
                                %conjugate.  It is normalized by multiplying
                                %by the width of the potential well (a).
          plot(xvec_norm,x_probden_norm,'k');%Plots the normalized probabilty
                                %density as a function of the normalized
                                %position.
          axis([0 max(xvec_norm) 0 2*numfuncs]);%Sets the limits of the axes
          grid on;%Draws vertical and horizontal grid lines
          xlabel('Normalized Position (x/a)');%Labels the x-axis
          ylabel('Normalized Probability Density');%Labels the y-axis
          title('Infinite Potential Well Probability Density vs. Time');
          
          %The following statements add text information to the plots.  The
          %first "text" statement and the "for" loop write the combined
          %wavefunction as "Psi(x,t)=" followed by two or more wavefunctions
          %(fraction 1) psi_1 plus (fraction 2)psi_2 + ... (up to however
          %many wavefunctions the user specified).
          text(0.01,1.9*numfuncs,...
              ['\Psi(x,t) = ',num2str(frac(1),'%4.2f'),'\psi_{' int2str(n(1)) '}']);
          for count=2:numfuncs;
              text(0.1*count,1.9*numfuncs,...
                  ['+',num2str(frac(count),'%4.2f'),'\psi_{' int2str(n(count)) '}']);
          end;
          
         %The following two "text" statements write the time and the period
         %on the plot as the time increments upward from 0 to T.
          text(0.01,1.75*numfuncs,['t= ',num2str(t,2),' sec']);
          text(0.6,1.9*numfuncs,['Period= ',num2str(T,2),' sec']);
          
          %The following three statements make three thick lines to show outline of
          %potential well on the position-density plot
          line([0 0],[0 2*numfuncs],'Linewidth',3,'Color', 'k');
          line([1 1],[0 2*numfuncs],'Linewidth',3,'Color', 'k');
          line([0 1],[0 0],'Linewidth',3,'Color', 'k');
          pause(0.02);
end;

