%Fig3_17.m
%Essential Electron Transport for Device Physics
%Calculate inelastic transmission through an arbitrary
%potential containing an Einstein phonon at one location.  The temperature
%is 0 K, so that the incident electron cannot initially absorb a phonon.

clear all
clf

FS      = 12;	%label FontSize 18
FSN     = 12;	%number FontSize 16
LW      = 1;	%LineWidth

% Change default axes fonts.
set(0,'DefaultAxesFontName', 'Times New Roman');
set(0,'DefaultAxesFontSize', FSN);

% Change default text fonts.
set(0,'DefaultTextFontname', 'Times New Roman');
set(0,'DefaultTextFontSize', FSN);
%hbar_txt=(['\fontname{MT Extra}h\fontname{Times New Roman}']);
%*******************
% Input parameters *
%*******************
% Phonon parameters
omega = 0.036;  % Phonon e (eV) 0.036 0.01
g0 = 0.15*1e-9; % Static delta barrier e, units of eVnm -0.15 0.15 0.0
g1 = 0.08*1e-9; % units of eVm 0.12 0.08 0.008
g1min=0.02*1e-9;
g1max=0.09*1e-9;
numg1=200;%10
g1array=linspace(g1min,g1max,numg1);

% Electron-phonon coupling constant (eV) 0.07 0.15 0.22
NumPhonons = 8;  % Number of phonons included in calculation 3 8
% Energy Parameters
Emin = pi*1e-6;  % Minimum injection energy (eV) pi*1e-6
Emax = 3.5*omega;% Maximum injection energy (eV)4.5
Enum = 200;      % Number of samples in energy r 1000
Vbias = 0;       % Voltage bias (eV)

% Spatial parameters
xnum = 100;      % Number x-axis points for potential and wave function plots 500

% Barrier parameters
Potential =   [ 0.0 0.0 0.0 0.0 0.0 0.0 0.0];% Potential energy segment (eV)
Length = 1e-9*[ 1.5 1.5 1.5 1.5 1.5 1.5 1.5];% Length of each segment (m)
Phonon =      [ 0.0 0.0 0.0 1.0 0.0 0.0 0.0];% 1 Indicates location of phonon

%*********************************************
% Physical constants and material parameters *
%*********************************************
hbar = 1.05457148e-34;      % Planck's constant (Js)
echarge = 1.60217646e-19;   % Electron charge (C)
m0 = 9.11e-31;              % Electron mass (kg)
meff = 0.07;                % Effective electron mass
mass=m0*meff;

% Include zero phonon
NumPhonons = NumPhonons + 1;

%******************
% Build potential *
%******************
x = linspace(0,sum(Length),xnum);
dx = x(2) - x(1);
xPhonon = zeros(1,xnum);
V = zeros(1,xnum);
BiasSlope = Vbias/(sum(Length) - Length(1) - Length(end));

for j = 1:xnum
    for k = 1:length(Length)
        if (x(j) >= sum(Length(1:k-1))) && (x(j) < sum(Length(1:k)))
            V(j) = Potential(k);
            if (x(j) > Length(1)) && (x(j) <= (sum(Length) - Length(end)))
                V(j) = V(j) - (x(j) - Length(1))*BiasSlope;
            elseif x(j) > (sum(Length) - Length(end))
                V(j) = V(j) - Vbias;
            end
        end        
    end
end

V(end) = Potential(end) - Vbias;

% Find phonon location
PhononLocation = 0;
for j = 1:length(Length)
    if Phonon(j) == 1
        break
    else
        PhononLocation = PhononLocation + Length(j);
    end
end

[z,PhononIndex] = min(abs(x - PhononLocation));
xPhonon(PhononIndex) = 1;

clear z PhononIndex PhononLocation

%******************************************************
% Calculate Transmission and Reflection               *
%******************************************************
% No phonons if g1=0
if g1 == 0
    NumPhonons = 1;
end

%Create injection energy array
E = linspace(Emin,Emax,Enum);

% %loop over values of g1
for nn=1:numg1
    g1=g1array(nn);
