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)
kb_T = 0.025;
% Contact coupling coefficients (eV)
gamma_1 = 0.0005;
gamma_2 = gamma_1;
% Energy level
E = 0.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);


% Reference no. of electrons in channel
N_0 = 0;

fermi(-0.25, -0.2, kb_T)

voltages = linspace(0, 1, 101);

% Terminal Voltages
V_G = 0;
V_S = 0;

for n = 1:length(voltages)
    % Set varying drain voltage
    V_D = voltages(n);

    mu_1 = mu - V_S;
    mu_2 = mu - V_D;

    % Laplace potential, does not change as solution is found (eV)
    U_L = -q * ((C_S*V_S + C_G*V_G + C_D*V_D) / C_E);

    % Poisson potential must change, assume 0 initially (eV)
    U_P = 0;

    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));
        % drain Fermi function
        f_2 = 1 / (1 + exp((E + U_L + U_P - mu_2) / kb_T));

        N(n) = dE * sum( )

        tmpU_P = U_0 * 
    end

end


%%Plotting commands