%Infinite potential well probability density plotter - this program allows the
%user to see the probability density in position and momentum space for 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 position probability density as well as the momentum probability
%density.

%The final portion of the program allows the user to find the integrated
%probability between user-selected limits of position or momentum, but
%requires the use of a numerical-integration function such as MATLAB’s
%“trapz” function. 

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
xvec=0:xinc:a-xinc;%Position vector (0 is left side of well, "a" is right side)
xvec_norm=xvec/a;%Normalizes position vector by dividing by width of well
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*xvec)*exp(-i*En*t/hbar); %Makes position wavefunction
x_probden=conj(psi_x).*psi_x;%Finds position probability density
x_probden_norm=x_probden*a;%Normalizes position probability density

lambda_max=2*a;%Longest wavelength is twice the width of the well
kmin=2*pi/lambda_max;%Minimum wavenumber is 2pi over maximum wavelength
kmax=4*n*kmin;%Sets maximum wavenumber to 4n times the minimum wavenumber
kinc=0.001*kmax;%Sets wavenumber increment
kvec=-kmax+1000:kinc:kmax;%Makes wavenumber vector  
pvec=hbar*kvec;%Relates momentum vector(pvec) to wavenumber vector (kvec)
pvec_norm=pvec/(hbar*pi/a);%Normalizes momentum vector

pn=hbar*kn;%Relates quantized momentum to quantized wavenumber
p_probden=(hbar/(pi*a))*((pn^2-pvec.^2).^-2).*(2*(pn^2)...
    -2*pn^2*cos(pvec*a/hbar)*(-1)^n);%Finds momentum probability density
p_probden_norm=p_probden/(a/(pi*hbar));%Normalizes momentum prob density

figure('units','normalized','outerposition',[0.1 0.1 0.8 0.8]); %Opens and
                                            %sizes figure window
subplot(1,2,1);%Refers to left plot in horizontal group of 2 plots
plot(xvec_norm,x_probden_norm,'k');%Plots normalized position probability density
axis([0 1 0 1.2*max(x_probden*a)],'square');%Sets axis limits
grid on;%Draws vertical and horizontal grid lines
title(['Position Probability Density']);%Writes title on plot
xlabel('Normalized Position $(x/a)$','Interpreter','latex');%Labels x-axis
ylabel('Normalized Probability Density');%Labels y-axis

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

text(0.02,1.15*max(x_probden_norm),['m= ',num2str(m,3),' kg']);
                %Writes mass on plot
text(0.6,1.15*max(x_probden_norm),['a= ',num2str(a,3),' m']);
                %Writes well width on plot
text(0.02,1.1*max(x_probden_norm),['n= ',num2str(n)]);
                %Writes quantum number on plot
text(0.6,1.1*max(x_probden*a),['t= ',num2str(t,3),' s']);
                %Writes time on plot
h1=text(0.02,1.05*max(x_probden_norm),'Area under curve = 1.0');
                %Writes area under curve
hold on;%Allows later additions to this plot

subplot(1,2,2);%Refers to right plot of 2 horizontal plots
plot(pvec_norm,p_probden_norm,'k');%Plots normalized momentum prob den
axis([min(pvec_norm) max(pvec_norm) 0 1.2*max(p_probden_norm)],'square');
                                    %Sets axis limits
grid on;%Draws vertical and horizontal grid lines
title(['Momentum Probability Density']);%Adds title to plot
xlabel('Normalized Momentum $(p/(\hbar\pi/a))$','Interpreter','latex');
                        %Labels x-axis
ylabel('Normalized Probability Density');%Labels y-axis

text(0.95*min(pvec_norm),1.15*max(p_probden_norm),['m= ',num2str(m,3),' kg']);
                        %Writes mass on plot
text(0.4*max(pvec_norm),1.15*max(p_probden_norm),['a= ',num2str(a,3),' m']);
                        %Writes well width on plot
text(0.95*min(pvec_norm),1.1*max(p_probden_norm),['n= ',num2str(n)]);
                        %Writes quantum number on plot
text(0.4*max(pvec_norm),1.1*max(p_probden_norm),['t= ',num2str(t,3),' s']);
                        %Writes time on plot