%******************************************************
% Initialize arrays
%******************************************************
Transmission = zeros(Enum,1);       % Transmission array
ElasticTrans = zeros(Enum,1);       % Elastic-only transmission array
TransCoeffs = zeros(Enum,NumPhonons + 1); % Transmission coefficients array    
%******************************************************
% Loop over injection energies
for j = 1:Enum
    
    P = eye(2*NumPhonons);  % Reset P
    ElasticP = eye(2);
        
    for k = 1:xnum-1  % Loop through x-axis

        % Create k arrays
        for n = 1:NumPhonons+2
            k1(n) = sqrt(2*mass*echarge*(E(j)-(n-2)*omega - V(k)))/hbar;
            k2(n) = sqrt(2*mass*echarge*(E(j)-(n-2)*omega - V(k+1)))/hbar;
        end

        if xPhonon(k) == 1
            % Create inelastic step matrix
            Pstep = zeros(2*NumPhonons,2*NumPhonons+4);

            for n = 1:NumPhonons      
 AlphaN = 1i*mass*echarge*g1*sqrt(n-1)*ones(1,2)/(k1(n+1)*hbar^2);
  BetaN = 1i*mass*echarge*g1*sqrt(n)*ones(1,2)/(k1(n+1)*hbar^2);

              Pstep(2*n-1,2*n-1:2*n) = AlphaN;
                Pstep(2*n,2*n-1:2*n) = -AlphaN;
            Pstep(2*n-1,2*n+3:2*n+4) = BetaN;
              Pstep(2*n,2*n+3:2*n+4) = -BetaN;

Pstep(2*n-1,2*n+1)= (k1(n+1)+k2(n+1))/(2*k1(n+1))+(1i*mass*echarge*g0)/(hbar^2*k1(n+1));
Pstep(2*n-1,2*n+2)= (k1(n+1)-k2(n+1))/(2*k1(n+1))+(1i*mass*echarge*g0)/(hbar^2*k1(n+1));
Pstep(2*n,2*n+1)  = (k1(n+1)-k2(n+1))/(2*k1(n+1))-(1i*mass*echarge*g0)/(hbar^2*k1(n+1));
Pstep(2*n,2*n+2)  = (k1(n+1)+k2(n+1))/(2*k1(n+1))-(1i*mass*echarge*g0)/(hbar^2*k1(n+1));
            end

            Pstep = Pstep(:,3:2*NumPhonons+2);
        else
            % Create elastic step matrix
            Pstep = zeros(2*NumPhonons);
            for n = 1:NumPhonons
            Pstep(2*n-1,2*n-1) = (k1(n+1) + k2(n+1))/(2*k1(n+1));
                Pstep(2*n,2*n) = Pstep(2*n-1,2*n-1);
              Pstep(2*n,2*n-1) = (k1(n+1) - k2(n+1))/(2*k1(n+1));
              Pstep(2*n-1,2*n) = Pstep(2*n,2*n-1);
            end
        end
        
        ElasticStep = zeros(2);
        % Calculate elastic step matrix
        if xPhonon(k) == 0
            ElasticStep(1,1) = k1(2) + k2(2);
            ElasticStep(2,2) = ElasticStep(1,1);
            ElasticStep(1,2) = k1(2) - k2(2);
            ElasticStep(2,1) = ElasticStep(1,2);
            ElasticStep = ElasticStep/(2*k1(2));
        else
        ElasticStep(1,1) = k1(2) + k2(2) + 1i*2*echarge*mass*g0/(hbar^2);
        ElasticStep(2,2) = ElasticStep(1,1);
        ElasticStep(1,2) = k1(2) - k2(2) + 1i*2*echarge*mass*g0/(hbar^2);
        ElasticStep(2,1) = ElasticStep(1,2);
        ElasticStep = ElasticStep/(2*k1(2));
        end
        
        % Create Propagation Matrix
        Pprop = zeros(2*NumPhonons);
        for n = 1:NumPhonons
  Pprop(2*n-1,2*n-1) = exp(-1i*k2(n+1)*dx);
      Pprop(2*n,2*n) = exp(1i*k2(n+1)*dx);
        end
        
        % Create inelastic propagation amtrix
        ElasticProp = zeros(2);
        ElasticProp(1,1) = exp(-1i*k2(2)*dx);
        ElasticProp(2,2) = exp(1i*k2(2)*dx);

        % Multiply out transfer matrices
        P = P*Pstep*Pprop;
        ElasticP = ElasticP*ElasticStep*ElasticProp;
    end

    % Move bn's into T matrix
    for n = 1:NumPhonons
        P(:,2*n) = zeros(2*NumPhonons,1);
        P(2*n,2*n) = -1;
    end

    % Create A matrix
    A = zeros(2*NumPhonons,1);
    A(1,1) = 1;

    % Solve system
    T = P\A;%note that P\A is faster that inv(P)*A

    % Store largest n that has a propagating mode leaving the barrier
    nMax = 0;
    for n = 1:NumPhonons
        if E(j) > (n-1)*omega + V(end)
            nMax = n;
        end
    end

    % Calculate k values for transmission
    kin = sqrt(2*mass*echarge*(E(j) - V(1)))/hbar;
    kout = zeros(1,nMax);
    for n = 1:nMax
        kout(n) = sqrt(2*mass*echarge*(E(j) - (n-1)*omega - V(end)))/hbar;
    end
    
    % Add transmission coefficients for each of propogating modes, n
    for n = 1:nMax
