%Fig3_20.m
%Essential Electron Transport for Device Physics
clear all
clc
clf
warning off

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 = .036;           % Phonon energy (eV)
g0 = 0*1e-9;            % Static delta barrier energy (eVm)
g1 = 0.1*1e-9;           % Electron-phonon coupling constant (eVm) 0.1 0.5
NumPhonons = 3;         % Number of phonons included in calculation
NumPhononsPlot = 4;     % +1 Number of channels shown in movie when inelastic scattering occurs (must be >= 1)

% Bias Parameters
Bias = 0.0;             % Potential bias applied to device (V)
VShift = 0.0;           % Offset for entire potential profile (V)

% Wavenumber Parameters
kCenter = 2e8;          % Center wavenumber of gaussian (m^-1) (2.1432e8 = 25 meV)
sigmak = 0.5e8;         % Broadening Full-Width-Half-Maximum (eV) (sigmax = 1/sigmak) 0.5
NumGauss = 200;         % Number of Gaussian frequency components

% Spatial Parameters
Numx = 170;             % Number of points in x-axis for potential and wave function plots 151

% Material Parameters
meff = .07;             % Effective mass coefficient
NDoping = 1e18;         % Carrier concentration (cm^-3)
epsr1 = 13.1;           % Relative permittivity of outer regions
epsr2 = 13.1;           % Relative permittivity of barrier regions

% Temporal Parameters
tMin = 0;               % Start time
RunTimeMultiplier = 2;  % Run time of movie in terms of the time it takes for the center of the wave packet to traverse the screen
                        % e.g. Setting equal to 1 will end movie with center of wave packet at right edge of screen
Numt = 350;             % Number of samples of time range 250

% Barrier Parameters
Potential = [0.0 -0.075 -0.075 -0.025];    % Potential in each section
Length =    [70 15 15 70]*1e-9; % Length of each section
Phonon =    [0 0 1 0];          % Phonon locations. A '1' puts a delta barrier at the left edge of the section. Put a '2' in last column to put delta at right edge of domain.

%*************************
% Constants and matrices *
%*************************
% Universal Constants
hbar = 1.0545715e-34;			% Reduced Planck's Constant (J s)
m0 = 9.109382e-31;              % Bare Electron Mass (kg)
echarge = 1.6021764e-19;		% Electron Charge (C)
eps0 = 8.854188e-12;            % Permittivity of Free Space (F/m)

% Modify Parameters
m = meff*m0;					% Effective electron mass (kg)
eps1 = epsr1*eps0;              % Permittivity of Ohmic Contact
eps2 = epsr2*eps0;              % Permittivity of Ohmic Contact
NDoping = NDoping*1e6;          % Modify doping concentration to m^-3
NumPhonons = NumPhonons + 1;    % Add zero phonon

% Calculate Parameters
tMax = RunTimeMultiplier*m*sum(Length)/(hbar*kCenter);
ECenter = hbar^2*kCenter^2/(2*m*echarge)
ESigmak = hbar^2*sigmak^2/(2*m*echarge)

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

% Initialize Matrices
kAmps = zeros(1, NumGauss);
psi = zeros(2*NumPhonons, Numx, NumGauss);
Psi = zeros(Numt, Numx, NumPhonons);

%*********************************************
% Set up arrays and matrices for calculation *
%*********************************************
% Create k array
k = linspace(kCenter - 3*sigmak, kCenter + 3*sigmak, NumGauss);

% Warn user if negative k values present
if k(1) < 0
    disp('Warning: Negative values of k are in Gaussian wave packet and will be ignored')
end

% Create energy array
Energy = hbar^2*k.*k/(2*m*echarge);

% Calculate energy amplitudes
kAmps = sqrt(sigmak/sqrt(pi))*exp(-((kCenter - k).*(kCenter - k))/(2*sigmak^2))/(sigmak*sqrt(2*pi));

% Create Time Array
Time = linspace(tMin, tMax, Numt);

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

