\documentclass{article}

\usepackage{graphicx}
\usepackage{setspace}
\usepackage{listings}
\usepackage{color}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{float}
\usepackage[justification=centering]{caption}
\usepackage{subcaption}
\usepackage[margin=0.75in]{geometry}
\usepackage[parfill]{parskip} 
\usepackage{fancyhdr}
\usepackage{titlesec}
\usepackage{lipsum}
\usepackage{enumitem}
\usepackage{cancel}

\titleformat{\subsection}[runin]{\normalfont \large \bfseries}
    {\thesubsection}{1em}{}

\renewcommand{\thesubsection}{\indent(\alph{subsection})}

    \definecolor{dkgreen}{rgb}{0,0.6,0}
    \definecolor{gray}{rgb}{0.5,0.5,0.5}
    \definecolor{mauve}{rgb}{0.58,0,0.82}


    \definecolor{mygreen}{RGB}{28,172,0} % color values Red, Green, Blue
    \definecolor{mylilas}{RGB}{170,55,241}


    \lstset{language=Matlab,%
        %basicstyle=\color{red},
        breaklines=true,%
        morekeywords={matlab2tikz},
        keywordstyle=\color{blue},%
        morekeywords=[2]{1}, keywordstyle=[2]{\color{black}},
        identifierstyle=\color{black},%
        stringstyle=\color{mylilas},
        commentstyle=\color{mygreen},%
        showstringspaces=false,%without this there will be a symbol in the places where there is a space
        numbers=left,%
        numberstyle={\tiny \color{black}},% size of th_sume numbers
        numbersep=9pt, % this defines how far the numbers are from the text
        emph=[1]{for,end,break},emphstyle=[1]\color{red}, %some words to emphasise
        %emph=[2]{word1,word2}, emphstyle=[2]{style},    
        xleftmargin=2em
    }
    %\lstset{basicstyle=\small,
    %        keywordstyle=\color{mauve},
    %        identifierstyle=\color{dkgreen},
    %        stringstyle=\color{gray},
    %        numbers=left,
    %        xleftmargin=5em
    %        }


\title{ECE 456 - Problem Set 2}
\date{2021-03-01}
\author{David Lenfesty \\ lenfesty@ualberta.ca
    \and Phillip Kirwin \\ pkirwin@ualberta.ca}
    
\pagestyle{fancy}

\fancyhead[L]{\textbf{ECE 456} - Problem Set 2}
\fancyhead[R]{David Lenfesty and Phillip Kirwin}
\fancyfoot[C]{Page \thepage}
\renewcommand{\headrulewidth}{1pt}
\renewcommand{\footrulewidth}{1pt}

\begin{document}
    \doublespacing
    \maketitle
    \thispagestyle{empty}
    \singlespacing
    \newpage
    hello
    \section*{Problem 1}

        \begin{enumerate}[align=left,leftmargin=*,labelsep=1em,label=\bfseries(\alph*)] 
            \item %1a
            Code: 
            \begin{lstlisting}[language=Matlab]
clear all;
%physical constants in MKS units

hbar = 1.054e-34;
q = 1.602e-19;
m = 9.110e-31;

%generate lattice

N  = 100;                   %number of lattice points
n  = [1:N];                 %lattice points
a  = 1e-10;                 %lattice constant
x  = a * n;                 %x-coordinates
t0 = (hbar^2)/(2*m*a^2)/q;  %encapsulating factor
L  = a * (N+1);             %total length of consideration

%set up Hamiltonian matrix 

U = 0*x; %0 potential at all x
main_diag  = diag(2*t0*ones(1,N)+U,0); %create main diagonal matrix
lower_diag = diag(-t0*ones(1,N-1),-1); %create lower diagonal matrix
upper_diag = diag(-t0*ones(1,N-1),+1); %create upper diagonal matrix

H = main_diag + lower_diag + upper_diag; %sum to get Hamiltonian matrix

[eigenvectors,E_diag] = eig(H); %"eigenvectors" is a matrix wherein each 
                                %column is an eigenvector
                                %"E_diag" is a diagonal matrix where the
                                %corresponding eigenvalues are on the
                                %diagonal.
                                
E_col = diag(E_diag); %folds E_diag into a column vector of eigenvalues

% return eigenvectors for the 1st and 50th eigenvalues

