%Finite potential well transcendental equation graph - This program allows
%the user to find the allowed values of energy for a quantum particle
%trapped in a finite potential well.  The user specifies the depth (V0) and
%width (a) of the well as well as the mass of the particle.  To find the
%allowed energies, this program solves a transcendental equation using a
%graphical technique.  Two functions, which in this program are called the
%"kappa/k"function and the "tangent" (for even solutions) or "negative
%cotangent" functions (for odd solutions) are both graphed on the same plot,
%and the intersections of the two curves can be used to find the allowed
%energies.

clear all;%Clears any existing data

options.WindowStyle='normal'; %Allows user to manipulate windows while
                                %waiting for user input
options.Interpreter='tex'; %Allows use of LaTeX notation in titles

prompt={'Depth of Well (V_0) in eV','Width of Well (a) in meters'};
                            %Prompts user for input of well parameters
title1='Finite Well Parameters';%Title of dialog box
dims=[1,60];%Dimensions of dialog box
definput={'60','2.5e-10'};%Default inputs for well depth and width
params=inputdlg(prompt,title1,dims,definput,options);%Gets user inputs for
    V0_eV=str2num(char(params(1))); %depth of potential well in electron-volts
    a=str2num(char(params(2))); %width of potential well in meters

V0norm=V0_eV/60;%Normalizes depth V0 to 60 eV
anorm=a/2.5e-10;%Normalizes width a to 2.5e-10 meters
m_suggest=(9.11e-31)/(V0norm*anorm^2);%Based on the user's choice of the
                %depth and width of the well, this mass will give the most
                %interesting solutions (but it's only suggested to the
                %user, who may enter any desired value).

prompt=['Mass (suggested range for this well is ',...
        num2str(0.1*m_suggest,1),' to ',num2str(10*m_suggest,1),' kg)'];
                        %Prompts user for mass and suggests range that will
                        %give interesting solutions
        title2='Mass of Particle';%Title of dialog box
        dims=[1,60];%Dimensions of dialog box
        definput={num2str(m_suggest)};%Default input for mass
        params=inputdlg(prompt,title2,dims,definput,options);%Gets user
                                            %input for mass of particle
        m=str2num(char(params(1))); %mass of particle in kg

%The following statements allow the user to select the type of solution to
%display - the choices are Even solutions, Odd solutions, or Both types.
list = {'Even','Odd','Both Even and Odd'};%Choices given to user
[indx,tf] = listdlg( 'PromptString','Type of Solution',...
    'ListString',list,'SelectionMode','single','ListSize',[140,80]);
V0=V0_eV*1.602e-19;%Converts V0 from eV to Joules
hbar=1.05457e-34;%Planck's modified constant
numsteps=2000;%Number of values to calculate
E=linspace(1e-6*V0,2*V0,numsteps);%Makes a vector of energy values starting
    %at V0/10^6 and ending at 2V0 with 2000 linearly spaced steps
k=sqrt((2*m/hbar^2)*E);%Makes a wavenumber vector proportional to the
                        %energy vector
kappa=sqrt((2*m/hbar^2)*(V0-E));%Makes a modified wavenumber vector needed
                                 %to find solutions
maxval_even=3*max(k*(a/2)/pi);%Sets maximum value at which the tangent
                                %function will be calculated for even
                                %solutions
ka_even(1,:)=linspace(0.001,pi/2-0.001,numsteps);%Makes the first portion
                    %of a linearly spaced vector from near zero to
                    %slightly less than pi/2.  This vector will be used as
                    %the argument in the tangent function to find solutions
                    %(the small offset of 0.001 prevents the tangent
                    %function from blowing up).
for count=2:round(maxval_even);%This loop makes the remaining portions of
                    %the ka_even vector, from pi/2 to pi, and from pi to
                    %3pi/2,etc.  The reason that the ka_even vector is
                    %calculated in chunks of pi/2 is that at values of
                    %pi/2, 3pi/2, 5pi/2, etc. the tangent function switches
                    %from large positive to large negative values, and we
                    %don't want that transition to be shown on the plot. So
                    %the first chunk stops just short of pi/2, the second
                    %chunk starts just after pi/2 and stops just short of
                    %pi, the third chunk picks up just after pi and stops
                    %just short of 3pi/2, and so forth.
    ka_even(count,:)=pi/2+ka_even(count-1,:);%Makes the later chunks
