%Fig7_12ab.m
%Essential Electron Transport for Device Physics
%Plot RPA dielectric function at finite temperature including phonons
clear all;
clf;
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);
%hbar1=['\fontname{MT Extra}h\fontname{Times}'];

colormap(jet);

hbar = 1.05457159e-34;  % Planck's constant (Js)
nDoping1 = 1e18;        % n-type doping concentration (cm^-3) 1e18
nDoping = nDoping1*1e6; % Doping concentration to m^-3
m2 = 0.07;              % Effective electron mass coefficient
m0 = 9.109382e-31;      % Bare electron mass (kg)
wLO=36.3;               % longitudinal optic phonon energy (meV)
wTO=33.3;               % transverse optic phonon energy (meV)

a0=0.529177e-10;        % Bohr radius (m)

einf=11.1;              % high frequency dielectric constant
m = m2*m0;              % Effective electron mass (kg)
echarge = 1.6021764e-19;                        % Electron charge (C)
EFermi = hbar^2*(3*pi^2*nDoping)^(2/3)/(2*m);   % Unit :J
EFermi1 = (EFermi*1e3)/echarge;                 % Fermi energy in meV
wLO2=(wLO/EFermi1)^2;
wTO2=(wTO/EFermi1)^2;
gamma=0.01;              %Energy broadening (GAMMA / Ef) 0.03 0.01
eye=complex(0,1);        %square root of minus one

kFermi = (3*pi^2*nDoping)^(1/3);    % Fermi wave vector (m^-1)
T = 100;                            % Absolute temperature 300, 10, 4.2 (K)
kB = 1.3806505e-23/echarge;         % Boltzmann constant (eV/K)
beta = 1/(kB*T);                    % Inverse thermal energy (eV^-1)
Energymax=100;                      % Electron energy maximum (meV)
Energystep=Energymax/60;            % Electron energy step (meV)
ymax=Energymax/EFermi1;             % max value of energy loss, omega/Ef 18
%max value of scattered wave vector q/kf set at ymax value
xmax=(sqrt(ymax))*(1+sqrt(2));%approximately 15, 6 for n=1e17, 1e18
mu = CalculateChemicalPotential(T, nDoping, m, EFermi1/1e3);
x1=[gamma/2:gamma/2:xmax];%[0.01:.005:xmax] scattered wave vector step at gamma/2
y=[-ymax:gamma/8:ymax]+eye*gamma;% energy loss step is gamma
f1 = zeros(length(y),length(x1));
%*************************************************************************
% Calculate the response and Im_Xe
%*************************************************************************
factor = 4*pi*8.85e-12;
E_q_factor=hbar^2/(2*m*echarge);
Im_Xe_factor=2*echarge^2*m^2/(hbar^4*beta/echarge);
Energy = y*EFermi1;         %energy loss in meV

    x = 0.2;                    %scattered wave vector q/k_F
    q = x*kFermi;               %wave vector in m^-1
    E_q = E_q_factor*q^2;
    Energy1 = (real(Energy)/1e3-E_q).^2/(4*E_q);
    Energy2 = (real(Energy)/1e3+E_q).^2/(4*E_q);
    Im_Xe = (Im_Xe_factor/(q^3))*log((1+exp(-beta*(Energy1-mu)))./...
        (1+exp(-beta*(Energy2-mu))));
    
    Re_Xe = -imag(hilbert(Im_Xe));  %use MATLAB Hilbert transform function
                                    %may require Signal Processing Toolbox
    
    epsi = ((y/(y-eye*gamma)*(Re_Xe+(eye*Im_Xe)))/factor)...
        +(einf*((y.^2-wLO2)./(y.^2-wTO2)));
   
figure(1)    
plot(real(y)*EFermi1,real(epsi),'b');
hold on;
plot(real(y)*EFermi1,imag(epsi),'r');
%axis([-Energymax Energymax -max(1.1*real(epsi)) max(1.1*real(epsi))]);
axis([-Energymax Energymax -150 150]);
xttl=['Energy loss, $\hbar\omega$ (meV)'];
xlabel(xttl,'Interpreter','latex');
ylabel('Real (blue) and Imaginary (red) dielectric function');
ttl1=['\rmFig7.12, GaAs, \itn\rm_0=',num2str(nDoping1,'%4.1e'),' cm^{-3}, \itk\rm_F=',...
num2str(kFermi*1e-3,'%4.1e'),' m^{-1}, \itq\rm/\itk\rm_F=',num2str(x,'%5.2f'),...
', \itm\rm^*_e=',num2str(m2,'%5.2f'),'\times\itm\rm_0, \gamma=',...
num2str(gamma,'%5.2f'),'\times\itE\rm_F, \itT\rm=',num2str(T,'%5.0f'),' K'];
title(ttl1);
hold off;
    
figure(2)    
plot(real(y)*EFermi1,imag(epsi./(abs(epsi).^2)),'b');
hold on;
plot(real(y)*EFermi1,100*imag(epsi./(abs(epsi).^2)),'r');
%axis([-Energymax Energymax -max(1.1*real(epsi)) max(1.1*real(epsi))]);
axis([-Energymax Energymax -25 25]);
xttl=['Energy loss, $\hbar\omega$ (meV)'];
xlabel(xttl,'Interpreter','latex');
ylabel('Dielectric loss function, -Im(1/\epsilon)');
title(ttl1);
hold off;