phi_1  = eigenvectors(:,1);
phi_50 = eigenvectors(:,50);

% find the probability densities of position for 1st and 50th eigenvectors

P_1 = phi_1 .* conj(phi_1);
P_50 = phi_50 .* conj(phi_50);

% Find first N analytic eigenvalues
E_col_analytic = (1/q) * (hbar^2 * pi^2 * n.*n) / (2*m*L^2);

% Plot the probability densities for 1st and 50th eigenvectors

figure(1); clf; h = plot(x,P_1,'kx',x,P_50,'k-');
grid on; set(h,'linewidth',[2.0]); set(gca,'Fontsize',[18]);
xlabel('POSITION [m]'); ylabel('PROBABILITY DENSITY [1/m]');
legend('n=1','n=50');

% Plot numerical eigenvalues
figure(2); clf; h = plot(n,E_col,'kx'); grid on;
set(h,'linewidth',[2.0]); set(gca,'Fontsize',[18]);
xlabel('EIGENVALUE NUMBER'); ylabel('ENERGY [eV]');
axis([0 100 0 40]);

% Add analytic eigenvalues to above plot

hold on;
plot(n,E_col_analytic,'k-');
legend({'Numerical','Analytical'},'Location','northwest');

            \end{lstlisting}
            %
            \begin{figure}[H]
                \centering
                \begin{subfigure}{0.4\linewidth}
                    \centering
                    \includegraphics[width=\textwidth]{q1a_1and50.png}
                    \caption{}
                    \label{fig:q3_ne}
                \end{subfigure}%
                \begin{subfigure}{0.4\linewidth}
                    \centering
                    \includegraphics[width=\textwidth]{q1a_eigenvals.png}
                    \caption{Comparison of first 101 numerical and analytic eigenvalues.}
                    \label{fig:q3_current}
                \end{subfigure}
                \caption{(a) Probability densities for \(n=1\) and \(n=50\). 
                            (b) Comparison of first 101 numerical and analytic eigenvalues.}
            \end{figure}
            \item %1b
            \begin{enumerate}[align=left,leftmargin=*,labelsep=0em,label=(\roman*)] 
                \item %1bi
                The analytical solution is

                \begin{equation}
                    \phi(x) = A sin \left( \frac{n \pi}{L} x \right)
                    \label{eq:analytical}
                \end{equation}

                In order to normalise this equation it must conform to the following:

                \begin{equation}
                    \int_{0}^{L} \left| \phi(x) \right| ^2 dx = 1.
                    \label{eq:normalise}
                \end{equation}

                We use the following identity:

                \begin{equation}
                    \int\sin^2{(a x)}\:dx = \frac{1}{2} x -\frac{1}{4a}\sin{(2a x)}.
                    \label{eq:integral_ident}
                \end{equation}

                Given that $sin$ of a real value is always real, we can disregard the norm operation, and directly relate (\ref{eq:analytical}) 
                to the above identity.
                Evaluating the integral gives us the following relationship

                \begin{equation*}
                    \frac{1}{A^2} = \frac{1}{2} L - \cancel{\frac{L}{4n \pi} sin \left( \frac{2n \pi}{L} L  \right)} - \cancel{\frac{1}{2} \cdot 0} + \cancel{\frac{L}{4n \pi} sin(0)}
                \end{equation*}

                From this, we find:

                \begin{equation*}
                    A = \sqrt{\frac{2}{L}}
                \end{equation*}

                \item %1bii

                Starting with the normalization condition for the numerical case:
                \begin{align}
                    a\sum_{l=1}^N\left|\phi_l\right|^2 = a \\
                    a\sum_{l=1}^N\left|B\sin{\left(\frac{n\pi}{L}x_l\right)}\right|^2 = a,
                    \label{eq:numerical_norm}
                \end{align}
                recalling that \(x = al\), and allowing \(a \to 0\) while holding \(L\) 
                constant implies that \(N \to \infty\), since \(a=\frac{L}{N}\). 
                An integral is defined as the limit of a Riemann sum as follows:
                \begin{equation}
                    \int_c^df(x)\:dx \equiv \lim_{n \to \infty}\sum_{i=1}^n\Delta x\cdot f(x_i).
                    \label{eq:Riemann_sum}
                \end{equation}
                where \(\Delta x = \frac{d-c}{n}\) and \(x_i = c + \Delta x \cdot i\). In our case,
                \(n = N\), \(i = l\), \(c = 0\), \(d = L\), and \(\Delta x = a\), \(x_i = x_l\), 
                \(f(x) = \left|B\sin{\left(\frac{n\pi}{L}x\right)}\right|^2\). Therefore we can write
                \begin{equation*}
                    \int_0^L \left|B\sin{\left(\frac{n\pi}{L}x\right)}\right|^2\:dx 
                    = \lim_{N \to \infty}\sum_{l=1}^N a \cdot \left|B\sin{\left(\frac{n\pi}{L}x_l\right)}\right|^2
                    = a.
                \end{equation*}
                
                Using (\ref{eq:integral_ident}), we have
                \begin{equation*}
                    \int_0^L \left|B\sin{\left(\frac{n\pi}{L}x\right)}\right|^2\:dx 
                    = \frac{1}{2}L - \cancel{\frac{L}{4n\pi}\sin{\left(\frac{2n\pi}{L}L\right)}}
                        - 0 + 0
                    = \frac{a}{B^2}.
                \end{equation*}

                This means that B must be
                \begin{equation*}
                    B = \sqrt{\frac{2 a}{L}} = \sqrt{a} \times A.
                \end{equation*}

            \end{enumerate}

            \item %1c
            \begin{enumerate}[align=left,leftmargin=*,labelsep=0em,label=(\roman*)]
                \item %1ci

                % TODO check that B is expressed correctly in all these equations (alpha, not ell, wait for pdf gen)
                From the base form of $\phi _\ell = Bsin \left( \frac {n \pi}{L} a \ell \right)$, we can see that
                $\phi _ {\ell + 1}$ and $\phi _ {\ell - 1}$ correspond to the trigonometric identities
                $sin(a + B) = sin(a)cos(B) + cos(a)sin(B)$ and
                $sin(a + B) = sin(a)cos(B) + cos(a)sin(B)$, respectively, where
                $a = \frac{n \pi a \ell}{L}$ and $B = \frac{n \pi a}{L}$.

                Plugging these identities into equation (7) from the assignment and simplifying, we get to this equation:

                \begin{equation*}
                    -t_0 B sin \left( \frac{n \pi a \ell}{L} \right) + 2 t_0 \phi _{\ell} - t_0 sin \left( \frac{n \pi a \ell}{L} \right)
                \end{equation*}

                At this point, we notice that $\phi _ \ell = B sin \left( \frac {n \pi}{L} a \ell \right)$,
                so we can factor it out.

                With some minor rearranging, this leaves us with the final expression for $E$:

                \begin{equation}
                    E = 2t_0 \left( 1 - cos \left( \frac{n \pi a}{L} \right) \right)
                    \label{eq:numerical_E}
                \end{equation}
                
                \item %1cii

                \begin{figure}[H]
                    \centering
                    \includegraphics[width=\textwidth]{q1cii.png}
                    \caption{Analytic solution to numeric system, plotted.}
                \end{figure}

                We can see here that the "predicted" numerical response matches nearly exactly the actual
                calculated numerical solution.

                \item %1ciii

                Applying the approximation $cos(\Theta) = 1 - \frac{\Theta^2}{2}$ for small $\Theta$,
                on equation (\ref{eq:numerical_E}), we get the following expression

                \begin{equation*}
                    E = 2 t_0 \left( \frac{n^2 \pi^2 a^2}{2 L^2} \right)
                \end{equation*}

                Fully substituting the known value of $t_0$, and we can get our final expression for $E$:

                \begin{equation*}
                    E = \frac{ \hbar^2 n^2 \pi^2 }{ 2 m L^2 }
                \end{equation*}

                \item %1civ

                \begin{figure}[H]
                    \centering
                    \begin{subfigure}{0.5\textwidth}
                        \centering
                        \includegraphics[width=\textwidth]{q1civ_fig1.png}
                        \caption{Probability densities for \(n=1\)\\and \(n=50\).}
                    \end{subfigure}%
                    \begin{subfigure}{0.5\textwidth}
                        \centering
                        \includegraphics[width=\textwidth]{q1civ_fig2.png}
                        \caption{Comparison of first 101 numerical and analytic eigenvalues.}
                    \end{subfigure}
                \end{figure}

                With the decreased lattice spacing, and increased number of points, we can see the numerical solution
                of eigenvalues matches the analytical solution much closer. As well, we can see that the probability
                density function appears to be "squeezed" in the centre, for the $n = 50$ case. The \(n=50\) case is 
                now a constant-amplitude wave, which corresponds to the expected analytic result - in contrast to the 
                plot in section (a), which has a low-frequency envelope around it.

            \end{enumerate}
            \item %1d
            \begin{enumerate}[align=left,leftmargin=*,labelsep=0em,label=(\roman*)]
                \item %1di

                In order to modify the computations to those for a particle in a "ring"
                we simply had to add $-t_0$ elements as the "corner" elements of the
                hamiltonian operator array:
                \begin{lstlisting}[language=Matlab]
