%plot for double barrier 
% Minor modifications WM Kaminsky 9/29/2006
%   -- Fixed typo in electron mass.  m_e = 9.1094e-31 (in kg), not
%    9.0195e-31 kg
%   -- Fixed misleading title for second figure.  It now reads "Tunneling
%   through a Double Barrier" not "Tunnneling through a Barrier" because,
%   indeed, it's just the transmission plot in figure 1 with transmission
%   now plotted on a semilog spacing.

clear T R;
global ENG k1a k2a exp_k1a exp_k2a V_0 X AA BB
%
%input of potential in region -a<x<-b and b<x<a
%
V_A= input ( 'What is V_0 in eV?   ');
AA=input( 'What is a in Angstroms?  ');
BB=input('What is b in Angstroms?  ');
mass=input('What is the mass in electron masses?  ');
%
%
%Useful constants so that energies are correct
%
hbar=1.0546e-34;   % Rationalized Planck constant [J*s]
m_e=9.1094e-31;    % Electron rest mass [kg]
JpereV=1.6022e-19; % Number of joules per electron volt
Ang=1e-10;         % Angstrom in meters
% rtwell and the basic program are in reduced units
% as given below.
E_units=2*m_e*mass*BB^2*Ang^2/hbar^2;
%
% Now put V_0 into the reduced units and set range of Energy scale
%
T = [];
R = [];
V_0=V_A*JpereV*E_units;
ERG=[0:2*V_0/5000:2*V_0];
   for  ERG1= ERG,
    COEF=rtwell(ERG1);  % This calls the matrix inversion program
   T=[T,COEF(8)];
   R=[R,COEF(1)];
    end
%
% Transmission and reflection coefficients are now plotted
%
figure
ERG_eV=ERG/(E_units*JpereV);            % This sets Energies back to eV's
subplot(211), plot(ERG_eV,abs(T).^2);
ylabel('T');
title(['Tunneling through a Double Barrier, V= ',num2str(V_A),' eV A=' num2str(AA) ' B= ' num2str(BB)])
subplot(212), plot(ERG_eV,abs(R).^2);
ylabel('R');
xlabel('E in units of eV');
figure
semilogy(ERG_eV, abs(T).^2)
xlabel('E in units of eV');
ylabel('T');
title(['Tunneling through a Double Barrier, V= ',num2str(V_A),' eV'])