%  FILE: MOS_band.m
%  Band-bending for MOS Capacitor
%     as a funcion of reverse bias

% Constants
eps_0 =   8.85e-14        ; % Units: F/cm
kToq  =   0.0259          ; % Units: V
q     =   1.6e-19         ; % Units: C
cm    =   1.0e4           ; % Units: micron  

% Parameters for Silicon
eps_si =  11.8 * eps_0    ; % Units: F/cm
n_i    =   1.5e10         ; % Units: 1/cm^3
N_a    =   1.0e16         ; % Units: 1/cm^3  <--- INPUT DOPING DENSITY

 
% Parameters for Oxide
eps_ox =   3.9 * eps_0     ; % Units: F/cm
d_ox   =  10.0             ; % Units: nm     <--- INPUT OXIDE THICKNESS
 
 
Phi_F  =  kToq * log( N_a/n_i)   ;
C_ox   =  eps_ox / (d_ox * 1e-7) ;
W_max  =  2 * sqrt( (eps_si*Phi_F) / (q*N_a) ) ;
Q_d    =  q * N_a * W_max  ;
V_T    =  Q_d / C_ox + 2 * Phi_F  ;


V  = linspace(0,1,11);          % Define Range of Depletion Bias Voltages

A = (C_ox/eps_si)*(C_ox/2/q/N_a) ;  a =1/(2*A) ;

V_ox = -a + sqrt(a*a + V/A) ;
V_si = V - V_ox             ;

W = sqrt( 2*eps_si*V/q/N_a) ;

X_si = linspace(0,500,51) ;
X_ox = [-10 0] ;
X_m  = linspace(-50,-10,5) ;
XX   = [X_m X_si];

axis( [-100 400 -1 .1]); hold on
plot( [-100 400] , [ 0 0] ) ; hold on
plot( [0 400] , [-Phi_F -Phi_F], ':') ; hold on
plot( [0 400] , [-2*Phi_F -2*Phi_F], ':') ; hold on
plot( [0 0] , [-1 .1] ) ; plot( [-10 -10] , [-1 .1] ) ;

for iV =1:11

Vm = V(iV)     ;
Vi = V_ox(iV)  ;
Vs = V_si(iV)  ;
WW = W(iV)*cm*1000 + 0.1 ;

xV = linspace(0,WW,15) ;
vV = Vs * (xV/WW - 1) .* (xV/WW - 1) ;


plot( [-80 -10]  , [-Vm -Vm] ) ; hold on
plot( X_ox , [ -Vs-Vi , -Vs] ) ; hold on
plot( xV   , -vV) ; hold on

end