%%%% Calculation of complex wavefunctions and their amplitudes for a model
%%%% system of a particle moving towards a rectangular potential barrier.
%%%% The kinetic energy is defined in terms of the de Broglie wavelength,
%%%% the barrier height relative to the kinetic energy and the length of
%%%% the barrier relative to the de Broglie wavelength. Negative x values
%%%% represent the space of the particle moving towards the barrier as well
%%%% as being backscattered, positive x values larger than the barrier
%%%% length represent the space of the particle moving away from the
%%%% barrier after tunneling. For the depiction of the wavefunction, solid
%%%% lines are the real part, dashed lines the imaginary component.
%%%% Figure 2 shows the transmission probability T as a function of the
%%%% relative potential barrier height (upper) and length (lower), if all
%%%% other parameters are kept unchanged. The input values are marked by a
%%%% dotted vertical line.

%%%% General solution of the tunneling problem 
%%%% psi_beforebarrier = A1*exp(i*k1*x)+A2*exp(-i*k1*x)
%%%% psi_inbarrier     = A3*exp(k2*x)+A4*exp(-k2*x)
%%%% psi_afterbarrier  = A5*exp(i*k1*x)

%%%% Coefficients relate to the nomenclature found in Atkins' Physical
%%%% Chemistry, 10th ed., Oxford University Press as follows:
%%%% A1 <-> A, A2 <-> B, A3 <-> C, A4 <-> D, A5 <-> A´
%%%% k1 <-> k, k2 <-> kappa

%%%% Version 1.0 2020/05/14 by B. Corzilius bjoern.corzilius@uni-rostock.de
%%%% distribution according to Creative Commons CC-BY 4.0 license

clear all

%% define symbolic variables

syms A2(A1,k1,k2,L) A3(A1,k1,k2,L) A4(A1,k1,k2,L) A5(A1,k1,k2,L)

