Added matlab listings and fixed matlab colouring

This commit is contained in:
David Lenfesty 2021-02-08 21:03:25 -07:00
parent 15b3adfb1d
commit ca5b0a6de1

View File

@ -16,29 +16,35 @@
\definecolor{gray}{rgb}{0.5,0.5,0.5} \definecolor{gray}{rgb}{0.5,0.5,0.5}
\definecolor{mauve}{rgb}{0.58,0,0.82} \definecolor{mauve}{rgb}{0.58,0,0.82}
%\lstset{language=Matlab,%
% %basicstyle=\color{red}, \definecolor{mygreen}{RGB}{28,172,0} % color values Red, Green, Blue
% breaklines=true,% \definecolor{mylilas}{RGB}{170,55,241}
% morekeywords={matlab2tikz},
% keywordstyle=\color{blue},%
% morekeywords=[2]{1}, keywordstyle=[2]{\color{black}}, \lstset{language=Matlab,%
% identifierstyle=\color{black},% %basicstyle=\color{red},
% stringstyle=\color{mylilas}, breaklines=true,%
% commentstyle=\color{mygreen},% morekeywords={matlab2tikz},
% showstringspaces=false,%without this there will be a symbol in the places where there is a space keywordstyle=\color{blue},%
% numbers=left,% morekeywords=[2]{1}, keywordstyle=[2]{\color{black}},
% numbersty_sumle={\tiny \color{black}},% size of th_sume numbers identifierstyle=\color{black},%
% numbersep=9pt, % this defines how far the numbers are from the text stringstyle=\color{mylilas},
% emph=[1]{for,end,break},emphstyle=[1]\color{red}, %some words to emphasise commentstyle=\color{mygreen},%
% %emph=[2]{word1,word2}, emphstyle=[2]{style}, showstringspaces=false,%without this there will be a symbol in the places where there is a space
%} numbers=left,%
\lstset{basicstyle=\small, numberstyle={\tiny \color{black}},% size of th_sume numbers
keywordstyle=\color{mauve}, numbersep=9pt, % this defines how far the numbers are from the text
identifierstyle=\color{dkgreen}, emph=[1]{for,end,break},emphstyle=[1]\color{red}, %some words to emphasise
stringstyle=\color{gray}, %emph=[2]{word1,word2}, emphstyle=[2]{style},
numbers=left,
xleftmargin=5em xleftmargin=5em
} }
%\lstset{basicstyle=\small,
% keywordstyle=\color{mauve},
% identifierstyle=\color{dkgreen},
% stringstyle=\color{gray},
% numbers=left,
% xleftmargin=5em
% }
\title{ECE 456 - Problem Set 1} \title{ECE 456 - Problem Set 1}
\date{2021-02-06} \date{2021-02-06}
@ -230,6 +236,8 @@ ylabel('Current [A]');
\subsection*{(c)} \subsection*{(c)}
Code for this section can be found in appendix A.
\subsubsection*{(i)} \subsubsection*{(i)}
\begin{figure}[H] \begin{figure}[H]
\centering \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. 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)} \subsection*{(b)}
\begin{figure}[H] \begin{figure}[H]
@ -703,6 +782,175 @@ ylabel('Current [A]');
% 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} \end{document}