end;

%The following statements make the same kind of vector (ka_odd) for the odd
%solutions, which use the (negative) cotangent function instead of the
%tangent function.  The cotangent function transitions from large positive
%to large negative values at pi, 2pi, 3pi, etc., so these values are
%avoided.
maxval_odd=2*max(k*(a/2)/pi);%Sets the maximum value at which the cotangent
                            %function will be calculated for odd solutions
ka_odd(1,:)=linspace(0.001,pi-0.001,numsteps);%Makes the first portion
                    %of a linearly spaced vector from near zero to
                    %slightly less than pi.  This vector will be used as
                    %the argument in the cotangent function to find solutions
                    %(the small offset of 0.001 prevents the cotangent
                    %function from blowing up).
for count=2:round(maxval_odd);%This loop makes the remaining portions of the
                    %vector, from pi to 2pi, and from 2pi to 3pi, etc.
    ka_odd(count,:)=pi+ka_odd(count-1,:);%Makes the later chunks of ka_odd
end;

figure('units','normalized','outerposition',[0.1 0.1 0.8 0.8]); %Opens and
                                            %sizes figure window
subplot(1,2,1);%Refers to left subplot of 2 horizontal plots
plot(k*(a/2)/pi,real(kappa)./k,'k');%Plots the wavenumber function kappa/k
                                    %as a solid black line
grid on;%Draws vertical and horizontal grid lines
hold on;%Allows multiple curves on same graph
xplotmax=1.1*max(k*(a/2)/pi);%Sets maximum value of x-axis for plot
if (indx==1);%If the user wants to see even solutions
    plot(ka_even(1,:)/pi,tan(ka_even(1,:)),'--k');%Plots the tangent of the
                                    %first chunk of the ka_even vector
                                    %on the y-axis as a function of ka_even
    for count=2:round(maxval_even);%This loop counts through all chunks
                                    %of the ka_even vector
        plot(ka_even(count,:)/pi,tan(ka_even(count,:)),'--k');%Plots the
                                %later chunk of the tangent of ka_even
                                %vector on the y-axis as a function of
                                %ka_even
     
    end;

yplotmax1=10*real(kappa(round(numsteps/3)))./k(round(numsteps/3));%The
                    %y-axis plot limit will be set to 10 times the value
                    %of the kappa/k curve at a point that's between x=0
                    %and the point at which the kappa/k curve hits the x-axis.
axis([0 xplotmax 0 abs(yplotmax1)],'square');%Sets the axis limits
titletext='Finite Well Even Solutions';%Adds a title
end;

if (indx==2);%If the user selects odd solutions
    plot(ka_odd(1,:)/pi,-cot(ka_odd(1,:)),'-.k');%Plots the first chunk
                            %of the negative cotangent of the ka_odd vector
                            %on the y-axis as a function of ka_odd
    for count=2:round(maxval_odd);%Loop for plotting the later chunks
                                  %of the negative cotangent function
        plot(ka_odd(count,:)/pi,-cot(ka_odd(count,:)),'-.k');%Plots the
                                  %later chunks of the negative cotangent
    end;

    yplotmax1=10*real(kappa(round(numsteps/3)))./k(round(numsteps/3));
                %Sets the y-axis limit on the plot to 10 times the value of
                %the kappa/k function at a point between x=0 and the oint
                %at which the kappa/k curve hits the x-axis (that's why the
                %index is numsteps/3).
    axis([0 xplotmax 0 abs(yplotmax1)],'square');%Sets the axis limites
    titletext='Finite Well Odd Solutions';%Adds a title
end;
     