A2(A1,k1,k2,L) = (A1*k1^4 - A1*k2^4 + A1*k1*k2^3*2i + A1*k1^3*k2*2i - 2*A1*k1^4*exp(2*L*k2) + 2*A1*k2^4*exp(2*L*k2) + A1*k1^4*exp(4*L*k2) - A1*k2^4*exp(4*L*k2) - A1*k1*k2^3*exp(4*L*k2)*2i - A1*k1^3*k2*exp(4*L*k2)*2i)/(k1^4*exp(4*L*k2) - 2*k2^4*exp(2*L*k2) - 2*k1^4*exp(2*L*k2) + k2^4*exp(4*L*k2) + k1^4 + k2^4 + 2*k1^2*k2^2 + 12*k1^2*k2^2*exp(2*L*k2) + 2*k1^2*k2^2*exp(4*L*k2));
A3(A1,k1,k2,L) = (A1*k1*1i - (k1*(A1*k1^4 - A1*k2^4 + A1*k1*k2^3*2i + A1*k1^3*k2*2i - 2*A1*k1^4*exp(2*L*k2) + 2*A1*k2^4*exp(2*L*k2) + A1*k1^4*exp(4*L*k2) - A1*k2^4*exp(4*L*k2) - A1*k1*k2^3*exp(4*L*k2)*2i - A1*k1^3*k2*exp(4*L*k2)*2i)*1i)/(k1^4*exp(4*L*k2) - 2*k2^4*exp(2*L*k2) - 2*k1^4*exp(2*L*k2) + k2^4*exp(4*L*k2) + k1^4 + k2^4 + 2*k1^2*k2^2 + 12*k1^2*k2^2*exp(2*L*k2) + 2*k1^2*k2^2*exp(4*L*k2)))/(k2 + (k2*(k1^2*exp(2*L*k2) - k2^2*exp(2*L*k2) + k1*k2*exp(2*L*k2)*2i))/(k1^2 + k2^2));
A4(A1,k1,k2,L) = -((k1^2*exp(2*L*k2)*(A1*k1*1i - (k1*(A1*k1^4 - A1*k2^4 + A1*k1*k2^3*2i + A1*k1^3*k2*2i - 2*A1*k1^4*exp(2*L*k2) + 2*A1*k2^4*exp(2*L*k2) + A1*k1^4*exp(4*L*k2) - A1*k2^4*exp(4*L*k2) - A1*k1*k2^3*exp(4*L*k2)*2i - A1*k1^3*k2*exp(4*L*k2)*2i)*1i)/(k1^4*exp(4*L*k2) - 2*k2^4*exp(2*L*k2) - 2*k1^4*exp(2*L*k2) + k2^4*exp(4*L*k2) + k1^4 + k2^4 + 2*k1^2*k2^2 + 12*k1^2*k2^2*exp(2*L*k2) + 2*k1^2*k2^2*exp(4*L*k2))))/(k2 + (k2*(k1^2*exp(2*L*k2) - k2^2*exp(2*L*k2) + k1*k2*exp(2*L*k2)*2i))/(k1^2 + k2^2)) - (k2^2*exp(2*L*k2)*(A1*k1*1i - (k1*(A1*k1^4 - A1*k2^4 + A1*k1*k2^3*2i + A1*k1^3*k2*2i - 2*A1*k1^4*exp(2*L*k2) + 2*A1*k2^4*exp(2*L*k2) + A1*k1^4*exp(4*L*k2) - A1*k2^4*exp(4*L*k2) - A1*k1*k2^3*exp(4*L*k2)*2i - A1*k1^3*k2*exp(4*L*k2)*2i)*1i)/(k1^4*exp(4*L*k2) - 2*k2^4*exp(2*L*k2) - 2*k1^4*exp(2*L*k2) + k2^4*exp(4*L*k2) + k1^4 + k2^4 + 2*k1^2*k2^2 + 12*k1^2*k2^2*exp(2*L*k2) + 2*k1^2*k2^2*exp(4*L*k2))))/(k2 + (k2*(k1^2*exp(2*L*k2) - k2^2*exp(2*L*k2) + k1*k2*exp(2*L*k2)*2i))/(k1^2 + k2^2)) + (k1*k2*exp(2*L*k2)*(A1*k1*1i - (k1*(A1*k1^4 - A1*k2^4 + A1*k1*k2^3*2i + A1*k1^3*k2*2i - 2*A1*k1^4*exp(2*L*k2) + 2*A1*k2^4*exp(2*L*k2) + A1*k1^4*exp(4*L*k2) - A1*k2^4*exp(4*L*k2) - A1*k1*k2^3*exp(4*L*k2)*2i - A1*k1^3*k2*exp(4*L*k2)*2i)*1i)/(k1^4*exp(4*L*k2) - 2*k2^4*exp(2*L*k2) - 2*k1^4*exp(2*L*k2) + k2^4*exp(4*L*k2) + k1^4 + k2^4 + 2*k1^2*k2^2 + 12*k1^2*k2^2*exp(2*L*k2) + 2*k1^2*k2^2*exp(4*L*k2)))*2i)/(k2 + (k2*(k1^2*exp(2*L*k2) - k2^2*exp(2*L*k2) + k1*k2*exp(2*L*k2)*2i))/(k1^2 + k2^2)))/(k1^2 + k2^2);
A5(A1,k1,k2,L) = -(exp(-L*k1*1i)*((k2*exp(-L*k2)*((k1^2*exp(2*L*k2)*(A1*k1*1i - (k1*(A1*k1^4 - A1*k2^4 + A1*k1*k2^3*2i + A1*k1^3*k2*2i - 2*A1*k1^4*exp(2*L*k2) + 2*A1*k2^4*exp(2*L*k2) + A1*k1^4*exp(4*L*k2) - A1*k2^4*exp(4*L*k2) - A1*k1*k2^3*exp(4*L*k2)*2i - A1*k1^3*k2*exp(4*L*k2)*2i)*1i)/(k1^4*exp(4*L*k2) - 2*k2^4*exp(2*L*k2) - 2*k1^4*exp(2*L*k2) + k2^4*exp(4*L*k2) + k1^4 + k2^4 + 2*k1^2*k2^2 + 12*k1^2*k2^2*exp(2*L*k2) + 2*k1^2*k2^2*exp(4*L*k2))))/(k2 + (k2*(k1^2*exp(2*L*k2) - k2^2*exp(2*L*k2) + k1*k2*exp(2*L*k2)*2i))/(k1^2 + k2^2)) - (k2^2*exp(2*L*k2)*(A1*k1*1i - (k1*(A1*k1^4 - A1*k2^4 + A1*k1*k2^3*2i + A1*k1^3*k2*2i - 2*A1*k1^4*exp(2*L*k2) + 2*A1*k2^4*exp(2*L*k2) + A1*k1^4*exp(4*L*k2) - A1*k2^4*exp(4*L*k2) - A1*k1*k2^3*exp(4*L*k2)*2i - A1*k1^3*k2*exp(4*L*k2)*2i)*1i)/(k1^4*exp(4*L*k2) - 2*k2^4*exp(2*L*k2) - 2*k1^4*exp(2*L*k2) + k2^4*exp(4*L*k2) + k1^4 + k2^4 + 2*k1^2*k2^2 + 12*k1^2*k2^2*exp(2*L*k2) + 2*k1^2*k2^2*exp(4*L*k2))))/(k2 + (k2*(k1^2*exp(2*L*k2) - k2^2*exp(2*L*k2) + k1*k2*exp(2*L*k2)*2i))/(k1^2 + k2^2)) + (k1*k2*exp(2*L*k2)*(A1*k1*1i - (k1*(A1*k1^4 - A1*k2^4 + A1*k1*k2^3*2i + A1*k1^3*k2*2i - 2*A1*k1^4*exp(2*L*k2) + 2*A1*k2^4*exp(2*L*k2) + A1*k1^4*exp(4*L*k2) - A1*k2^4*exp(4*L*k2) - A1*k1*k2^3*exp(4*L*k2)*2i - A1*k1^3*k2*exp(4*L*k2)*2i)*1i)/(k1^4*exp(4*L*k2) - 2*k2^4*exp(2*L*k2) - 2*k1^4*exp(2*L*k2) + k2^4*exp(4*L*k2) + k1^4 + k2^4 + 2*k1^2*k2^2 + 12*k1^2*k2^2*exp(2*L*k2) + 2*k1^2*k2^2*exp(4*L*k2)))*2i)/(k2 + (k2*(k1^2*exp(2*L*k2) - k2^2*exp(2*L*k2) + k1*k2*exp(2*L*k2)*2i))/(k1^2 + k2^2))))/(k1^2 + k2^2) + (k2*exp(L*k2)*(A1*k1*1i - (k1*(A1*k1^4 - A1*k2^4 + A1*k1*k2^3*2i + A1*k1^3*k2*2i - 2*A1*k1^4*exp(2*L*k2) + 2*A1*k2^4*exp(2*L*k2) + A1*k1^4*exp(4*L*k2) - A1*k2^4*exp(4*L*k2) - A1*k1*k2^3*exp(4*L*k2)*2i - A1*k1^3*k2*exp(4*L*k2)*2i)*1i)/(k1^4*exp(4*L*k2) - 2*k2^4*exp(2*L*k2) - 2*k1^4*exp(2*L*k2) + k2^4*exp(4*L*k2) + k1^4 + k2^4 + 2*k1^2*k2^2 + 12*k1^2*k2^2*exp(2*L*k2) + 2*k1^2*k2^2*exp(4*L*k2))))/(k2 + (k2*(k1^2*exp(2*L*k2) - k2^2*exp(2*L*k2) + k1*k2*exp(2*L*k2)*2i))/(k1^2 + k2^2)))*1i)/k1;