% Modify hamiltonian for circular boundary conditions
H(1, N) = -t0;
H(N, 1) = -t0;
                \end{lstlisting}
                
                \begin{figure}[H]
                    \centering
                    \begin{subfigure}{0.5\textwidth}
                        \centering
                        \includegraphics[width=\textwidth]{q1di_fig1.png}
                        \caption{Probability densities for \(n=4\)\\and \(n=5\).}
                    \end{subfigure}%
                    \begin{subfigure}{0.5\textwidth}
                        \centering
                        \includegraphics[width=\textwidth]{q1di_fig2.png}
                        \caption{Comparison of first 101 numerical and analytic eigenvalues.}
                    \end{subfigure}
                \end{figure}

                \item %1dii

                The energy levels for eigenvalues number 4 and 5 are both $0.06eV$.
                These eigenstates are degenerate because they both have the same 
                eigenvalue/energy (0.06).

                \item %1diii

                \begin{figure}[H]
                    \centering
                    \includegraphics[width=\textwidth, angle=-90]{q1div.jpg}
                    \caption{Sketch of appropriate energy levels}
                \end{figure}

                \item %1div

                Plugging the valid levels for $n$ into equation (10) from the assignment (and dividing
                by the requisite $q$), we get the energy levels of $0eV, 0.0147eV, and 0.0589eV$, for the
                eigenvalue numbers $n = 0, 1, and 2$, respectively. These match very closely, to within
                acceptable margin of the numerical results from part (ii).


            \end{enumerate}

        \end{enumerate}


    \section*{Problem 2}

        \begin{enumerate}[align=left,leftmargin=*,labelsep=1em,label=\bfseries(\alph*)] 
           
            \item %2a
            Since \(t_0 = \frac{\hbar^2}{2ma^2}\), and \(U_l = -\frac{q^2}{4\pi\varepsilon_0r_l}
            + \frac{l_o(l_o + 1)\hbar^2}{2mr_l^2}\), the middle diagonal elements will have values
            \begin{equation*}
                \hat{H}_{ll} = \frac{\hbar^2}{ma^2}-\frac{q^2}{4\pi\epsilon_0r_l} + \frac{l_o(l_o + 1)\hbar^2}{2mr_l^2},
            \end{equation*}
            and the upper and lower diagonals elements will have values
            \begin{equation*}
                \hat{H}_{l(l\pm1)} = -\frac{\hbar^2}{2ma^2}.
            \end{equation*}
            \item %2b
            Homogenous boundary conditions imply that the corner entries of \(\hat{H}\) will be 0.
            \item %2c
            Code:
            \begin{lstlisting}
