function InelasticWavefunctions
clear all;
clc
warning off

%Essential Electron Transport for Device Physics
%****************************************************
%                                                   *
%                   Main program                    *
%                                                   *
%****************************************************
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 = .25;      % Phonon energy (eV) 0.25
g0 = 0*1e-9;      % Static delta barrier energy (eVm)
g1 = .01*1e-9;    % Electron-phonon coupling constant (eVm)
NumPhonons = 3;   % Number of phonons included in calculation 1, 3

% Energy parameters
Energy = .058;    % Energy Value for Plotting of Wavefunction 0.1

% Bias parameters
VBias = 0;        % Voltage bias (V)
VShift = 0;       % Shifts entire potential profile (V)

% Spatial Parameters
Numx = 1200;      % Number of points for potential and wave function plots

% Material parameters
meff = 0.07;      % Effective electron mass

% Barrier Parameters
  % Potential Energy of Each Segment (eV)
Potential = [0 -.5 -.5 -.5 -.3];
  % Length of Each Segment (m)
Length =    [20 20.0 5.0 5.0 20]*1e-9; %10 5/5.5 5 5 10
  % 1 Indicates which barrier phonon is located at, 
  % 2 in last entry for phonon at last barrier
Phonon =    [0 0 1 0 0];

%*************************
% Constants and matrices *
%*************************
% Physical constants
hbar = 1.05457148e-34;          % Reduced Planck's constant (Js)
echarge = 1.60217646e-19;       % Electron charge (C)
m0 = 9.109382e-31;              % Electron mass (kg)
eps0 = 8.854188e-12;            % Permittivity of free space (F/m)

% Modify some parameters
m = meff*m0;					% Effective electron mass (kg)
NumPhonons = NumPhonons + 1;    % Include zero phonon

%******************
% Build potential *
%******************
% Make new barrier length array
NumBarrier = length(Length);
Length = [0 Length];
Steps = Length*triu(ones(NumBarrier+1));

% Create x-axis
xMax = Steps(end);
x = linspace(0, xMax, Numx);
dx = mean(diff(x));

% Get potential profile
V = PotentialProfile(x, Steps, Potential, VShift, VBias);

% If any parameter indicates no phonons, clear all phonon related parameters
if sum(Phonon) == 0 || g1 == 0 || omega == 0 || NumPhonons == 1
    g1 = 0;
    omega = 0;
    NumPhonons = 1;
end

% Find phonon locations
xPhonon = FindPhonon(Steps, Phonon, x);

%**************************
% Calculate wave functions *
%**************************    
% Get wave function coefficients at right edge of system
T = InelasticTransmission(V, xPhonon, Energy, dx, m, omega, g0, g1, NumPhonons);

% Calculate wave function
psi = InelasticWavefunction(T, V, xPhonon, Energy, dx, m, omega, g0, g1, NumPhonons);

PsiPlot = zeros(NumPhonons, Numx);
for PhononIndex = 1:NumPhonons
    PsiPlot(PhononIndex, :) = psi(2*PhononIndex-1,:) + psi(2*PhononIndex,:);
end

Psi2Plot = conj(PsiPlot).*PsiPlot;

% Normalize wave functions
for PhononIndex = 1:NumPhonons
    PsiPlot(PhononIndex,:) = PsiPlot(PhononIndex,:)*omega/max(max(abs(real(PsiPlot(PhononIndex,:)))),max(abs(imag(PsiPlot(PhononIndex,:)))));
    Psi2Plot(PhononIndex,:) = Psi2Plot(PhononIndex,:)/max(Psi2Plot(PhononIndex,:));
end
    
%***************
% Plot results *
%***************
% Renormalize x-axis for plotting
xMax = xMax*1e9;
x = x*1e9;

%*****************************************************************
% Plot potential with overlayed wave functions
ttl1=['\rmFig3.18, \itE\rm=',num2str(Energy),' eV, \itm\rm^*_e=',num2str(meff),...
        '\times\itm\rm_0, h\omega=',num2str(omega),...
        ' eV, \itg\rm_1=',num2str(g1*1e9),' eV nm']; 

