%Fig7_13ab.m
%Essential Electron Transport for Device Physics
%Inelastic electron scattering rate at finite temperature
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=200;                      % Electron energy maximum (meV) 300
Energystep=Energymax/160;            % Electron energy step (meV)/60
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:ymax]+eye*gamma;% energy loss step is gamma
f1 = zeros(length(y),length(x1));
%*************************************************************************
% Calculate the response and Im_Xe map
%*************************************************************************
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
    for ii=1:length(x1)
    x = x1(ii);
    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
    
    epsi = (y/(y-eye*gamma)*(Re_Xe+(eye*Im_Xe)))/factor+(einf*((y.^2-wLO2)./(y.^2-wTO2)));
    f1(:,ii)=-imag(1./epsi)./x; %loss function in the integrand weighted by 1/x
end
%*************************************************************************
% Bose distribution
%*************************************************************************
Bose_mask = 1.0+(1./(exp(beta*real(y)*EFermi1/1000)-1.0));
Bose_indicator = Bose_mask'*ones(1,length(x1));

Emax_tmp = 0.1:Energystep:Energymax;
spec_e = zeros(length(Emax_tmp),1);
%*************************************************************************
% Integration for electron scattering rate
%*************************************************************************
for ii = 1:length(Emax_tmp)
    Emax = Emax_tmp(ii);

    E_final = Emax-real(y)*EFermi1;
    f_fermi = 1./(exp(beta*(E_final/1000-mu))+1.0);

    Emax1 = Emax/EFermi1;

    yplot=Emax1-(x1-sqrt(Emax1)).^2;  %parabola of integration
    yplot2=Emax1-(x1+sqrt(Emax1)).^2; %lower boundary of integration
%*************************************************************************
%electron-only scattering
%*************************************************************************
  Fermi_indicator = (1.0-f_fermi')*ones(1,length(x1));
   f2 = f1.*Bose_indicator.*Fermi_indicator;
%*************************************************************************
%quasi-particle scattering
%*************************************************************************
%     Fermi_indicator = (f_fermi')*ones(1,length(x1));
%     f2 = f1.*(Bose_indicator-Fermi_indicator);
%*************************************************************************
    f2(isnan(f2))=0;%Array elements that are NaN

    [X1,Y1]=meshgrid(x1,y);
    sumindicator=Y1<=ones(length(y),1)*yplot&Y1>=ones(length(y),1)*yplot2; %find values of sum indicator for each column
    Integral=sum(sum(f2(sumindicator)))/(sqrt(Emax1));  
    spec_e(ii)=Integral;
end
spec = spec_e;

%*************************************************************************
figure(1);
surf(x1,real(y)*EFermi1,log10(abs(f1)));
                az = 0;           %set up the angle of view
                el = 90;          %looking from directly above
                view(az, el);     %set the view
                shading interp;   %do color interpolation
       colorbar;
       axis([0 xmax -Energymax Energymax]);
       hold on;
       plot(x1, yplot*EFermi1,'r','LineWidth',2);
       plot(x1, yplot2*EFermi1,'r','LineWidth',2);
       hold off;
       xlabel('Wave vector, \itq\rm/\itk\rm_F');
       ylabel('Energy loss, $\hbar\omega$ (meV)', 'Interpreter', 'latex');
ttl1=['\rmFig7.13a, GaAs, log_{10}(\itS\rm(\itq\rm, \omega )/\itq\rm), \itn\rm_0=',...
    num2str(nDoping1,'%4.1e'),' cm^{-3}, \itk\rm_F=',...
    num2str(kFermi*1e-3,'%4.1e'),' m^{-1}, \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, \mu=',...
    num2str(mu*1.0e3,'%5.1f'),' meV'];
title(ttl1);
%*************************************************************************
figure(2);
surf(x1,real(y)*EFermi1,log10(abs(f2)));
                az = 0;           %set up the angle of view
                el = 90;          %looking from directly above
                view(az, el);     %set the view
                shading interp;   %do color interpolation
       colorbar;
       axis([0 xmax -Energymax Energymax]);
       hold on;
       plot(x1, yplot*EFermi1,'b','LineWidth',2);
       plot(x1, yplot2*EFermi1,'b','LineWidth',2);
       hold off;
       xlabel('Wave vector, \itq\rm/\itk\rm_F');
       ylabel('Energy loss, $\hbar\omega$ (meV)', 'Interpreter', 'latex');
ttl2=['\rmFig7.13b, GaAs, log_{10}(\itS\rm(\itq\rm, \omega )(1-\itf\rm)\itg\rm/\itq\rm), \itn\rm_0=',...
    num2str(nDoping1,'%4.1e'),' cm^{-3}, \itk\rm_F=',num2str(kFermi*1e-3,...
    '%4.1e'),' m^{-1}, \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, \mu=',...
    num2str(mu*1.0e3,'%5.1f'),' meV'];
title(ttl2);
%*************************************************************************
dx1 = x1(2)-x1(1);
dy1 = y(2)-y(1);
prefactor=(2.0*m2*dx1*dy1*EFermi1*0.001*echarge)/(pi*a0*kFermi*hbar);
%*************************************************************************