%% set parameters

a1         = 1;                                     % set scaling factor of a1
wavelength = 1;                                     % de Broglie wavelength of particle
potential  = 2;                                     % height of potential relative to kinetic energy of particle
length     = 0.2 * wavelength;                      % length of barrier relative to de Broglie Wavelength

%% calculate derived parameters

wavenumber = 2*pi/wavelength;                       % wavenumber k of linear motion                     
kinetic    = wavenumber^2/8/pi^2;                   % kinetic energy of particle outside of potential
barrier    = potential * kinetic;                   % potential energy of barrier

if potential == 1
    kappa  = 1e-8;
else
    kappa  = sqrt(2*(barrier-kinetic))*2*pi;        % tunneling decay coefficient
end

%% compute wavefunction amplitudes

a2 = double(vpa(A2(a1,wavenumber,kappa,length),12));
a3 = double(vpa(A3(a1,wavenumber,kappa,length),12));
a4 = double(vpa(A4(a1,wavenumber,kappa,length),12));
a5 = double(vpa(A5(a1,wavenumber,kappa,length),12));

%% compare analytical transmission probablity with amplitude ratio

gamma = abs(a5)^2/abs(a1)^2
T     = (1+(exp(kappa*length)-exp(-kappa*length))^2/(16*(kinetic/barrier)*(1-(kinetic/barrier))))^-1

