% --------------------------------------------------------------------------------
%  MATLAB code used to generate the figures in the book:
%
%    "Quantum Transport: Atom to Transistor," by Supriyo Datta
%    published by Cambridge University Press, May 2005
%      (ISBN-10: 0521631459 | ISBN-13: 9780521631457)
%    http://www.cambridge.org/uk/catalogue/catalogue.asp?isbn=0521631459
%
% THIS FILE FOR: Chapter 7, Figure 7.2.5
%
% --------------------------------------------------------------------------------
% Copyright (c) 2005  Supriyo Datta
% --------------------------------------------------------------------------------
% Modified for SDSMT-NANO 705, P. Ahrenkiel, 2017

clear all

%input
Vg=.20;%Volts
Lz=5;%m
meff_m0=0.3;%electron effective-mass ratio
epsr=5;
Ec_ox=3.;%conduction band edge of oxide
Ec=0.;%conduction band edge of channel

Nox=50;Nc=100;%# of points in each oxide and in channel

%Constants (all MKS, except energy which is in eV)
hc=1240;%eV.nm
m0c2=511e3;%eV
q=1.602e-19;%C
kT=0.025;%eV

%hbar=1.06e-34;
eps0=8.85E-12/q/1e9;%1/(V.nm)
mc2=meff_m0*m0c2;
mu=0;
n0=2*pi*mc2*kT/(hc^2);%1/nm^3

%inputs
a=Lz/Nc;
t0=(hc^2)/2/(2*pi)^2/mc2/(a^2);
e0=q*a/eps0;
Np=Nox+Nc+Nox;
XX=a*((1:1:Np)-(Np+1)/2);
Ec=[Ec_ox*ones(Nox,1);Ec*ones(Nc,1);Ec_ox*ones(Nox,1)];

%Hamiltonian matrix
T=(2*t0*diag(ones(1,Np)))-(t0*diag(ones(1,Np-1),1))-(t0*diag(ones(1,Np-1),-1));

%dielectric constant matrix
D2=epsr*((2*diag(ones(1,Np)))-(diag(ones(1,Np-1),1))-(diag(ones(1,Np-1),-1)));
iD2=inv(D2);

Ubdy=-epsr*[Vg;zeros(Np-2,1);Vg];
%Ubdy=-epsr*[0;zeros(Np-2,1);Vg];;%for asymmetric bias
U0=iD2*Ubdy;

%self-consistent calculation
U1=1e-6*ones(Np,1);UU=U1;change=1;
while change>1e-3
	U1=U1+(0.1*(UU-U1));
	[P,D]=eig(T+diag(Ec)+diag(U1));D=diag(D);
	rho=log(1+exp((mu-D)./kT));rho=P*diag(rho)*P';
	n=2*n0*diag(rho);
	UU=U0+(iD2*e0*n);
	change=max((abs(UU-U1)));
	U=Ec+U1;%self-consistent band profile
end

%electron density in channel per cm2	
ns=1e-4*sum(sum(n.*[zeros(Nox,1);ones(Nc,1);zeros(Nox,1)]))
nn=1e21*n./a;%electron density per cm3
Fn=mu*ones(Nc+Nox+Nox,1);

hold on
h=plot(XX,nn,'g');
%h=plot(XX,Ec,'g');
%h=plot(XX,Ec+U1,'g');
%h=plot(XX,Fn,'g:');
set(h,'linewidth',[2.0])
set(gca,'Fontsize',[14])
xlabel('z (nm)')
%ylabel('E (eV) ')
ylabel('n (cm^-3)')
%axis([0 8 -.5 4])
%axis([0 8 0 15e18])
grid on