Fowler Nordheim Plot Simulations (F-N Plot)

Last few days I was trying to do a simulation on Fowler Nordheim tunneling. It has different name. Mostly it is known as F-N plot or J/E^2 vs. 1/E plot. Here is the matlab coding for this simulation.



%% Electric Field
% $$ V_t=V_{fb}+2(\psi_b)+ $$
Ub=.42;
Vfb=4.1-(4.05+.56+Ub);        %Flate band voltage   4.1-(4.05+1.12/2+.6)        
Vg=0:0.1:10;                %For positive gate voltage
%Vg(202:402)=0:-0.1:-20;
%Vgneg=-20:0.1:0;              %For negative gate voltage
Vsub=-1;                     %For positive sunstrate voltage
%Vsub(202:402)=-1;
                 %Vsubng=-1;                    %For negative sunstrate voltage
doxfs=7.8e-7;                 % Dielectric thickness
E=(Vg-Vfb-Vsub)./doxfs;
%% Current Density J
%A constant
mo=9.1093e-31;
mox=0.5*mo;
%moxmo=mox/mo
q=1.602e-19; %C
Si_SiO=3.2; %eV
h=4.135e-15;   %eV.s
%A=1.54e-6;
%B=6.87e7;
%A=(mo/mox)*(q^3)/((8*pi*h*Si_SiO)*q^2)
%B Constant
%B=(8*pi*(sqrt(2*mox))*Si_SiO^(1.5))*q/(3*q*h)  %V/m
%Alternative A and B calculation          
% Ref: On tunneling in metal-oxide-silicon structures by z. A.  Weinberg
%IBM Thomas J.  Watson Research Center,  Yorktown Heights,  New York 10598
A=(mo/mox)*(1/Si_SiO)*1.54e-6         %A/{V}^2
B=sqrt(mox/mo)*(Si_SiO^(3/2))*6.83e7  %V/cm
% Current Density Finial Calculation
y=(A.*exp(-B./(E)));
J=y.*(E.^2);
I=J*1.05e-7;
%J=(A*E.^2).*exp(-B./E);
%y=log(J./E.^2);
%x=(1e-7.*E).^-1;
x=1./(E);
%y=log(A)-(B./E);
%y=abs(yy);
xx=transpose(x);
yy=transpose(y);
plot(x,y)
ylabel('J/{E}^2,(Amp/V)')
xlabel('1/{E},({10}^{-7}cm/MV)')
title('FN Plot')