function [Transmission, TransCoeffs, WaveAmps] = InelasticTransmissionSingleE(V, xPhonon, Energy, dx, m, omega, g0, g1, NumPhonons)
% [Transmission, TransCoeffs, WaveAmps] = InelasticTransmissionSingleE(V, xPhonon, Energy, dx, m, omega, g0, g1, NumPhonons)
% 
% Description:
% This function calculates the transmission through a system containing a
% single Einstein phonon for an electron having energy given by 'Energy'.
% Only those inelastic channels that are propagating out of the right hand
% side of the system contribute to the total transmission.  The
% wavefunction amplitudes are also output, providing a means for
% back-calculating the wavefunctions.
%
% Outputs:
% Transmission - Total transmission of the system for
%                each enegy value in 'Energy'
% TransCoeffs  - Vector containing the transmission coefficients for each
%                of the inelastic channels.  The first column is the total
%                transmission, with each of the following columns
%                containing the transmission due to no phonons, 1 phonon,
%                2 phonons, etc.
% WaveAmps     - Save the output amplitudes of the wavefunctions in the
%                same order as TransCoeffs less the first row.
%
% Inputs:
% V          - Spatial vector containing the potential profile of the system (V)
% xPhonon    - Spatial vector containing a 1 at the location of the delta
%              and zeros elsewhere
% Energy     - Energy of injected electron in eV
% dx         - Step in x-axis values in m
% m          - Electron mass in kg
% omega      - Energy of the phonon in eV
% g0         - Static delta potential in eV
% g1         - Electron-phonon coupling constant in eV
% NumPhonons - Maximum number of phonons that are allowed to be excited
%              plus one for the elastic channel

%*************************
% Constants and Matrices *
%*************************
% Physical Constants
hbar = 1.05457148e-34;          % Reduced Planck's Constants (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
    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 Wavefunction 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