%Fig3_8.m %Essential Electron Transport for Device Physics %resonant tunneling through 3 barriers clear clf; FS = 12; %label fontsize 18 FSN = 12; %number fontsize 16 LW = 1; %linewidth % Change default axes fonts. set(0,'DefaultAxesFontName', 'Times'); set(0,'DefaultAxesFontSize', FSN); % Change default text fonts. set(0,'DefaultTextFontname', 'Times'); set(0,'DefaultTextFontSize', FSN); hbar_txt=['\fontname{MT Extra}h\fontname{Arial}']; t=cputime; nbarrier=3; %number of potential barriers bx=0.4e-9; %barrier width (m) wx=0.6e-9; %well width (m) V0=0.5; %barrier energy (eV) N=140+1; %number of samples of potential for j=1:1:20 %set up position and potential array dL(j)=wx/1000; V(j)=0; dL(j+20)=bx/20; V(j+20)=V0; %barrier 1 dL(j+40)=wx/20; V(j+40)=0; dL(j+60)=bx/20; V(j+60)=V0; %barrier 2 dL(j+80)=wx/20; V(j+80)=0; dL(j+100)=bx/20; V(j+100)=V0; %barrier 3 dL(j+120)=wx/20; V(j+120)=0; end dL(N)=wx/20; V(N)=0; xp(1)=dL(1); for j=2:N;xp(j)=xp(j-1)+dL(j);end; %set up x-plot axis xpmax=max(xp); Emax=0.6; %maximum particle energy (eV) npoints=141; %number of points in energy plot hbar=1.0545715e-34; %Planck's constant (J s) eye=complex(0.,1.); %square root of -1 m0=9.109382e-31; %bare electron mass (kg) meff=1.0; %effective electron mass / m0 m=meff*m0; %effective electron mass (kg) echarge=1.6021764e-19; %electron charge (C) % %************************************************************************** %Calculate electron wave function %************************************************************************** %loop energy E1=[Emax/npoints:Emax/npoints:Emax]; %set up energy array of npoints psi2=zeros(npoints,N);%initialize psi2 to zero C=zeros(1,N);D=zeros(1,N); for j=1:npoints bigP=[1,0;0,1]; %default value of matrix bigP for i=1:N k(i)=sqrt(2*echarge*m*(E1(j)-V(i)))/hbar;%wave vector at energy E1 end for n=1:(N-1) %multiply out propagation matrix a=(1+k(n+1)/k(n)); b=(1-k(n+1)/k(n)); c=(exp(-eye*k(n)*dL(n))); d=(exp(eye*k(n)*dL(n))); p(1,1)=a*c; p(1,2)=b*c; p(2,1)=b*d; p(2,2)=a*d; bigP=0.5*bigP*p; end C(1)=1/bigP(1,1); %transmission amplitude B(1)=bigP(2,1)*C(1); %reflection amplitude A(1)=1; %incident amplitude bigP=[1,0;0,1]; %default value of matrix bigP for n=1:(N-1) %multiply out propagation matrix a=(1+k(n+1)/k(n)); b=(1-k(n+1)/k(n)); c=(exp(-eye*k(n)*dL(n))); d=(exp(eye*k(n)*dL(n))); p(1,1)=a*c; p(1,2)=b*c; p(2,1)=b*d; p(2,2)=a*d; bigP=0.5*bigP*p; Prop=inv(bigP)*[A(1);B(1)]; C(n)=Prop(1); D(n)=Prop(2); end for n=1:(N-1) psi(n)=(C(n)+D(n)); %calculate wave function end psi(N)=psi(N-1); %dummy value for N-th value of wave function psi2(j,:)=conj(psi).*psi; %calculate modulus wave function squared end psi2max=max(max(psi2)); %find maximum value of psi2 %************************************************************************** %Plot modulus of wave function squared as function of position and energy %************************************************************************** ttl1=['\rmFig3.8 (',num2str(nbarrier),'-barriers), m^*_e='... ,num2str(meff),'\timesm_0']; figure(1)%surface plot with lighting surfl(xp,E1,psi2); shading interp; colormap(copper); axis square; xlabel('Position, x (m)'),ylabel('Energy, E (eV)'),zlabel('|\psi(x)|^2'); title(ttl1); axis([ 0 xpmax 0 Emax 0 psi2max]); ttl2=['\rmFig3.8, -ln(|\psi|^2) contour (',num2str(nbarrier),... '-barriers), m^*_e= ',num2str(meff),'\timesm_0']; figure(2)%contour plot of -ln(psi2) contour(xp,E1,-log(psi2),40); axis square; xlabel('Position, x (m)'),ylabel('Energy, E (eV)'),zlabel('|\psi(x)|^2'); title(ttl2); axis([ 0 xpmax 0 Emax]); e=cputime-t; e