%Fig2_6
%Essential Electron Transport for Device Physics
%close all
clear all;
clf;
%plotting parameters + fontsizes
FS      = 12;	%label fontsize 18
FSN     = 12;	%number fontsize 16
LW      = 1;	%linewidth

% Change default axes fonts.
set(0,'DefaultAxesFontName', 'Times');
set(0,'DefaultAxesFontSize', FSN);

% Change default text fonts.
set(0,'DefaultTextFontname', 'Times');
set(0,'DefaultTextFontSize', FSN);

hbar=1.054571596e-34;										
echarge=1.602176462e-19;									
baremass=9.10938188e-31;
meff = 0.07;
me = meff*baremass;

FS = 14;
LW = 1.5;
MS = 8;

numsol = 4;

length = 10e-9; %length of potential well [m] 10, use 20 with step change in potential
N  = 0.3e3+1;   %number of spatial samples
dx = length/(N-1);
x  = linspace(0,length,N);

v = zeros(1,N); %rectangular potential with ininfite barrier energy
%step change in potential
%v(1,1:N/2)=echarge*0.3;%step change in potential 0.3

[E,psi] = solve_sch(length,N,v,me,numsol);%solve for rectangular potential

EeV = E/echarge;
psi = [zeros(numsol,1) psi zeros(numsol,1)];
psi0 = psi(1,:);          %ground-state wave function
dpsi0 = gradient(psi0)/dx;  %numberically evaluate gradient
E0 = E(1,1);            %ground-state energy eigenvalue

vn = v-E0;                      %subtract E0 to create vn potential
%create supersymmetric pair potential vp numerically
vp = -vn+(hbar^2/me)*(dpsi0./psi0).^2;
%create supersymmetric pair potential vp analytically for rectangular well
%vp = -vn+(pi^2*hbar^2/(length^2*me))*((cot(x*pi/length)).^2);

[En,psin] = solve_sch(length,N,vn,me,numsol);%solve for vn potential
[Ep,psip] = solve_sch(length,N,vp,me,numsol);%solve for vp potential

%break symmetry by applying uniform electric field
Nstark = 50;%100+1 total potential drop is Nstark * E_0

Astark = linspace(0,Nstark,Nstark);
E1n = zeros(1,Nstark);          E1p = zeros(1,Nstark);
Enstark = zeros(numsol,Nstark); Epstark = zeros(numsol,Nstark);
psinstark=zeros(numsol,N);      psipstark=zeros(numsol,N);

for ii = 1:Nstark
    vstark = (ii-1)*E0*(x-length/2)/length;

    vnstark = vn+vstark;
    vpstark = vp+vstark;

    [Enstarktemp,psinstark] = solve_sch(length,N,vnstark,me,numsol);
    [Epstarktemp,psipstark] = solve_sch(length,N,vpstark,me,numsol);

    diagEn = diag(En).';
    diagEp = diag(Ep).';
    diagEnstark = diag(Enstarktemp).';
    diagEpstark = diag(Epstarktemp).';
    
    E1n(ii) = diagEn(2);
    E1p(ii) = diagEp(1);
    Enstark(:,ii) = diagEnstark;
    Epstark(:,ii) = diagEpstark;
end
%**************************************************************************
%**************** plot results ********************************************
%**************************************************************************
disp(diag(En/echarge))
disp(diag(Ep/echarge))
X = ['E0 = ',num2str(E0/echarge),' eV'];disp(X)

ttl=['\rmm_{eff} = ',num2str(meff),'\timesm_0, L = ',...
    num2str(length*1e9),' nm, E = 0 V/m'];
ttl2=['\rm Fig2.6 m_{eff} = ',num2str(meff),'\timesm_0, L = ',...
    num2str(length*1e9),' nm, E_0 = ',num2str(E0/echarge),...
    ' eV, \Delta_{max} = ',num2str(Nstark*E0/echarge),' eV'];
