clear all;

% Physical constants in MKS units

hbar = 1.054e-34;
q = 1.602e-19;

% Energy parameters in eV; included are the single-electron charging
% energy U0; the kBT product; the equilibrium Fermi level mu; and the
% energy level cal_E, which is short-form for "calligraphic E"

% Please especially note that ALL ENERGY VARIABLES IN THIS CODE
% ARE IN eV (NOT joules); the equations from class must be adjusted
% accordingly, multiplying or dividing appropriate terms by a factor of q

U0 = 0.25;
kBT = 0.025;
mu = 0;
cal_E = 0.2;

% Capacitance parameters

alpha_G = 0.5;
alpha_D = 0.5;
alpha_S = 1 - alpha_G - alpha_D;

% Energy grid in eV, from -1 eV to 1 eV

NE = 501;
E = linspace(-1,1,NE);
dE = E(2) - E(1);

% Gamma parameters, in eV

gamma_1 = 0.005;
gamma_2 = 0.005;
gamma = gamma_1 + gamma_2;

% Lorentzian density of states, normalized so that its integral is unity

D = (gamma/(2*pi))./((E-cal_E).^2+(gamma/2)^2);
D = D./(dE*sum(D));

% Reference number of electrons in the channel, assumed to be zero in
% this code

N0 = 0;

% Voltage values to consider for the final plots

NV = 101;
VV = linspace(0,1,NV);
dV = VV(2) - VV(1);

% Loop over voltage values and compute number of electrons and current
% for each voltage value in a self-consistent manner

for count = 1:NV

    % Set terminal voltages

    VG = 0;
    VD = VV(count);
    VS = 0;

    % Values of mu1 and mu2; notice that the usual factor of q multiplying
    % the voltages is omitted, because in this code, energy is in eV

    mu1 = mu - VS;
    mu2 = mu - VD;

    % Value of Laplace potential in eV

    UL = - (alpha_G*VG) - (alpha_D*VD) - (alpha_S*VS);

    % Initial value of Poisson part in eV
    UP = 0;

    % Iterate until self-consistent potential is achieved by monitoring
    % the Poisson part (the Laplace part does not change)

    dUP = 1;
    while dUP > 1e-6

        % Compute source and drain Fermi functions

        f1 = 1./(1+exp((E + UL + UP - mu1)./kBT));
        f2 = 1./(1+exp((E + UL + UP - mu2)./kBT));

        % Compute number of channel electrons

        N(count) = dE*sum( ((gamma_1/gamma).*f1 + (gamma_2/gamma).*f2).*D );

        % Newly calculated Poisson part of self-consistent potential

        UPnew = U0*( N(count) - N0 );

        % Change in Poisson part between iterations

        dUP = abs(UP - UPnew);

        %New guess for next iteration, found by adding a fraction of the
        %difference between iterations to the old guess

        UP = UP + 0.1*(UPnew - UP);

    end

    % Compute the current in A after the self-consistent potential
    % has been achieved; notice the extra factor of q preceding the
    % equation, which is needed since the gammas are in eV

    I(count) = q*(q/hbar)*(gamma_1*gamma_2)/(gamma)*dE*sum((f1-f2).*D);

end

% Plotting commands, including lines to modify the linewidth
% and Fontsize, just to make the plots look nicer; you don’t
% need to worry about how these work

figure(1); h = plot(VV,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(VV,I,'k'); grid on;
set(h,'linewidth',[2.0]); set(gca,'Fontsize',[18]);
xlabel('DRAIN VOLTAGE [V]'); ylabel('CURRENT [A]');