%Fig3_16b.m
%Essential Electron Transport for Device Physics
%Calculate inelastic transmission through an arbitrary
%potential containing an Einstein phonon at one location.  The temperature
%is 0K, so that the incident electron cannot 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.01, g1=0.07, meff=0.07, Lb=2, LW=2.95+0.01+2.95, V0=0.3, Emax=30
%omega=0.01, g1=0.07, meff=0.07, Lb=1.5, LW=2.95+0.01+2.95, V0=0.3, Emax=30
%%omega=0.01, g1=0.22, meff=1.0, Lb=0.4, LW=0.295+0.01+0.295, V0=1, Emax=150
omega = 0.036;   % Phonon Energy (eV) 0.036 0.01
g0 = -0.15*1e-9; % Static Delta Barrier Energy, units of eVnm -0.15 0.15 0.0
g1 = 0.08*1e-9;  % units of eVm 0.12 0.08 0.008
%g1c is 0.14 eV nm in GaAs
% Electron-Phonon Coupling Constant (eV) 0.07 0.15 0.22
NumPhonons = 3;  % Number of Phonons Included in Calculation 8

% Energy Parameters
Eplot = 1.3774;  % Energy Value for Plotting of Wavefunction
EMin = pi*1e-5;  % Minimum Injection Energy (eV)
EMax = 3.5*omega;% Maximum Injection Energy (eV)4.5 8.5
NumE = 500;      % Number of Samples in Energy Range 1000
Vbias = 0;       % Voltage Bias (eV)

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

% Material Parameters

% Barrier Parameters
Potential =   [ 0.0 0.0 0.0 0.0 0.0 0.0 0.0];  % Potential Energy of Each 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 Which Barrier Phonon is Located at

%*************************
% Constants and Matrices *
%*************************
% Physical Constants
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;
%convert g to units of eV
%gconv=sqrt(meff*m0*echarge/(2*hbar));
% gconv=sqrt(2*meff*m0*echarge)/hbar;
% g0=g0*gconv;
% g1=g1*gconv;

% Include Zero Phonon
NumPhonons = NumPhonons + 1;

% Initialize Arrays
Transmission = zeros(NumE,1);       % Transmission Array
ElasticTrans = zeros(NumE,1);       % Elastic Only Transmission Array
TransCoeffs = zeros(NumE,NumPhonons + 1); % Transmission Coefficients Array

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

for j = 1:Numx
    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 Due to System *
%******************************************************
% No Phonons if g1=0
if g1 == 0
    NumPhonons = 1;
end

%Create Injection Energy Array
E = linspace(EMin,EMax,NumE);

% Find E value closest to Eplot
Ediff = (E - Eplot).^2;
[z,EplotIndex] = min(Ediff);
clear z

% Loop Over Injection Energies
for j = 1:NumE
    
    P = eye(2*NumPhonons);  % Reset P
    ElasticP = eye(2);
        
    for k = 1:Numx-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 Matrix
        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 Propogating Modes
    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);

    % Save data for wavefunction plots
    if j == EplotIndex
        TPlot = T;
    end
    
    % Solve elastic system
    ElasticTrans(j) = 1/abs(ElasticP(1,1))^2;
end

%**************************
% Calculate Wavefunctions *
%**************************
% Eliminate b values from T matrix
for n = 1:NumPhonons
    TPlot(2*n) = 0;
end

% Initialize arrays and matrices for wavefunction calculation
Psi = zeros(2*NumPhonons,Numx);
Psi(:,end) = TPlot;
PsiPlot = zeros(NumPhonons,Numx);

