%NANO 705 - Spring 2014
%Homework 1 - Problem 2

clear all

%Constant
q=1.602e-19;%electron charge in C

%Parameters
kT=0.025;%thermal energy
eps0_m_mu0=0;%chemical potential of isolated channel
mu0_m_mus=0.40;%source chemical potential
Up=0.128;%single-electron charging energy
Cfac=0.9;%laplace factor Cfac=1/(1+Cs/Cg)
Neps=10.0;% Number of degenerate states in level eps
r=5.0e12;%transport rate in s^-1
Vg=0;%gate bias

%pre-calculation
eps0_m_mus=mu0_m_mus+eps0_m_mu0;
f0=1./(1+exp(eps0_m_mu0/kT));%initial fermi function for level
fs0=1./(1+exp((eps0_m_mus)/kT));

'------------------------------'
I0=q*r*Neps*(fs0-f0)%display initial current
N0=Neps*f0%display initially filled states

NV=201;VV=linspace(0,2.5,NV);

for iV=1:NV
	Vg=VV(iV);  
    UL=Cfac*Vg;
    U=-UL;%Self-consistent field
    dU=1;%to ensure at least one pass
    alpha=0.1;%convergence factor

    %iterate
    iter=0;
    while dU>1e-10 && iter<10000
        eps_m_mus=eps0_m_mus+U;
        f=1./(exp((eps_m_mus)/kT)+1);
        dN=Neps*(f-fs0);
        Unew=Up*dN-UL;
        dU=abs(U-Unew);
        U=(1-alpha)*U+alpha*Unew;
        iter=iter+1;
    end%while
 
    N(iV)=Neps*(1./(exp((eps_m_mus)/kT)+1))%filled states in channel
    
end
hold on;
p=plot(VV,N);
set(p,'linewidth',[4.0])
set(gca,'Fontsize',[14])
grid on
xlabel('V_g (V)')
ylabel('N')
%legend('eps-mu','eps(0)-mu')
grid on