figure(1)
% Plot deltas
DeltaLocs = x(logical(xPhonon));
for DeltaIndex = 1:length(DeltaLocs)
line([DeltaLocs(DeltaIndex) DeltaLocs(DeltaIndex)], [-5 5], 'Color', 'r')
end
hold on

% Plot potential
plot(x, V, 'k','Linewidth',LW)

% Plot wave functions
for PhononIndex = 1:NumPhonons
plot(x,(Energy - (PhononIndex - 1)*omega)*ones(1, Numx),'g')
plot(x,real(PsiPlot(PhononIndex,:))*omega + Energy - (PhononIndex - 1)*omega,'b','Linewidth',LW)
plot(x,imag(PsiPlot(PhononIndex,:))*omega + Energy - (PhononIndex - 1)*omega,'r','Linewidth',LW)
end
hold off

title(ttl1)
xlabel('Position, {\itx} (nm)')
ylabel('Re(blue), Im(red), \psi(\itx\rm) and {\it V} (eV)')
xlim([0 xMax])
ylim([1.2*min(0,min(min(V),Energy - (NumPhonons - .5)*omega)) 1.2*max(Energy + omega/2)])

%*****************************************************************
% Plot potential with overlayed probability densities
figure(2)
% Plot deltas
DeltaLocs = x(logical(xPhonon));
for DeltaIndex = 1:length(DeltaLocs)
    line([DeltaLocs(DeltaIndex) DeltaLocs(DeltaIndex)], [-5 5], 'Color', 'r')
end
hold on

% Plot potential
plot(x, V, 'k','Linewidth',LW)

% Plot wave function intensities
for PhononIndex = 1:NumPhonons
    plot(x, (Energy - (PhononIndex - 1)*omega)*ones(1, Numx),'g')
    plot(x, Psi2Plot(PhononIndex,:)*omega*.75 + Energy - (PhononIndex - 1)*omega,'k','Linewidth',LW)
end
hold off

title(ttl1)
xlabel('Position, {\it x} (nm)')
ylabel('|\psi(\itx\rm)|^2 and {\it V} (eV)')
xlim([0 xMax])
ylim([1.2*min(0,min(min(V),Energy - (NumPhonons - 1)*omega)) 1.2*(Energy + omega)])
end

%****************************************************
%                                                   *
%                   Subroutines                     *
%                                                   *
%****************************************************
function [Potential] = PotentialProfile(x, Steps, V, VShift, VBias)
%*************************
% Constants and Matrices *
%*************************
% Calculated Parameters
Numx = max(size(x));                % Number of Samples in x-Axis
NumBarrier = max(size(Steps)) - 1;  % Number of Barriers

d2 = Steps(end);        %total width of potential profile (m)

%***************************
% Create Potential Profile *
%***************************
Potential = zeros(1, Numx);
for xIndex = 1:Numx
    if x(xIndex) < 0
        Potential(xIndex) = 0.0;
    end
    if x(xIndex) >= 0.0
        if x(xIndex) < d2
            for BarrIndex = 1:NumBarrier 
                if x(xIndex) >= Steps(BarrIndex)
                    if x(xIndex) < Steps(BarrIndex + 1)
                        Potential(xIndex) = V(BarrIndex); 
                    end 
                end 
            end
            % Lower potential dut to bias
            Potential(xIndex) = Potential(xIndex) - VBias*x(xIndex)/d2;
        end 
    end
    if x(xIndex) >= d2
        if VBias == 0
            Potential(xIndex) = V(end);
        else
            Potential(xIndex) = -VBias;
        end
    end
end

% Offset Potential Profile
Potential = Potential + VShift;
end

function [xPhonon] = FindPhonon(Length, Phonon, x)
if sum(Phonon) == 0
    % If no phonon, return empty array
    xPhonon = zeros(1,length(x));
    
