From feca86b40b3ded108fb9256f60ce737817de89c3 Mon Sep 17 00:00:00 2001
From: pkirwin <pkirwin@ualberta.ca>
Date: Sat, 6 Feb 2021 15:20:27 -0700
Subject: [PATCH] saved changes

---
 PS1/q1b.m | 45 +++++++++++++++++++++++++++++++++++++--------
 1 file changed, 37 insertions(+), 8 deletions(-)

diff --git a/PS1/q1b.m b/PS1/q1b.m
index 86d49ab..496bb32 100644
--- a/PS1/q1b.m
+++ b/PS1/q1b.m
@@ -12,12 +12,11 @@ mstar = 0.25 * 9.11e-31;
 % Single-charge coupling energy (eV)
 U_0 = 0.25;
 % (eV)
-kb_T = 0.025;
+kBT = 0.025;
 % Contact coupling coefficients (eV)
 gamma_1 = 0.0005;
 gamma_2 = gamma_1;
-% Energy level
-E = 0.2;
+gamma_sum = gamma_1 + gamma_2;
 % Capacitive gate coefficient
 a_G = 0.5;
 % Capacitive drain coefficient
@@ -31,6 +30,11 @@ mu = 0;
 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);
 
 
 % Reference no. of electrons in channel
@@ -60,16 +64,41 @@ for n = 1:length(voltages)
     dU_P = 1;
     while dU_P > 1e-6
         % source Fermi function
-        f_1 = 1 / (1 + exp((E + U_L + U_P - mu_1) / kb_T));
+        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) / kb_T));
+        f_2 = 1 / (1 + exp((E + U_L + U_P - mu_2) / kBT));
 
-        N(n) = dE * sum( )
+        % Update channel electrons against potential
+        N(n) = dE * sum( ((gamma_1/gamma_sum) .* f1 + (gamma_2/gamma_sum) .* f2) .* D);
 
-        tmpU_P = U_0 * 
+        % 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
 
+    I(n) = q * (q/hbar) * (gamma_1 * gamma_1 / gamma_sum) * dE * sum((f1-f2).*D);
+
 end
 
 
-%%Plotting commands
\ No newline at end of file
+%%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]');
\ No newline at end of file