% Calculate Poisson parameters
d2 = Steps(NumBarrier+1);  %Total width of potential profile (m)
d1Max = -(eps1*d2)/(2.0*eps2)+sqrt(Bias*eps1/(echarge*NDoping)+(eps1*d2)^2/(2.0*eps2)^2);  % Width of accumulation/depletion region

% Set x-axis limits
xMin = -(d1Max + 2.0e-9);  %lower cutoff in x-direction, depends on NDoping                        
xMax = d2 + d1Max + 2.0e-9;%upper cutoff in x-direction, depends on NDoping

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

% Locate phonon in space
xPhonon = FindPhonons(Steps, Phonon, x);

% Calculate potential profile
V = PotentialProfileWithPoisson(x, Steps, Potential, VShift, Bias, NDoping, eps1, eps2);

% Fix end of potential if no bias
if Bias == 0
    for xIndex = 1:Numx
        if V(Numx + 1 - xIndex) == Potential(end)
            break
        else
            V(Numx + 1 - xIndex) = Potential(end);
        end
    end
end

%**************************
% Calculate Wavefunctions *
%**************************
for kIndex = 1:NumGauss
         % Calculate output wave function amplitudes
        [a, TransCoeffs, WaveAmps] = InelasticTransmissionSingleE(V, xPhonon, Energy(kIndex), dx, m, omega, g0, g1, NumPhonons);
        
        a = InelasticWavefunction(WaveAmps, V, xPhonon, Energy(kIndex), dx, m, omega, g0, g1, NumPhonons);
        % Get total wave functions
        [psi(:,:,kIndex)] = InelasticWavefunction(WaveAmps, V, xPhonon, Energy(kIndex), dx, m, omega, g0, g1, NumPhonons);
        
        % Clear variables for next run
        clear a TransCoeffs WaveAmps
end

%****************************
% Calculate Current Density *
%****************************
for PhononIndex = 1:NumPhonons
    for TimeIndex = 1:Numt
        for kIndex = 1:NumGauss
            % Assemble Total Wavefunctions
            Psi(TimeIndex, :, PhononIndex) = Psi(TimeIndex, :, PhononIndex) + kAmps(kIndex)*(psi(2*PhononIndex - 1, :, kIndex) + psi(2*PhononIndex, :, kIndex))*exp(-i*(Energy(kIndex) - (PhononIndex - 1)*omega)*echarge*Time(TimeIndex)/hbar);
        end
    end
end

PsiMag = Psi.*conj(Psi);

%********************************************
% Calculate Width and Amplitude Versus Time *
%********************************************
sigmax = sqrt(sigmak^-2 + hbar^2*sigmak^2*Time.*Time/(4*m^2));

%***************
% Create Movie *
%***************
for PhononIndex = 1:NumPhonons
    if omega == 0
  Maxes(PhononIndex) = max(max(PsiMag(:, :, PhononIndex)))/0.5;
    else
  Maxes(PhononIndex) = max(max(PsiMag(:, :, PhononIndex)))/(0.75*omega);
    end
end

% Set up movie parameters
avifps = ceil(Numt/10);
if omega == 0
    MovieMax = max(ECenter + 0.5,1.2*max(V));
    MovieMin = min(0, 1.2*min(V));
else
    MovieMax = max(ECenter + omega,1.2*max(V));
    MovieMin = min(1.2*min(Potential), 1.2*(ECenter - (NumPhononsPlot - 1)*omega));
end

% Open movie file

avimov = VideoWriter('Test1.avi');
open(avimov);
%avimov = avifile('Test1.avi','fps',avifps);

% Record movie
%*******************************************************************
ttl=(['\rmFig3.20, ',hbar_txt,'\omega=',...
num2str(omega),' eV, \itm\rm^*_e=',num2str(meff),...
'\times\itm\rm_0, \itg\rm_0=',num2str(g0*1e9),' eV nm, \itg\rm_1=',...
num2str(g1*1e9),' eV nm, \itn\rm_{max}=',num2str(NumPhonons-1)]);
figure(1)
for TimeIndex = 1:Numt
    clf
    plot(x*1e9, V, 'k','LineWidth',LW)
    title(ttl);
    hold on
    if (g0 ~= 0 || g1 ~= 0) && sum(Phonon) ~= 0
        DeltaLoc = find(xPhonon == 1);
        for DeltaIndex = 1:sum(xPhonon)