else
    if Phonon(end) == 2
    % Put last phonon on end of Phonon
        Phonon(end) = 0;
        Phonon = [Phonon 1];
    end
    
    % Must find phonon
    xPhonon = zeros(1,length(x));

    % Find phonon locations
    PhononLocations = Length(logical(Phonon));
    
    for PhononIndex = 1:length(PhononLocations)
        [z,PhononLocIndex] = min(abs(x - PhononLocations(PhononIndex)));
        
% Adjust phonon so that potential change occur between PhononLocIndex and PhononLocIndex  1
        if x(PhononLocIndex) > PhononLocations(PhononIndex)
            PhononLocIndex = PhononLocIndex - 1;
        end
    
        xPhonon(PhononLocIndex) = 1;
    end
end
end

function [WaveAmps] = InelasticTransmission(V, xPhonon, Energy, dx, m, omega, g0, g1, NumPhonons)
%*************************
% Constants and Matrices *
%*************************
% Physical Constants
hbar = 1.05457148e-34;          % Reduced Planck's constant (Js)
echarge = 1.60217646e-19;       % Electron charge (C)

% Calculated Parameters
Numx = max(size(V));            % Number of x-axis samples

% Initialize Transmission and TransCoeffs
Transmission = 0;
TransCoeffs = zeros(1,NumPhonons + 1);

%******************************************************
% Calculate transmission and reflection due to system *
%******************************************************
% If any parameter indicates no phonons, clear all phonon related parameters
if sum(xPhonon) == 0 || g1 == 0 || omega == 0 || NumPhonons == 1
    g1 = 0;
    omega = 0;
    NumPhonons = 1;
end

% Loop Over Injection Energies
P = eye(2*NumPhonons);  % Reset P
for xIndex = 1:Numx-1  % Loop through x-axis

    % Create k arrays
    k1 = zeros(1, NumPhonons + 2);
    k2 = zeros(1, NumPhonons + 2);
    for PhononIndex = 1:NumPhonons+2
        k1(PhononIndex) = sqrt(2*m*echarge*(Energy-(PhononIndex-2)*omega - V(xIndex)))/hbar;
        k2(PhononIndex) = sqrt(2*m*echarge*(Energy-(PhononIndex-2)*omega - V(xIndex+1)))/hbar;
    end
    if xPhonon(xIndex) == 1
        % Create Inelastic Step Matrix
        Pstep = zeros(2*NumPhonons,2*NumPhonons+4);

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

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

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

        Pstep = Pstep(:,3:2*NumPhonons+2);
    else
        % Create Elastic Step Matrix
        Pstep = zeros(2*NumPhonons);
        for PhononIndex = 1:NumPhonons
            Pstep(2*PhononIndex-1,2*PhononIndex-1) = 0.5*(1 + k2(PhononIndex+1)/k1(PhononIndex+1));
            Pstep(2*PhononIndex,2*PhononIndex) = Pstep(2*PhononIndex-1,2*PhononIndex-1);
            Pstep(2*PhononIndex,2*PhononIndex-1) = 0.5*(1 - k2(PhononIndex+1)/k1(PhononIndex+1));
            Pstep(2*PhononIndex-1,2*PhononIndex) = Pstep(2*PhononIndex,2*PhononIndex-1);
        end
    end            

    % Create propagation matrix
    Pprop = zeros(2*NumPhonons);
    for PhononIndex = 1:NumPhonons
        Pprop(2*PhononIndex-1,2*PhononIndex-1) = exp(-i*k2(PhononIndex+1)*dx);
        Pprop(2*PhononIndex,2*PhononIndex) = exp(i*k2(PhononIndex+1)*dx);
    end

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

% Isolate P matrix elements needed to solve for c's
P = P(1:2:2*NumPhonons,1:2:2*NumPhonons);

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

% Solve system
T = P\A;

% Save wave function amplitudes
WaveAmps = T;

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

% Calculate k values for transmission
kin = sqrt(2*m*echarge*(Energy - V(1)))/hbar;
kout = zeros(1,nMax);
for PhononIndex = 1:nMax
    kout(PhononIndex) = sqrt(2*m*echarge*(Energy - (PhononIndex-1)*omega - V(end)))/hbar;
