diff --git a/PS1/doc.tex b/PS1/doc.tex index 2f7e5e6..f22ff84 100644 --- a/PS1/doc.tex +++ b/PS1/doc.tex @@ -16,29 +16,35 @@ \definecolor{gray}{rgb}{0.5,0.5,0.5} \definecolor{mauve}{rgb}{0.58,0,0.82} -%\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,% -% numbersty_sumle={\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}, -%} -\lstset{basicstyle=\small, - keywordstyle=\color{mauve}, - identifierstyle=\color{dkgreen}, - stringstyle=\color{gray}, - numbers=left, - xleftmargin=5em - } + +\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=5em +} +%\lstset{basicstyle=\small, +% keywordstyle=\color{mauve}, +% identifierstyle=\color{dkgreen}, +% stringstyle=\color{gray}, +% numbers=left, +% xleftmargin=5em +% } \title{ECE 456 - Problem Set 1} \date{2021-02-06} @@ -230,6 +236,8 @@ ylabel('Current [A]'); \subsection*{(c)} + Code for this section can be found in appendix A. + \subsubsection*{(i)} \begin{figure}[H] \centering @@ -443,6 +451,77 @@ ylabel('Current [A]'); Below is our code. Note that some variable names are different from those in the example code. + + \begin{lstlisting}[language=Matlab] +% Thermo-electric current + +% Physical constants +hbar = 1.054e-34; +q = 1.602e-19; + +% Parameters (eV) +kBT_1 = 0.025; +kBT_2 = 0.026; +mu = 0; +gamma_1 = 0.005; +gamma_2 = gamma_1; +gamma_sum = gamma_1 + gamma_2; + +% Channel energy levels, varying between -0.25eV and 0.25eV +epsilon = linspace(-0.25, 0.25, 101); +depsilon = epsilon(2) - epsilon(1); + +% Energy grid +E = linspace(-1, 1, 501); +dE = E(2) - E(1); + +% Contact fermi functions +f_1 = 1 ./ (1 + exp((E - mu)./kBT_1)); +f_2 = 1 ./ (1 + exp((E - mu)./kBT_2)); + +% Iterate through channel energy levels +for n = 1:length(epsilon) + % Compute energy level density functions - integral normalized to unity + D = (gamma_sum./(2*pi))./((E-epsilon(n)).^2+((gamma_sum./2).^2)); + D = D./(dE*sum(D)); + + % Compute number of channel electrons + + N(n) = dE*sum( ((gamma_1./gamma_sum).*f_1 + (gamma_2./gamma_sum).*f_2).*D ); + + % Compute the current in Amps; factor of q to resolve units + + I(n) = q*(q/hbar)*dE*sum((f_1 - f_2).*D.*gamma_1.*gamma_2./gamma_sum); + + %plot f_1 - f_2 and D/2 + + if (abs(epsilon(n) + 0.05) <= depsilon / 2) & (epsilon(n) <= 0) + figure(3); + h = plot(f_1-f_2, E, 'x', D/2500, E, 'k-'); + set(gca, 'Fontsize', [18]); + axis([-0.01 0.02 -1 1]); + xlabel('f1(E) - f2(E), D(E)/2500'); + ylabel('ENERGY [eV]'); + legend('f1-f2', 'D(E)/2500'); + title('CHANNEL LEVEL = -0.05 eV'); + elseif (abs(epsilon(n)) <= depsilon / 2) & (epsilon(n) <= 0) + figure(4); + h = plot(f_1-f_2, E, 'x', D/2500, E, 'k-'); + set(gca, 'Fontsize', [18]); + axis([-0.01 0.02 -1 1]); + xlabel('f1(E) - f2(E), D(E)/2500'); + ylabel('ENERGY [eV]'); + legend('f1-f2', 'D(E)/2500'); + title('CHANNEL LEVEL = 0 eV'); + end + +end + + +% Final plots +figure(1); + \end{lstlisting} + \subsection*{(b)} \begin{figure}[H] @@ -700,9 +779,178 @@ ylabel('Current [A]'); \end{equation*} % We conclude there are 3 levels. - -% TODO appendix for Part C code + +\appendix + + \newpage + \section{Question 1b Code} + + \begin{lstlisting}[language=Matlab] +clear all; + +%% Constants + +% Physical constants +hbar = 1.052e-34; +q = 1.602e-19; +%epsilon_0 = 8.854e-12; +%epsilon_r = 4; +%mstar = 0.25 * 9.11e-31; + +% Single-charge coupling energy (eV) +U_0 = 0.25; +% (eV) +kBT = 0.025; +% Contact coupling coefficients (eV) +gamma_1 = 0.005; +gamma_2 = gamma_1; +gamma_sum = gamma_1 + gamma_2; +% Capacitive gate coefficient +a_G = 0.5; +% Capacitive drain coefficient +a_D = 0.5; +a_S = 1 - a_G - a_D; + +% Central energy level +mu = 0; + +% Energy grid, from -1eV to 1eV +NE = 501; +E = linspace(-1, 1, NE); +dE = E(2) - E(1); +% TODO name this better +cal_E = 0.2; + +% Lorentzian density of states, normalized so the integral is 1 +D = (gamma_sum / (2*pi)) ./ ( (E-cal_E).^2 + (gamma_sum/2).^2 ); +D = D ./ (dE*sum(D)); + +% Reference no. of electrons in channel +N_0 = 0; + +voltages = linspace(0, 1, 101); +dV = voltages(2) - voltages(1); + +% Terminal Voltages +V_G = 0; +V_S = 0; + +for n = 1:length(voltages) + % Set varying drain voltage + V_D = voltages(n); + + % Shifted energy levels of the contacts + mu_1 = mu - V_S; + mu_2 = mu - V_D; + + % Laplace potential, does not change as solution is found (eV) + % q is factored out here, we are working in eV + U_L = - (a_G*V_G) - (a_D*V_D) - (a_S*V_S); + + % Poisson potential must change, assume 0 initially (eV) + U_P = 0; + + % Assume large rate of change + dU_P = 1; + + % Run until we get close enough to the answer + while dU_P > 1e-6 + % source Fermi function + f_1 = 1 ./ (1 + exp((E + U_L + U_P - mu_1) ./ kBT)); + % drain Fermi function + f_2 = 1 ./ (1 + exp((E + U_L + U_P - mu_2) ./ kBT)); + + % Update channel electrons against potential + N(n) = dE * sum( ((gamma_1/gamma_sum) .* f_1 + (gamma_2/gamma_sum) .* f_2) .* D); + + % Re-update Poisson portion of potential + tmpU_P = U_0 * ( N(n) - N_0); + dU_P = abs(U_P - tmpU_P); + + % Unsure why U_P is updated incrementally, perhaps to avoid oscillations? + %U_P = tmpU_P; + U_P = U_P + 0.1 * (tmpU_P - U_P); + end + + % Calculate current based on solved potential. + % Note: f1 is dependent on changes in U but has been updated prior in the loop + I(n) = q * (q/hbar) * (gamma_1 * gamma_1 / gamma_sum) * dE * sum((f_1-f_2).*D); + + if (abs(V_D-0.0) <= dV/2) + figure(3); title('VD = 0.0 V'); + subplot(2,3,1); plot(f_1,E,'k-'); axis([-0.1 1.1 -1 1]); + xlabel('f1(E+U)'); ylabel('ENERGY [eV]'); title('VD = 0.0 V'); + subplot(2,3,2); plot(D/100,E,'k-'); axis([-0.1 1.1 -1 1]); + xlabel('D(E)/100'); ylabel('ENERGY [eV]'); title('VD = 0.0 V'); + subplot(2,3,3); plot(f_2,E,'k-'); axis([-0.1 1.1 -1 1]); + xlabel('f2(E+U)'); ylabel('ENERGY [eV]'); title('VD = 0.0 V'); + subplot(2,3,5); plot(f_1-f_2,E,'--',D/100,E,'k-'); axis([-0.1 1.1 -1 1]); + xlabel('f1(E+U)-f2(E+U), D(E)/100'); ylabel('ENERGY [eV]'); title('VD = 0.0 V'); + elseif (abs(V_D-0.2) <= dV/2) + figure(4); title('VD = 0.2 V'); + subplot(2,3,1); plot(f_1,E,'k-'); axis([-0.1 1.1 -1 1]); + xlabel('f1(E+U)'); ylabel('ENERGY [eV]'); title('VD = 0.2 V'); + subplot(2,3,2); plot(D/100,E,'k-'); axis([-0.1 1.1 -1 1]); + xlabel('D(E)/100'); ylabel('ENERGY [eV]'); title('VD = 0.2 V'); + subplot(2,3,3); plot(f_2,E,'k-'); axis([-0.1 1.1 -1 1]); + xlabel('f2(E+U)'); ylabel('ENERGY [eV]'); title('VD = 0.2 V'); + subplot(2,3,5); plot(f_1-f_2,E,'--',D/100,E,'k-'); axis([-0.1 1.1 -1 1]); + xlabel('f1(E+U)-f2(E+U), D(E)/100'); ylabel('ENERGY [eV]'); title('VD = 0.2 V'); + elseif (abs(V_D-0.5) <= dV/2) + figure(5); title('VD = 0.5 V'); + subplot(2,3,1); plot(f_1,E,'k-'); axis([-0.1 1.1 -1 1]); + xlabel('f1(E+U)'); ylabel('ENERGY [eV]'); title('VD = 0.5 V'); + subplot(2,3,2); plot(D/100,E,'k-'); axis([-0.1 1.1 -1 1]); + xlabel('D(E)/100'); ylabel('ENERGY [eV]'); title('VD = 0.5 V'); + subplot(2,3,3); plot(f_2,E,'k-'); axis([-0.1 1.1 -1 1]); + xlabel('f2(E+U)'); ylabel('ENERGY [eV]'); title('VD = 0.5 V'); + subplot(2,3,5); plot(f_1-f_2,E,'--',D/100,E,'k-'); axis([-0.1 1.1 -1 1]); + xlabel('f1(E+U)-f2(E+U), D(E)/100'); ylabel('ENERGY [eV]'); title('VD = 0.5 V'); + elseif (abs(V_D-0.8) <= dV/2) + figure(6); title('VD = 0.8 V'); + subplot(2,3,1); plot(f_1,E,'k-'); axis([-0.1 1.1 -1 1]); + xlabel('f1(E+U)'); ylabel('ENERGY [eV]'); title('VD = 0.8 V'); + subplot(2,3,2); plot(D/100,E,'k-'); axis([-0.1 1.1 -1 1]); + xlabel('D(E)/100'); ylabel('ENERGY [eV]'); title('VD = 0.8 V'); + subplot(2,3,3); plot(f_2,E,'k-'); axis([-0.1 1.1 -1 1]); + xlabel('f2(E+U)'); ylabel('ENERGY [eV]'); title('VD = 0.8 V'); + subplot(2,3,5); plot(f_1-f_2,E,'--',D/100,E,'k-'); axis([-0.1 1.1 -1 1]); + xlabel('f1(E+U)-f2(E+U), D(E)/100'); ylabel('ENERGY [eV]'); title('VD = 0.8 V'); + elseif (abs(V_D-1.0) <= dV/2) + figure(7); title('VD = 1.0 V'); + subplot(2,3,1); plot(f_1,E,'k-'); axis([-0.1 1.1 -1 1]); + xlabel('f1(E+U)'); ylabel('ENERGY [eV]'); title('VD = 1.0 V'); + subplot(2,3,2); plot(D/100,E,'k-'); axis([-0.1 1.1 -1 1]); + xlabel('D(E)/100'); ylabel('ENERGY [eV]'); title('VD = 1.0 V'); + subplot(2,3,3); plot(f_2,E,'k-'); axis([-0.1 1.1 -1 1]); + xlabel('f2(E+U)'); ylabel('ENERGY [eV]'); title('VD = 1.0 V'); + subplot(2,3,5); plot(f_1-f_2,E,'--',D/100,E,'k-'); axis([-0.1 1.1 -1 1]); + xlabel('f1(E+U)-f2(E+U), D(E)/100'); ylabel('ENERGY [eV]'); title('VD = 1.0 V'); + end + +end + + +%%Plotting commands + +figure(1); +h = plot(voltages, N,'k'); +grid on; +set(h,'linewidth',[2.0]); +set(gca,'Fontsize',[18]); +xlabel('Drain voltage [V]'); +ylabel('Number of electrons'); + +figure(2); +h = plot(voltages, I,'k'); +grid on; +set(h,'linewidth',[2.0]); +set(gca,'Fontsize',[18]); +xlabel('Drain voltage [V]'); +ylabel('Current [A]'); + \end{lstlisting} + \end{document} \ No newline at end of file