%% Calculate transmission probability function

dimtrans         = 1000;
height           = logspace(-2,2,dimtrans) * kinetic;
width            = logspace(-2,2,dimtrans) * wavelength;
transh(dimtrans) = 0;
transw(dimtrans) = 0;

for n = 1:dimtrans
    kap  = sqrt(2*(height(n)-kinetic))*2*pi;        
    transh(n) = (1+(exp(kap*length)-exp(-kap*length))^2/(16*(kinetic/height(n))*(1-(kinetic/height(n)))))^-1;
    transw(n) = (1+(exp(kappa*width(n))-exp(-kappa*width(n)))^2/(16*(kinetic/barrier)*(1-(kinetic/barrier))))^-1;
end

figure(2)
clf
subplot(2,1,1)
hold on
gtrans(1,1) = semilogx(height/kinetic,transh)
gtrans(1,2) = scatter(barrier/kinetic,T)
gtrans(1,3) = xline(barrier/kinetic)
set(gca,'XScale','log')
gtrans(1,1).LineWidth=1;
gtrans(1,1).Color='b';
gtrans(1,2).LineWidth=1;
gtrans(1,2).MarkerEdgeColor='r';
gtrans(1,3).Color='r';
gtrans(1,3).LineWidth=1;
gtrans(1,3).LineStyle=':';
xlim([1e-2 1e+2])
ylim([0 1])
xlabel('{\it V} /{\it E_{kin}}')
ylabel('{\it T}')
hold off

subplot(2,1,2)
hold on
gtrans(2,1) = semilogx(width,transw)
gtrans(2,2) = scatter(length,T)
gtrans(2,3) = xline(length)
set(gca,'XScale','log')
gtrans(2,1).LineWidth=1;
gtrans(2,1).Color='b';
gtrans(2,2).LineWidth=1;
gtrans(2,2).MarkerEdgeColor='r';
gtrans(2,3).Color='r';
gtrans(2,3).LineWidth=1;
gtrans(2,3).LineStyle=':';
xlim([1e-2 1e+2])
ylim([0 1])
xlabel('{\it L} /{\it\lambda}')
ylabel('{\it T}')
hold off



%% define space variables

dimx=250;

xleft      = linspace(-2*wavelength,0,dimx);
xbox       = linspace(0,length,dimx);
xright     = linspace(length,length+2*wavelength,dimx);

%% calculate wavefunctions

psileftfwd(dimx) = 0;
psileftbwd(dimx) = 0;
psiboxfwd(dimx)  = 0;
psiboxbwd(dimx)  = 0;
psiright(dimx)   = 0;

for n=1:dimx
    psileftfwd(n) = a1*exp(1i*wavenumber*xleft(n));
    psileftbwd(n) = a2*exp(-1i*wavenumber*xleft(n));
    psiboxfwd(n)  = a3*exp(kappa*xbox(n));
    psiboxbwd(n)  = a4*exp(-kappa*xbox(n));
    psiright(n)   = a5*exp(1i*wavenumber*xright(n));