end

% Add Transmission coefficients for propogating modes
for PhononIndex = 1:nMax
   Transmission = real(Transmission + kout(PhononIndex)*T(PhononIndex)*conj(T(PhononIndex))/kin);
   TransCoeffs(1,PhononIndex+1) = real(kout(PhononIndex)*T(PhononIndex)*conj(T(PhononIndex))/kin);
end
TransCoeffs(1,1) = Transmission;
end

function [Psi] = InelasticWavefunction(WaveAmps, V, xPhonon, Energy, dx, m, omega, g0, g1, NumPhonons)
%*************************
% Constants and Matrices *
%*************************
% Physical Constants
hbar = 1.05457148e-34;          % Planck constant (Js)
echarge = 1.60217646e-19;       % Electron charge (C)

% Calculated parameters
Numx = max(size(V));            % Number of x-axis samples

% Find first x index where potential changes
xIndex = 1;
if ~isempty(find(V, 1))
    while V(xIndex) == V(1)
        xIndex = xIndex + 1;
    end
    xMin = xIndex - 1;
elseif sum(xPhonon) ~= 0
    xMin = find(xPhonon == 1, 1, 'first')
else
    xMin = Numx + 1;
end

%*************************************************
% Calculate time-independent wave function values *
%*************************************************
% If any parameter indicates no phonons, clear all phonon related parameters
if sum(xPhonon) == 0 || g1 == 0 || omega == 0 || NumPhonons == 1
    g1 = 0;
    omega = 0;
    NumPhonons = 1;
end

% Put d values into T matrix
T = zeros(2*NumPhonons,1);
T(1:2:2*NumPhonons,1) = WaveAmps;

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

% Back calculate wavefunction values
for xIndex = 1:Numx - 1
    RevxIndex = Numx - xIndex;
    
    % Create k arrays
    k1 = zeros(1, NumPhonons + 2);
    k2 = zeros(1, NumPhonons + 2);
    for PhononIndex = 1:NumPhonons+2
        k1(PhononIndex) = sqrt(2*m*echarge*(Energy-(PhononIndex-2)*omega - V(RevxIndex)))/hbar;
        k2(PhononIndex) = sqrt(2*m*echarge*(Energy-(PhononIndex-2)*omega - V(RevxIndex+1)))/hbar;
    end
        
    if xPhonon(RevxIndex) == 1
        % Create Inelastic Step Matrix
        Pstep = zeros(2*NumPhonons,2*NumPhonons+4);

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

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

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

        Pstep = Pstep(:,3:2*NumPhonons+2);
    else
        % Create Elastic Step Matrix
        Pstep = zeros(2*NumPhonons);
        for PhononIndex = 1:NumPhonons
            Pstep(2*PhononIndex-1,2*PhononIndex-1) = 0.5*(1 + k2(PhononIndex+1)/k1(PhononIndex+1));
            Pstep(2*PhononIndex,2*PhononIndex) = Pstep(2*PhononIndex-1,2*PhononIndex-1);
            Pstep(2*PhononIndex,2*PhononIndex-1) = 0.5*(1 - k2(PhononIndex+1)/k1(PhononIndex+1));
            Pstep(2*PhononIndex-1,2*PhononIndex) = Pstep(2*PhononIndex,2*PhononIndex-1);
        end
    end

    % Create propagation matrix
    Pprop = zeros(2*NumPhonons);
    for PhononIndex = 1:NumPhonons
        Pprop(2*PhononIndex-1,2*PhononIndex-1) = exp(-i*k2(PhononIndex+1)*dx);
        Pprop(2*PhononIndex,2*PhononIndex) = exp(i*k2(PhononIndex+1)*dx);
    end

    Psi(:,RevxIndex) = Pstep*Pprop*Psi(:,RevxIndex+1);
    
    if RevxIndex == xMin
        for PhononIndex = 1:NumPhonons
            if ~isreal(k1(PhononIndex + 1))
                Psi(2*PhononIndex - 1, RevxIndex) = 0;
            end
        end
    end
end
end