h2=text(0.95*min(pvec_norm),1.05*max(p_probden_norm),'Area under curve = 1.0');
                        %Writes area under curve
hold on;%Allows later additions to this plot

%Note: No thick lines showing the outline of the well on the
%momentum-density plot

%The remainder of this program allows the user to compute the integrated
%probability between user-specified limits in position or momentum.  This is
%very useful, but it requires the implementation of a numerical integration
%function such as MATLAB's "trapz" function.

pause(5);%Allows a few seconds for user to view plots

list = {'Position','Momentum'};%List of two options for user selection
[indx,tf] = listdlg('PromptString',...%Produces dialog for user selection
    'To find integrated probability, select Position or Momentum',...
    'ListString',list,'SelectionMode','single','ListSize',[300,80]);

if indx==1; %If the user selects Position
    Probfunc=x_probden_norm;%Function to be integrated (y-values) set to
                            %normalized position probability density
    rangevec=xvec_norm;%x-values set to position vector
    prompt={'Starting Position (0 to 1)','Ending Position (0 to 1)'};
            %Prompts user for input
    title1='Position Range';%Title of dialog box
    dims=[1,60];%Dimensions of dialog box
    definput={'0','1'};%Default inputs
    params=inputdlg(prompt,title1,dims,definput,options);%Gets user inputs
                                %for position integration range
    integ_start=str2num(char(params(1))); %Lowest x-value of integration range
    integ_end=str2num(char(params(2))); %Highest x-value of integration range
    startindex=1+round(1000*integ_start);%Converts lowest value of position
                                          %to index
    endindex=round(1000*integ_end);%Converts highest value of position to index
    indices=startindex:endindex;%Range of indices to be integrated
    
    subplot(1,2,1);%Refers to left plot in horizontal group of 2 plots
    area(rangevec(indices),Probfunc(indices),'FaceColor','k');%Makes plot
            %of Probfunc with filled-in area under curve
    Totalprob=trapz(rangevec(indices),Probfunc(indices));%Performs numerical
             %trapezoidal-rule integration to find total probability
    delete(h1);%Deletes original text showing area under entire curve
    text(0.02,1.05*max(x_probden_norm),...
        ['Integrated Probability = ',num2str(Totalprob,3)]);
                %Writes area under user-selected portion of curve
end;
    
if indx==2; %If the user selects Momentum
    Probfunc=p_probden_norm;%Function to be integrated (y-values) set to
                            %normalized momentum probability density
    rangevec=pvec_norm;%x-values set to momentum vector
    prompt={['Starting Momentum (',num2str(min(rangevec),2),' to ',...
        num2str(max(rangevec),2),')'],...
        ['Ending Momentum (',num2str(min(rangevec),2),' to ',...
        num2str(max(rangevec),2),')']};%Prompts user for input
    title1='Momentum Range';%Title of dialog box
    dims=[1,60];%Dimensions of dialog box
    definput={num2str(min(rangevec),2),num2str(max(rangevec),2)};%Default inputs
    params=inputdlg(prompt,title1,dims,definput,options);%Gets user inputs
                                %for momentum integration range
    integ_start=str2num(char(params(1))); %Lowest x-value of integration range
    integ_end=str2num(char(params(2))); %Highest x-value of integration range
    startindex=1001+round((1000/(4*n))*integ_start);%Converts lowest value
                                    %of momentum to index
    endindex=1000+round((1000/(4*n))*integ_end);%Converts highest value
                                    %of momentum to index
    indices=startindex:endindex;%Range of indices to be integrated
    
    subplot(1,2,2);%Refers to right plot in horizontal group of 2 plots
    area(rangevec(indices),Probfunc(indices),'FaceColor','k');%Makes plot
            %of Probfunc with filled-in area under curve
    Totalprob=trapz(rangevec(indices),Probfunc(indices));%Performs numerical
             %trapezoidal-rule integration to find total probability
    delete(h2);%Deletes original text showing area under entire curve
    text(0.95*min(pvec_norm),1.05*max(p_probden_norm),...
        ['Integrated Probability = ',num2str(Totalprob,3)]);
              %Writes area under user-selected portion of curve
end;