if (indx==3);%If the user wants to see both even and odd solutions
    plot(ka_even(1,:)/pi,tan(ka_even(1,:)),'--k');%Plots the first chunk of
                                    %the (even) tangent function                                            %
    for count=2:round(maxval_even);%Counts through the later chunks
        plot(ka_even(count,:)/pi,tan(ka_even(count,:)),'--k');%Plots the
                                    %later chunks of the tangent function
    end;
    
    plot(ka_odd(1,:)/pi,-cot(ka_odd(1,:)),'-.k');%Plots the first chunk of
                                    %the (odd) negative cotangent function
    for count=2:round(maxval_odd);%Counts through the later chunks
        plot(ka_odd(count,:)/pi,-cot(ka_odd(count,:)),'-.k');%Plots the 
                          %later chunks of the negative cotangent function
    end;

    yplotmax1=10*real(kappa(round(numsteps/3)))./k(round(numsteps/3));
                                %Finds the y-axis plot limit
    
    axis([0 xplotmax 0 abs(yplotmax1)],'square');%Sets the axis limits
    titletext='Finite Well Even and Odd Solutions';%Adds a title
end;

%The following statements add labels to the x- and y-axes and add a title
%to the left subplot
xlabel('$\frac{ka}{2\pi}$','Interpreter','latex','Fontsize',15);
ylabel('$\frac{\kappa}{k}$','Rotation',0,'Interpreter','latex','Fontsize',15);
title({titletext;['a= ',num2str(a),' m, V0= ',num2str(V0_eV),...
    ' eV, m= ',num2str(m),' kg']});

%The following statements provide an alternative graphical approach
%described in the book.  The x-axes are the same as in the approach used
%above, but the y-axes are (a/2)kappa (instead of kappa/k).  This approach
%will be shown in a second subplot to the right of the first.

subplot(1,2,2);%Refers to right subplot of 2 horizontal plots
plot(k*(a/2)/pi,(a/2)*real(kappa),'k');%Plots the kappa function
grid on;%Draws vertical and horizontal axes
hold on;%Allows multiple curves on same graph
yplotmax2=2*max((a/2)*real(kappa));%Finds maximum value of y-axis as twice
                                   %the largest value of the kappa function
axis([0 xplotmax 0 yplotmax2],'square');%Sets axis limits

if (indx==1);%If user wants to see even solutions
    plot(ka_even(1,:)/pi,ka_even(1,:).*tan(ka_even(1,:)),'--k');%Plots first
                                                %chunk of tangent function
    for count=2:maxval_even;%Counts through higher chunks
        plot(ka_even(count,:)/pi,ka_even(count,:).*tan(ka_even(count,:)),'--k');
                            %Plots higher chunks of tangent function
    end;
    titletext='Finite Well Even Solutions';%Adds title to plot
end;

if (indx==2);%user wants to see odd solutions
    plot(ka_odd(1,:)/pi,ka_odd(1,:).*-cot(ka_odd(1,:)),'-.k');%Plots first
                                    %chunk of negative cotangent function
    for count=2:maxval_odd;%Counts through higher chunks
        plot(ka_odd(count,:)/pi,ka_odd(count,:).*-cot(ka_odd(count,:)),'-.k');
                %Plots higher chunks of negative cotangent function
    end;
    titletext='Finite Well Odd Solutions';%Adds title to plot
end;

if (indx==3);%If user wants to see both even and odd solutions
    plot(ka_even(1,:)/pi,ka_even(1,:).*tan(ka_even(1,:)),'--k');%Plots
                           %the first chunk of the (even) tangent function
    for count=2:maxval_even;%Counts through the later (even) chunks
        plot(ka_even(count,:)/pi,ka_even(count,:).*tan(ka_even(count,:)),'--k');
                    %Plots the later (even) chunks of the tangent function
    end;
    plot(ka_odd(1,:)/pi,ka_odd(1,:).*-cot(ka_odd(1,:)),'-.k');%Plots the
                  %first chunk of the (odd) negative cotangent function
    for count=2:maxval_odd;%Counts through the later (odd) chunks
        plot(ka_odd(count,:)/pi,ka_odd(count,:).*-cot(ka_odd(count,:)),'-.k');
         %Plots the later (odd) chunks of the negative cotangent function
    end;
    titletext='Finite Well Even and Odd Solutions';%Add title to plot
end;

%The following statements add labels to the x- and y-axes and add a title
%to the right subplot
xlabel('$\frac{ka}{2\pi}$','Interpreter','latex','Fontsize',15);
ylabel('$\frac{a}{2}\kappa$','Rotation',0,'Interpreter',...
    'latex','Fontsize',15)
title({titletext;['a= ',num2str(a),' m, V0= ',num2str(V0_eV),...
    ' eV, m= ',num2str(m),' kg']});