Transmission(j)=Transmission(j) + abs(kout(n)*T(2*n-1)*conj(T(2*n-1))/kin);
TransCoeffs(j,n+1)=abs(kout(n)*T(2*n-1)*conj(T(2*n-1))/kin);
    end
    TransCoeffs(j,1) = Transmission(j);%Total

    % Solve elastic system
    ElasticTrans(j) = 1/abs(ElasticP(1,1))^2;%No phonon elastic 
end

Trans1(:,nn)=TransCoeffs(:,2);
Trans2(:,nn)=TransCoeffs(:,3);
Trans3(:,nn)=TransCoeffs(:,4);
Trans4(:,nn)=TransCoeffs(:,5);

end

%***************
% Plot results *
%***************
% Create logarithmic transmissions
TransLog = -log(Transmission);

% Set delta magnitude
if max(V) == zeros(1,xnum)
    DeltaMax = 1;
else
    DeltaMax = 1.2*max(V)*sign(max(V));
end
%**************************************************************************
ttl=(['\rmFig3.17b, h\omega=',...
    num2str(omega),' eV, m_{eff}=',num2str(meff),...
    ', g_0=',num2str(g0*1e9),' eV nm, g_1=',...
    num2str(g1*1e9),' eV nm, N_p=',num2str(NumPhonons-1)]);
%**************************************************************************
% Plot transmission by itself
figure(1)
plot(E,ElasticTrans,'k--','Linewidth',LW);
hold on
plot(E,Transmission,'Linewidth',LW)
axis([Emin Emax 0 1.0])
title(ttl)
xlabel(['Energy, E (eV)'])
ylabel('Transmission')
legend('Elastic','Inelastic','Location','Southeast')
grid on

%**************************************************************************
% Plot transmission and potential for comparison
figure(2)
subplot(1,3,1)
plot(Transmission,E/omega,'Linewidth',LW)
title(ttl)
xlabel('Transmission')
ylabel('Normalized energy, E/ $\hbar\omega$', 'Interpreter', 'latex')
axis([0 1 1.2*min(Emin,min(V))/omega max(Emax/omega,max(xPhonon))])
grid on

subplot(1,3,2)
plot(x*1e9,DeltaMax*xPhonon/omega,'r','Linewidth',LW)
hold on
plot(x*1e9,V/omega,'k','Linewidth',LW)
hold off
%title('Potential (delta (red))')
xlabel('x, (nm)')
ylabel('V/ $\hbar\omega$', 'Interpreter', 'latex')
axis([x(1)*1e9 x(end)*1e9 1.2*min(Emin,min(V))/omega max(Emax/omega,max(xPhonon))])
grid on

subplot(1,3,3)
plot(TransLog,E/omega,'Linewidth',LW)
%title('ln T')
xlabel('-ln Transmission')
ylabel('Normalized energy, E/ $\hbar\omega$', 'Interpreter', 'latex')
axis([0 max(TransLog) 1.2*min(Emin,min(V))/omega max(Emax/omega,max(xPhonon))])
grid on 
%**************************************************************************
% Plot Channel Coefficients
figure(3)
plot(E,TransCoeffs,'Linewidth',LW)
title(ttl)
xlabel(['Electron energy, E (eV)'])
ylabel('Transmission coefficients')
axis([Emin Emax 0 1])
grid on

ymax=max(max(-log10(TransCoeffs)));
% Plot Channel Coefficients
figure(4)
plot(E,-log10(TransCoeffs),'Linewidth',LW)
title(ttl)
xlabel(['Electron energy, E (eV)'])
ylabel('-log_{10}(Transmission coefficients)')
axis([Emin Emax 0 ymax])
grid on

% x=linspace(Emin,Emax,Enum);
y=E;
x=linspace(g1min, g1max, numg1);
%**************************************************************************
ttl=(['\rmFig3.17a, h\omega=',...
    num2str(omega),' eV, m_{eff}=',num2str(meff),...
    ', g_0=',num2str(g0*1e9),' eV nm, N_p=',num2str(NumPhonons-1)]);
%**************************************************************************
figure(5)
surf(x*1e9,y,-log10(Trans1));
hold on
surf(x*1e9,y,-log10(Trans2));
surf(x*1e9,y,-log10(Trans3));
surf(x*1e9,y,-log10(Trans4));
%colorbar;
colormap(jet);%jet copper bone gray winter spring
shading interp;
title(ttl);
ylabel(['Electron energy, E (eV)']);
xlabel('g_1 (eV nm)');
zlabel('-log_{10}(Transmission coefficients)');
axis([g1min*1e9 g1max*1e9 0 Emax 0 1.1*max(max(-log10(Trans4)))])
%Rotate the view about the z-axis by XXº.
az = 80; %rotate x-y plane about z-axis by XXº
el = 20;
view(az, el);