end

psileft = psileftfwd + psileftbwd;
psibox  = psiboxfwd + psiboxbwd;

%% define potantial function (for visualization only)

xpot(3*dimx) = 0;
pot(3*dimx)  = 0;

for n=1:dimx
    xpot(n)        = xleft(n);
    pot(n)         = 0;
    xpot(dimx+n)   = xbox(n);
    pot(dimx+n)    = barrier/kinetic;
    xpot(2*dimx+n) = xright(n);
    pot(2*dimx+n)  = 0;
end

%% output of amplitudes and probabilities
disp(' ')
disp('List of wavefunction amplitudes:')
disp(['A  =  ',num2str(a1,3)])
disp(['B  =  ',num2str(a2,3)])
disp(['C  =  ',num2str(a3,3)])
disp(['D  =  ',num2str(a4,3)])
disp(['A´ =  ',num2str(a5,3)])
disp(' ')
disp('Transmission probabality:')
disp(['T  =  ',num2str(gamma,3),'  (',num2str(100*gamma,3),' %)'])
disp(' ')


%% plot potential function

figure(1)

subplot(3,1,1)

gpot = plot(xpot,pot);
gpot.Color = 'red';
gpot.LineWidth = 1;
gpot.LineStyle = '-';

gkin = yline(1);
gkin.Color = 'blue';
gkin.LineWidth = 1;
gkin.LineStyle = ':';

xlim([-2*wavelength +2*wavelength+length])
xlabel('{\it x} /{\it\lambda}')
ylabel('{\it V} /{\it E_{kin}}')

%% plot wavefunctions (unnormalized!)

subplot(3,1,2)

g(1,1) = plot(xleft,real(psileft));
g(1,2) = line(xleft,imag(psileft));
g(2,1) = line(xbox,real(psibox));
g(2,2) = line(xbox,imag(psibox));
g(3,1) = line(xright,real(psiright));
g(3,2) = line(xright,imag(psiright));
l(1,1) = xline(0);
l(1,2) = xline(length);

Colors=[0 0.4470 0.7410;0.8500 0.3250 0.0980;0.4660 0.6740 0.1880];
Styles=["-","--"];

for n=1:3
    for m=1:2
        g(n,m).Color     = Colors(n,:);
        g(n,m).LineStyle = Styles(m);
        g(n,m).LineWidth = 1;
    end
end

l(1).Color = 'black';
l(1).LineWidth = 1;
l(1).LineStyle = ':';

l(2).Color = 'black';
l(2).LineWidth = 1;
l(2).LineStyle = ':';

set(gca,'XLim',[-2*wavelength +2*wavelength+length])
set(gca,'YTick',0)
xlabel('{\it x} /{\it\lambda}')
ylabel('{\it\Psi}')


%% plot probability density (unnormalized!)

subplot(3,1,3)

g11 = plot(xleft,abs(psileft).^2);
g11.Color = [0 0.4470 0.7410];
g11.LineWidth = 1;
g11.LineStyle = '-';

g13 = line(xbox,abs(psibox).^2);
g13.Color = [0.8500 0.3250 0.0980];
g13.LineWidth = 1;
g13.LineStyle = '-';

g15 = line(xright,abs(psiright).^2);
g15.Color = [0.4660 0.6740 0.1880];
g15.LineWidth = 1;
g15.LineStyle = '-';

g17= xline(0);
g17.Color = 'black';
g17.LineWidth = 1;
g17.LineStyle = ':';

g18 = xline(length);
g18.Color = 'black';
g18.LineWidth = 1;
g18.LineStyle = ':';

set(gca,'XLim',[-2*wavelength +2*wavelength+length])
set(gca,'YLim',[0 inf],'YTick',0)
xlabel('{\it x} /{\it\lambda}')
ylabel('{\it\Psi}^2')
    