% 
% figure(1)
% plot(1e9*x(2:end-1),psin,'linewidth',LW)
% grid on
% title(ttl);
% xlabel('Position, $x$ (nm)','FontSize',14,'Interpreter','latex')
% ylabel('$\psi^{-}$','FontSize',14,'Interpreter','latex')
% 
% figure(2)
% plot(1e9*x(2:end-1),psip,'linewidth',LW)
% grid on
% title(ttl);
% xlabel('Position, $x$ (nm)','FontSize',14,'Interpreter','latex')
% ylabel('$\psi^{+}$','FontSize',14,'Interpreter','latex')
% 
% figure(3)
% plot(1e9*x,vn/echarge,'linewidth',LW)
% grid on
% hold on
% plot(1e9*x,(vn+vstark)/echarge,'r','linewidth',LW)
% title(ttl2);
% xlabel('Position, $x$ (nm)','FontSize',14,'Interpreter','latex')
% ylabel('Potential energy, $V^{-} (eV)$','FontSize',14,'Interpreter','latex')
% ylim([-Nstark*E0/echarge 2*Nstark*E0/echarge])
% hold off
% 
% figure(4)
% plot(1e9*x,vp/echarge,'linewidth',LW)
% grid on
% hold on
% plot(1e9*x,(vp+vstark)/echarge,'r','linewidth',LW)
% title(ttl2);
% xlabel('Position, $x$ (nm)','FontSize',14,'Interpreter','latex')
% ylabel('Potential energy, $V^{(+)} (eV)$','FontSize',14,'Interpreter','latex')
% ylim([-Nstark*E0/echarge 2*Nstark*E0/echarge])
% hold off
% 
% figure(5)
% plot(1e9*x(2:end-1),psinstark,':','linewidth',LW)
% grid on
% title(ttl2);
% xlabel('Position, $x$ (nm)','FontSize',14,'Interpreter','latex')
% ylabel('$\psi^{-}$','FontSize',14,'Interpreter','latex')
% 
% figure(6)
% plot(1e9*x(2:end-1),psipstark,':','linewidth',LW)
% grid on
% title(ttl2);
% xlabel('Position, $x$ (nm)','FontSize',14,'Interpreter','latex')
% ylabel('$\psi^{+}$','FontSize',14,'Interpreter','latex')

xplot=Astark*E0/echarge;%setup x-axis as potential energy perturbation, W
figure(1)
plot(xplot,Enstark/echarge,xplot,Epstark(1:3,:)/echarge,'--','linewidth',LW);
grid on
title(ttl2);
xlabel('$W$ ($e|\textbf{E}|L$) (eV)','FontSize',14,'Interpreter','latex')
ylabel('Eigenenergy, $E$ (eV)','FontSize',14,'Interpreter','latex')
xlim([0 max(xplot)])
% legend('- stark','+ stark');

% figure(8)
% plot(xplot,Enstark(2,:)/echarge,xplot,Epstark(1,:)/echarge,'linewidth',LW);
% grid on
% title(ttl2);
% xlabel('$W$ ($e|\textbf{E}|L$) (eV)','FontSize',14,'Interpreter','latex')
% ylabel('Eigenenergy, $E$ (eV)','FontSize',14,'Interpreter','latex')
% legend('- stark','+ stark');
% xlim([0 max(xplot)])
% 
% figure(9)
% plot(xplot,Enstark(3,:)/echarge,xplot,Epstark(2,:)/echarge,'linewidth',LW);
% grid on
% title(ttl2);
% xlabel('$W$ ($e|\textbf{E}|L$) (eV)','FontSize',14,'Interpreter','latex')
% ylabel('Eigenenergy, $E$ (eV)','FontSize',14,'Interpreter','latex')
% legend('- stark','+ stark');
% xlim([0 max(xplot)])
% 
% figure(7)
% plot(xplot,4*del2(Enstark/echarge),xplot,4*del2(Epstark/echarge),'--','linewidth',LW);
% grid on
% title(ttl2);
% xlabel('$W$ ($e|\textbf{E}|L$) (eV)','FontSize',14,'Interpreter','latex')
% ylabel('4*del2, $E^{(-)}$ solid line, $E^{(+)}$ dashed line','FontSize',14,'Interpreter','latex')
% xlim([0 max(xplot)])
% 


%%

function [E,psi] = solve_sch(length,N,v,me,numsol)
    hbar=1.054571596e-34;										
    echarge=1.602176462e-19;   
    dx = length/(N-1);
    H = zeros(N-2,N-2);
    const = hbar^2/(me*dx^2);

    for ii = 1:N-2
        H(ii,ii) = const+v(ii+1);
        if(ii<N-2)
            H(ii+1,ii) = -const/2;
            H(ii,ii+1) = -const/2;
        end
    end

    %select out (numsol=4) lowest energy states
    [psi,E] = eig(H);    
    psi=psi(:,1:numsol);
    E=E(1:numsol,1:numsol);
    
%     E = E/echarge;
    psi = psi.';
end