function [Psi] = InelasticWavefunction(WaveAmps, V, xPhonon, Energy, dx, m, omega, g0, g1, NumPhonons)
% [Psi] = InelasticWavefunction(WaveAmps, V, xPhonon, Energy, dx, m, omega, g0, g1, NumPhonons)
% 
% Description:
% This function calculates the wavefunctions for each inelastic channel
% given the transmission coefficients
%
% Outputs:
% Psi - Spatial vector containing the wavefunction values at each location
%       in 'x'. The first row for each channel contains the coefficients
%       for the forward-moving wave, and the second row for each channel
%       contains those for the wave moving towards negative x.  New
%       channels begin in odd rows.
%
% Inputs:
% WaveAmps   - Outgoing wave amplitudes at the right hand of the system for
%              each inelastic channel.
% 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)
eps0 = 8.854188e-12;            % Permittivity of Free Space (F/m)

% Calculated Parameters
Numx = max(size(V));            % Number of x-Axis Samples

%*************************************************
% Calculate Time-Independent Wavefunction 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;
PsiPlot = zeros(NumPhonons,Numx);

% Back Calculate Wavefunction Values
for xIndex = 1:Numx - 1
    RevxIndex = Numx - xIndex;
    
        % Create k arrays
    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);
end
end