clear all;
%physical constants in MKS units

hbar = 1.054e-34;
q = 1.602e-19;
m = 9.110e-31;
epsilon_0 = 8.854e-12;

%generate lattice

N  = 100;                   %number of lattice points
n  = [1:N];                 %lattice points
a  = 0.1e-10;                 %lattice constant
r  = a * n;                 %x-coordinates
t0 = (hbar^2)/(2*m*a^2)/q;  %encapsulating factor
L  = a * (N+1);             %total length of consideration

%set up Hamiltonian matrix 

U = -q^2./(4*pi*epsilon_0.*r) * (1/q); %potential at r in [eV]
main_diag  = diag(2*t0*ones(1,N)+U,0); %create main diagonal matrix
lower_diag = diag(-t0*ones(1,N-1),-1); %create lower diagonal matrix
upper_diag = diag(-t0*ones(1,N-1),+1); %create upper diagonal matrix

H = main_diag + lower_diag + upper_diag; %sum to get Hamiltonian matrix

[eigenvectors,E_diag] = eig(H); %"eigenvectors" is a matrix wherein each column is an eigenvector
                                %"E_diag" is a diagonal matrix where the
                                %corresponding eigenvalues are on the
                                %diagonal.
                                
E_col = diag(E_diag); %folds E_diag into a column vector of eigenvalues

