clear all;

% Physical constants in MKS units

hbar = 1.054e-34;
q = 1.602e-19;
epsilon_0 = 8.854e-12;
epsilon_r = 4;
mstar = 0.25*9.11e-31;

% MOSFET dimensions in meters

W = 1.0e-6;
L = 10.0e-9;
tox = 1.5e-9;

% Capacitances and capacitance parameters

CG = epsilon_r*epsilon_0*W*L/tox;
CS = 0.05*CG;
CD = 0.05*CG;
CE = CG + CS + CD;
alpha_G = CG / CE;
alpha_D = CD / CE;
alpha_S = (1 - alpha_G - alpha_D);

% Energy parameters in eV

kBT = 0.025;
U0 = q/CE;

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

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

% Escape velocity [m/s] and gamma values [eV]

vR = 1e5;
gamma_1 = hbar*vR/(q*L);
gamma_2 = gamma_1;
gamma = gamma_1 + gamma_2;

% Step-like density of states, stepping from 0 to a finite value
% at an energy of 0 eV

D0 = mstar*q*W*L/(pi*hbar*hbar);
D = D0*[zeros(1,251) ones(1,250)];

% Equilibrium Fermi level in eV

mu = -0.2;

% Reference number of channel electrons

f0 = 1./(1+exp((E - mu)./kBT));
N0 = dE*sum(D.*f0);

% Voltage values to consider for the final plots

NV = 61;
VV = linspace(0,0.6,NV);

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

for VG = [0.25 0.5]

   for count = 1:NV

      % Set terminal voltages

      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 );

         % Update Poisson part of self-consistent potential

         UPnew = U0*( N(count) - N0 );
         dUP = abs(UP - UPnew);
         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 

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

      if (VG == 0.5)

         switch VD
            case 0.0
               figure(2); title('VD = 0.0 V');
               subplot(2,3,1); plot(f1,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/1E4,E,'k-'); axis([-0.1 1.1 -1 1]);
               xlabel('D(E)/1E4'); ylabel('ENERGY  [eV]'); title('VD = 0.0 V');
               subplot(2,3,3); plot(f2,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(f1-f2,E,'--',D/1E4,E,'k-'); axis([-0.1 1.1 -1 1]);
               xlabel('f1(E+U)-f2(E+U), D(E)/1E4'); ylabel('ENERGY  [eV]'); title('VD = 0.0 V');
            case 0.05
               figure(3); title('VD = 0.05 V');
               subplot(2,3,1); plot(f1,E,'k-'); axis([-0.1 1.1 -1 1]);
               xlabel('f1(E+U)'); ylabel('ENERGY  [eV]'); title('VD = 0.05 V');
               subplot(2,3,2); plot(D/1E4,E,'k-'); axis([-0.1 1.1 -1 1]);
               xlabel('D(E)/1E4'); ylabel('ENERGY  [eV]'); title('VD = 0.05 V');
               subplot(2,3,3); plot(f2,E,'k-'); axis([-0.1 1.1 -1 1]);
               xlabel('f2(E+U)'); ylabel('ENERGY  [eV]'); title('VD = 0.05 V');
               subplot(2,3,5); plot(f1-f2,E,'--',D/1E4,E,'k-'); axis([-0.1 1.1 -1 1]);
               xlabel('f1(E+U)-f2(E+U), D(E)/1E4'); ylabel('ENERGY  [eV]'); title('VD = 0.05 V');
            case 0.1
               figure(4); title('VD = 0.1 V');
               subplot(2,3,1); plot(f1,E,'k-'); axis([-0.1 1.1 -1 1]);
               xlabel('f1(E+U)'); ylabel('ENERGY  [eV]'); title('VD = 0.1 V');
               subplot(2,3,2); plot(D/1E4,E,'k-'); axis([-0.1 1.1 -1 1]);
               xlabel('D(E)/1E4'); ylabel('ENERGY  [eV]'); title('VD = 0.1 V');
               subplot(2,3,3); plot(f2,E,'k-'); axis([-0.1 1.1 -1 1]);
               xlabel('f2(E+U)'); ylabel('ENERGY  [eV]'); title('VD = 0.1 V');
               subplot(2,3,5); plot(f1-f2,E,'--',D/1E4,E,'k-'); axis([-0.1 1.1 -1 1]);
               xlabel('f1(E+U)-f2(E+U), D(E)/1E4'); ylabel('ENERGY  [eV]'); title('VD = 0.1 V');
            case 0.2
               figure(5); title('VD = 0.2 V');
               subplot(2,3,1); plot(f1,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/1E4,E,'k-'); axis([-0.1 1.1 -1 1]);
               xlabel('D(E)/1E4'); ylabel('ENERGY  [eV]'); title('VD = 0.2 V');
               subplot(2,3,3); plot(f2,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(f1-f2,E,'--',D/1E4,E,'k-'); axis([-0.1 1.1 -1 1]);
               xlabel('f1(E+U)-f2(E+U), D(E)/1E4'); ylabel('ENERGY  [eV]'); title('VD = 0.2 V');
            case 0.3
               figure(6); title('VD = 0.3 V');
               subplot(2,3,1); plot(f1,E,'k-'); axis([-0.1 1.1 -1 1]);
               xlabel('f1(E+U)'); ylabel('ENERGY  [eV]'); title('VD = 0.3 V');
               subplot(2,3,2); plot(D/1E4,E,'k-'); axis([-0.1 1.1 -1 1]);
               xlabel('D(E)/1E4'); ylabel('ENERGY  [eV]'); title('VD = 0.3 V');
               subplot(2,3,3); plot(f2,E,'k-'); axis([-0.1 1.1 -1 1]);
               xlabel('f2(E+U)'); ylabel('ENERGY  [eV]'); title('VD = 0.3 V');
               subplot(2,3,5); plot(f1-f2,E,'--',D/1E4,E,'k-'); axis([-0.1 1.1 -1 1]);
               xlabel('f1(E+U)-f2(E+U), D(E)/1E4'); ylabel('ENERGY  [eV]'); title('VD = 0.3 V');
         end

      end

   end

   figure(1); plot(VV,ID,'k-'); grid on; 
   xlabel('DRAIN VOLTAGE  [V]'); ylabel('CURRENT  [A]'); hold on;

end