% Back Calculate Wavefunction Values
for j = 1:Numx - 1
    k = Numx - j;
    
        % Create k arrays
    for n = 1:NumPhonons+2
        k1(n) = sqrt(2*mass*echarge*(E(EplotIndex)-(n-2)*omega - V(k)))/hbar;
        k2(n) = sqrt(2*mass*echarge*(E(EplotIndex)-(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*echarge*g1*sqrt(n-1)*ones(1,2)/(k1(n+1));
            BetaN = 1i*echarge*g1*sqrt(n)*ones(1,2)/(k1(n+1));

          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

    % 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

    Psi(:,k) = Pstep*Pprop*Psi(:,k+1);
end

% Build and Normalize Wavefunctions
for k = 1:Numx
    for n = 1:NumPhonons
        PsiPlot(n,k) = Psi(2*n-1,k) + Psi(2*n,k);
    end
end

% Normalize Wavefunctions and Find Zero Values for Plots
WaveMax = zeros(1,NumPhonons);
for n = 1:NumPhonons
    if max(PsiPlot(n,:)) ~= 0
        WaveMax(n) = max(PsiPlot(n,:));
        PsiPlot(n,:) = 0.9*PsiPlot(n,:)*omega/(2*WaveMax(n));
    end
    PsiZero(n) = E(EplotIndex) - (n-1)*omega;
end

PsiMagPlot = 3.5*PsiPlot.*conj(PsiPlot);
    
%***************
% Plot Results *
%***************
% Create logarithmic transmissions
TransLog = -log(Transmission);

% Calculate Transmission Derivatives
TransDiff = diff(Transmission);
Trans2Diff = diff(TransDiff);

% Set delta magnitude
if max(V) == zeros(1,Numx)
    DeltaMax = 1;
else
    DeltaMax = 1.2*max(V)*sign(max(V));
end
ttl=(['\rmFig3.16b, 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)
clf
%subplot(1,2,1)
plot(E,ElasticTrans,'k--','Linewidth',LW);
hold on
plot(E,Transmission,'Linewidth',LW)
%axis([EMin EMax min(Transmission) 1.002])
axis([EMin EMax 0 1.0])
title(ttl)
xlabel(['Energy, E (eV)'])
ylabel('Transmission')
legend('Elastic','Inelastic','Location','Southeast')
grid on

%**************************************************************************
% Plot rotated transmission and potential for comparison
figure(2)
clf
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)
clf
plot(E,TransCoeffs,'Linewidth',LW)
title(ttl)
xlabel(['Energy, E (eV)'])
ylabel('Transmission coefficients')
axis([EMin EMax 0 1])
grid on

% Make Normalization Text
if NumPhonons == 1
Normalizations = [];
else
Normalizations = ['Multipliers for wavefunctions are '];
Normalizations2= ['Multipliers for densities are '];
    for n = 1:NumPhonons
Normalizations =[Normalizations num2str(abs(WaveMax(1)/WaveMax(n))) ', '];
Normalizations2=[Normalizations2 num2str(abs(WaveMax(1)/WaveMax(n))^2) ', '];
    end
end
ymax=max(max(-log10(TransCoeffs)));
% Plot Channel Coefficients
figure(4)
clf
plot(E,-log10(TransCoeffs),'Linewidth',LW)
title(ttl)
xlabel(['Energy, E (eV)'])
ylabel('-log_{10}(Transmission coefficients)')
axis([EMin EMax 0 ymax])
grid on



% %**************************************************************************
% return
% %**************************************************************************

% Plot Transmission Derivatives
figure(5)
clf
subplot(1,2,1)
plot(E(1:end-1)/omega,TransDiff,'Linewidth',2)
title('\rmFirst derivative of transmission')
xlabel('Normalized energy (E/\Omega)')
ylabel('dT/dE')
xlim([E(1)/omega E(end-1)/omega])
grid on
subplot(1,2,2)
plot(E(1:end-2)/omega,Trans2Diff,'Linewidth',2)
title('\rmSecond derivative of transmission')
xlabel('Normalized energy (E/\Omega)')
ylabel('d^2T/dE^2')
xlim([E(1)/omega E(end-1)/omega])
grid on

return

% Plot Wavefunctions
figure(5)
clf
subplot(1,3,1)
hold on
for n = 1:NumPhonons
    plot(x*1e9,real(PsiPlot(n,:) + PsiZero(n)),'b')
    plot(x*1e9,imag(PsiPlot(n,:) + i*PsiZero(n)),'r')
    plot(x*1e9,PsiZero(n)*ones(1,Numx),'g')
end
hold off
title({['Wavefunctions for inelastically scattered electron with injection energy of ',num2str(E(EplotIndex)),' eV'];['Wavefunctions offset so that the zero level is E_{inj} - n\Omega and normalized to prevent overlap'];Normalizations})
xlabel('x (nm)')
ylabel('Wavefunction, \psi (m^{-1/2})')
axis([x(1)*1e9 x(end)*1e9 min((E(EplotIndex)-(NumPhonons-0.5)*omega),0) 1.2*max((E(EplotIndex)+omega/2),max(V))])

subplot(1,3,2)
plot(x*1e9,DeltaMax*xPhonon,'r')
hold on
plot(x*1e9,V,'k')
hold off
title(['Background potential (black) and delta phonon (red) with a coupling constant of ',num2str(g1),' eV'])
xlabel('x (nm)')
ylabel('Potential energy (eV)')
axis([x(1)*1e9 x(end)*1e9 min((E(EplotIndex)-(NumPhonons-0.5)*omega),0) 1.2*max((E(EplotIndex)+omega/2),max(V))])

subplot(1,3,3)
hold on
for n = 1:NumPhonons
    plot(x*1e9,real(PsiMagPlot(n,:)) + PsiZero(n))
    plot(x*1e9,PsiZero(n)*ones(1,Numx),'g')
end
hold off
title({['Probability density for inelastically scattered electron with injection energy of ',num2str(E(EplotIndex)),' eV'];['Wavefunctions offset so that the zero level is E_{inj} - n\Omega and normalized to prevent overlap'];Normalizations2})
xlabel('x (nm)')
ylabel('Probability density, |\psi|^2 (m^{-1})')
axis([x(1)*1e9 x(end)*1e9 min((E(EplotIndex)-(NumPhonons-0.5)*omega),0) 1.2*max((E(EplotIndex)+omega/2),max(V))])