%Chap3Fig14
%Plots the probability and real part of the wave function for an electron
%tunneling through a rectangular barrier of energy V0 and width bx1
%position runs from -5nm to 5nm and the barrier starts at x = 0 

function Tun
clear;
clf;

hbar=1.0545715e-34;  		%Planck's constant (J s)
echarge=1.602176462e-19;	%electron charge (C)
m0=9.10938188e-31;  		%electron mass (kg)

meff=1.0;					%effective electron mass / m0						
mass=m0*meff;				%effective electron mass (kg)
E=0.9;           			%energy of electron(eV)
V0=1.0 ;            		%potential barrier energy(eV)   
bx1=0.5;          			%barrier width(nm)

npoints=1000;              	%defines number of points in plot
scale=npoints/10;
%pre-allocate
x =zeros(1,npoints);
y =zeros(1,npoints);
z1=zeros(1,npoints);z2=zeros(1,npoints);

for j=1:npoints
   x(j)=-5+j/scale;			%position, x (nm)
   psi=Tunnel(x(j),hbar,echarge,mass,E,V0,bx1);%wave function
   y(j)=(abs(psi))^2;		%absolute value of wave function squared
   z1(j)=real(psi);			%real part of wave function
   z2(j)=imag(psi);         %imaginary part of wave function
   end

figure(1);
plot(x,y);
ttl=['Chap3Fig14, E = ',num2str(E),' eV, V_0 = ',num2str(V0),...
    ' eV, L = ',num2str(bx1),' nm, m = ',num2str(meff),'m_0'] ; 
title(ttl)
xlabel('Position, x (nm)');
ylabel('| \psi(x)|^2');

figure(2);
plot(x,z1,'b')
hold on
plot(x,z2,'r');
hold off
title(ttl)
xlabel('Position, x (nm)');
ylabel('Re \psi(x) (blue), Im \psi(x) (red)');

%................. Wave function calculation.............................

function tun=Tunnel(z,hbar,echarge,mass,E,V0,bx1)

L=bx1*1e-9;
x=z*1e-9;         
k=(sqrt(2*mass*E*echarge))/hbar;       %wave vector (outside barrier)
ka=(sqrt(2*mass*(V0-E)*echarge))/hbar; %wave vector (inside barrier)
u=(k/ka);  
F=(4*i*u*exp(-i*k*L))/((exp(-ka*L)*(1+i*u)^2)-(exp(ka*L)*(1-i*u)^2));
C=0.5*(1+i*u)*F*exp(i*k*L)*exp(-ka*L);
D=0.5*(1-i*u)*F*exp(i*k*L)*exp(ka*L);
B=(-i*(1+u^2)*exp(i*k*L)*F*sinh(ka*L))/(2*u);

if (z<=0)
   wf1=exp(i*k*x)+B*exp(-i*k*x);
   tun=wf1;
end   
if (0<z)&(z<=L*1e9)
   wf2=C*exp(ka*x)+D*exp(-ka*x);
   tun=wf2;  
end
if (z>L*1e9)
   wf3=F*exp(i*k*x);
   tun=wf3;
end