% return eigenvectors for the 1st and 50th eigenvalues

phi_1  = eigenvectors(:,1);
phi_2 = eigenvectors(:,2);

% find the probability densities of position for 1st and 50th eigenvectors

P_1 = phi_1 .* conj(phi_1);
P_2 = phi_2 .* conj(phi_2);

% Plot the probability densities for 1st and 2nd eigenvectors

figure(1); clf; h = plot(r,P_1,'k-');
grid on; set(h,'linewidth',[2.0]); set(gca,'Fontsize',[18]);
xlabel('RADIAL POSITION [m]'); ylabel('PROBABILITY DENSITY [1/m]');
yticks([0.02 0.04 0.06 0.08 0.10 0.12]);
legend('n=1');
axis([0 1e-9 0 0.12]);

figure(2); clf; h = plot(r,P_2,'k-');
grid on; set(h,'linewidth',[2.0]); set(gca,'Fontsize',[18]);
xlabel('RADIAL POSITION [m]'); ylabel('PROBABILITY DENSITY [1/m]');
yticks([0.005 0.01 0.015 0.02 0.025 0.03 0.035 0.04]);
legend('n=2');
axis([0 1e-9 0 0.04]);
            \end{lstlisting}

            \begin{figure}[H]
                \centering
                \begin{subfigure}{0.5\textwidth}
                    \centering
                    \includegraphics[width=\textwidth]{q2c_fig1_1s.png}
                    \caption{1s probability density.}
                \end{subfigure}%
                \begin{subfigure}{0.5\textwidth}
                    \centering
                    \includegraphics[width=\textwidth]{q2c_fig2_2s.png}
                    \caption{2s probability density.}
                \end{subfigure}
            \end{figure}
            \item %2d
            For the 1s level, \(E = -13.4978\) eV.
            \item %2e
            \lipsum[1]
            \item %2f
            In the figure below we can see that the numerical and analytical results agree up to scaling by \(a\). The scale difference is expected,
            as discussed in Problem 1. From (d), we also expect agreement in the curve shapes because the numerical and analytical energies for the 1s
            level are very similiar. We can see that the peak value of the analytic result is very slightly higher than that of the numerical result,
            which corresponds to the analytical result for the energy being slightly greater in magnitude (\(-13.6\) eV versus \(-13.4978\) eV).
            \begin{figure}[H]
                \centering
                \includegraphics[width=\textwidth]{q2f_fig.png}
                \caption{Numerical result (black line) and analytical solution scaled by \(a\) (orange circles).}
            \end{figure}
            


        \end{enumerate}
        
    \section*{Problem 3}

        \begin{enumerate}[align=left,leftmargin=*,labelsep=1em,label=\bfseries(\alph*)]
           
            \item %3a
            \lipsum[1]
            \item %3b
            \begin{enumerate}[align=left,leftmargin=*,labelsep=0em,label=(\roman*)]
                \item %3bi
                \lipsum[1]
                \item %3bii
                \lipsum[1]
                \item %3biii
                \lipsum[1]
                \item %3biv
                \lipsum[1]
                \item %3bv
                \lipsum[1]
                \item %3bvi
                \lipsum[1]
                \item %3bvii
                \lipsum[1]
            \end{enumerate}
            \item %3c
            \lipsum[1]
            \item %3d
            \lipsum[1]

        \end{enumerate}
        
    \section*{Problem 4}

        \begin{enumerate}[align=left,leftmargin=*,labelsep=1em,label=\bfseries(\alph*)]
           
            \item %4a
            \lipsum[1]
            \item %4b
            \begin{enumerate}[align=left,leftmargin=*,labelsep=0em,label=(\roman*)]
                \item %4bi
                \lipsum[1]
                \item %4bii
                \lipsum[1]
                \item %4biii
                \lipsum[1]
            \end{enumerate}

        \end{enumerate}
        
    
\end{document}