line([x(DeltaLoc(DeltaIndex))*1e9 x(DeltaLoc(DeltaIndex))*1e9], [V(DeltaLoc(DeltaIndex)) .5], 'Color', 'r','LineWidth',2)
        end
    end
    if omega == 0
plot(x*1e9, ECenter + PsiMag(TimeIndex, :, PhononIndex)/Maxes(PhononIndex))
    else
        for PhononIndex = 1:NumPhononsPlot
plot(x*1e9, ECenter - (PhononIndex - 1)*omega + PsiMag(TimeIndex, :, PhononIndex)/Maxes(PhononIndex),'LineWidth',LW)
        end
    end
    hold off
    xlim([x(1)*1e9 x(end)*1e9])
    ylim([MovieMin MovieMax])
    axis square
xlabel(['Position, \itx\rm (nm)'])
ylabel('|\psi(\itx\rm)|^2, \itV\rm(\itx\rm) (eV)')        
    Mov(TimeIndex) = getframe;
    writeVideo(avimov,Mov(TimeIndex));
    %avimov = addframe(avimov,Mov(TimeIndex));
end

% Close movie file
%avimov = close(avimov);
close(avimov);
%**********************************************************************
%snapshot at Numt/16

figure(2)
    plot(x*1e9, V, 'k','LineWidth',LW)
    title(ttl);
    hold on
    if (g0 ~= 0 || g1 ~= 0) && sum(Phonon) ~= 0
        DeltaLoc = find(xPhonon == 1);
        for DeltaIndex = 1:sum(xPhonon)
line([x(DeltaLoc(DeltaIndex))*1e9 x(DeltaLoc(DeltaIndex))*1e9], [V(DeltaLoc(DeltaIndex)) .5], 'Color', 'r','LineWidth',2)
        end
    end
    if omega == 0
plot(x*1e9, ECenter + PsiMag(ceil(Numt/6), :, PhononIndex)/Maxes(PhononIndex))
    else
        for PhononIndex = 1:NumPhononsPlot
plot(x*1e9, ECenter - (PhononIndex - 1)*omega + PsiMag(ceil(Numt/6), :, PhononIndex)/Maxes(PhononIndex),'LineWidth',LW)
        end
    end
    hold off
    xlim([x(1)*1e9 x(end)*1e9])
    ylim([MovieMin MovieMax])
    axis square
xlabel(['Position, \itx\rm (nm)'])
ylabel('|\psi(\itx\rm)|^2, \itV\rm(\itx\rm) (eV)')        


%********************************************************************
%space time plots
%********************************************************************
xp=linspace(0,sum(Length)*1e9,Numx);
yp=linspace(tMin,tMax*1e12,Numt);
ttl=(['\rmFig3.20, ',hbar_txt,'\omega = ',...
num2str(omega),' eV, \itg\rm_0 = ',num2str(g0*1e9),' eV nm, \itg\rm_1 = ',...
num2str(g1*1e9),' eV nm, \itn\rm = ',num2str(0)]);

figure(3)
pcolor(xp,yp,(PsiMag(:,:,1)));
colormap(jet);shading interp;%winter cool autumn jet gray copper bone
xlabel(['Position, \itx\rm (nm)'])
ylabel('Time, \itt\rm (ps)') 
title(ttl);

ttl=(['\rmFig3.20, ',hbar_txt,'\omega = ',...
num2str(omega),' eV, \itg\rm_0 = ',num2str(g0*1e9),' eV nm, \itg\rm_1 = ',...
num2str(g1*1e9),' eV nm, \itn\rm = ',num2str(1)]);
figure(4)
pcolor(xp,yp,(PsiMag(:,:,2)));
colormap(jet);shading interp;
xlabel(['Position, \itx\rm (nm)'])
ylabel('Time, \itt\rm (ps)') 
title(ttl);

ttl=(['\rmFig.3.20, ',hbar_txt,'\omega = ',...
num2str(omega),' eV, \itg\rm_0 = ',num2str(g0*1e9),' eV nm, \itg\rm_1 = ',...
num2str(g1*1e9),' eV nm, \itn\rm = ',num2str(2)]);
figure(5)
pcolor(xp,yp,(PsiMag(:,:,3)));
colormap(jet);shading interp;
xlabel(['Position, \itx\rm (nm)'])
ylabel('Time, \itt\rm (ps)') 
title(ttl);

ttl=(['\rmFig.3.20, ',hbar_txt,'\omega = ',...
num2str(omega),' eV, \itg\rm_0 = ',num2str(g0*1e9),' eV nm, \itg\rm_1 = ',...
num2str(g1*1e9),' eV nm, \itn\rm = ',num2str(3)]);
figure(6)
pcolor(xp,yp,(PsiMag(:,:,4)));
colormap(jet);shading interp;
xlabel(['Position, \itx\rm (nm)'])
ylabel('Time, \itt\rm (ps)') 
title(ttl);

%located to the far right 
figure(7)
plot(yp,log10(PsiMag(:,ceil(0.9*Numx),1)),'k','LineWidth',LW);
hold on
plot(yp,log10(PsiMag(:,ceil(0.9*Numx),2)),'b','LineWidth',LW);
plot(yp,log10(PsiMag(:,ceil(0.9*Numx),3)),'g','LineWidth',LW);
plot(yp,log10(PsiMag(:,ceil(0.9*Numx),4)),'r','LineWidth',LW);
hold off
xlabel('Time, \itt\rm (ps)') 
ylabel(['log_{10}(|\psi(\itx\rm, \itt\rm)|^2)'])
title(ttl);

figure(8)
plot(yp,(PsiMag(:,ceil(0.9*Numx),1))./max(PsiMag(:,ceil(0.9*Numx),1)),'k','LineWidth',LW);
hold on
plot(yp,(PsiMag(:,ceil(0.9*Numx),2))./max(PsiMag(:,ceil(0.9*Numx),2)),'b','LineWidth',LW);
plot(yp,(PsiMag(:,ceil(0.9*Numx),3))./max(PsiMag(:,ceil(0.9*Numx),3)),'g','LineWidth',LW);
plot(yp,(PsiMag(:,ceil(0.9*Numx),4))./max(PsiMag(:,ceil(0.9*Numx),4)),'r','LineWidth',LW);
hold off
xlabel('Time, \itt\rm (ps)') 
ylabel(['Normalized |\psi(\itx\rm, \itt\rm)|^2'])
title(ttl);

%located at phonon
figure(9)
plot(yp,log10(PsiMag(:,ceil(0.5*Numx),1)),'k','LineWidth',LW);
hold on
plot(yp,log10(PsiMag(:,ceil(0.5*Numx),2)),'b','LineWidth',LW);
plot(yp,log10(PsiMag(:,ceil(0.5*Numx),3)),'g','LineWidth',LW);
plot(yp,log10(PsiMag(:,ceil(0.5*Numx),4)),'r','LineWidth',LW);
hold off
xlabel('Time, \itt\rm (ps)') 
ylabel(['log_{10}(|\psi(\itx\rm, \itt\rm)|^2)'])
title(ttl);

figure(10)
plot(yp,(PsiMag(:,ceil(0.5*Numx),1))./max(PsiMag(:,ceil(0.5*Numx),1)),'k','LineWidth',LW);
hold on
plot(yp,(PsiMag(:,ceil(0.5*Numx),2))./max(PsiMag(:,ceil(0.5*Numx),2)),'b','LineWidth',LW);
plot(yp,(PsiMag(:,ceil(0.5*Numx),3))./max(PsiMag(:,ceil(0.5*Numx),3)),'g','LineWidth',LW);
plot(yp,(PsiMag(:,ceil(0.5*Numx),4))./max(PsiMag(:,ceil(0.5*Numx),4)),'r','LineWidth',LW);
hold off
xlabel('Time, \itt\rm (ps)') 
ylabel(['Normalized |\psi(\itx\rm, \itt\rm)|^2'])
axis([0 max(yp) 0 1.1]);
title(ttl);
