%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize constants   %
%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all;
UENERGY_IN_JOULE = 1.602192E-19;
ULENGTH_IN_METER = 1.E-10;
UMASS_IN_KG = 1.66053E-27;
UTIME_IN_SECOND = ULENGTH_IN_METER*sqrt(UMASS_IN_KG/UENERGY_IN_JOULE);
UFREQ_IN_HZ = 1./UTIME_IN_SECOND/2/pi;
UFREQ_IN_THZ = UFREQ_IN_HZ / 1E12;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize Si perfect crystal LDOS   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load dos_si.out;
freq = 1:dos_si(2);
Si = sum(dos_si(4:3+dos_si(2)))

%%%%%%%%%%
% SiTc   %
%%%%%%%%%%
% 1000 supercell k
clf;
plot((freq-0.5)*dos_si(3)*UFREQ_IN_THZ,dos_si(4:3+dos_si(2))/Si/dos_si(3)/UFREQ_IN_THZ);
xlabel ('THz');
ylabel ('1/THz');
hold on;
load dos_SiTc.out;
plot(dos_SiTc(:,1),filtfilt([0.25 0.25 0.25 0.25],[1],dos_SiTc(:,2)),'o');
title ('<111> LDOS of Si-interstitial at TC site (circles) compared to that in perfect crystal (line)');
axis([0 80 0 0.32]);
PRINT;

%%%%%%%%%%%
% SiTsi   %
%%%%%%%%%%%
% 1000 supercell k
clf;
plot((freq-0.5)*dos_si(3)*UFREQ_IN_THZ,dos_si(4:3+dos_si(2))/Si/dos_si(3)/UFREQ_IN_THZ);
xlabel ('THz');
ylabel ('1/THz');
hold on;
load dos_SiTsi.out;
plot(dos_SiTsi(:,1),filtfilt([0.25 0.25 0.25 0.25],[1],dos_SiTsi(:,2)),'o');
title ('<111> LDOS of Si-interstitial at TSi site (circles) compared to that in perfect crystal (line)');
axis([ 0 40 0 0.27]);
PRINT;

%%%%%%%%%%%%%%%
% Si near Vc  %
%%%%%%%%%%%%%%%
clf;
plot((freq-0.5)*dos_si(3)*UFREQ_IN_THZ,dos_si(4:3+dos_si(2))/Si/dos_si(3)/UFREQ_IN_THZ);
xlabel ('THz');
ylabel ('1/THz');
hold on;
load dos_Vc.out;
plot(dos_Vc(:,1),filtfilt([0.25 0.25 0.25 0.25],[1],dos_Vc(:,2)),'o');
title ('<111> LDOS of nearest Si beside C-vacancy (circles) compared to that in perfect crystal (line)');
axis([ 0 35 0 0.25]);
PRINT;

%%%%%%%%%
% Si_c  %
%%%%%%%%%
clf;
plot((freq-0.5)*dos_si(3)*UFREQ_IN_THZ,dos_si(4:3+dos_si(2))/Si/dos_si(3)/UFREQ_IN_THZ);
xlabel ('THz');
ylabel ('1/THz');
hold on;
load dos_Si_c.out;
plot(dos_Si_c(:,1),filtfilt([0.25 0.25 0.25 0.25],[1],dos_Si_c(:,2)),'o');
title ('<111> LDOS of Si-antisite (circles) compared to that in perfect crystal (line)');
axis([ 0 35 0 0.45]);
PRINT;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize C perfect crystal LDOS   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load dos_c.out;
freq = 1:dos_c(2);
C = sum(dos_c(4:3+dos_c(2)))

%%%%%%%%%
% CTc   %
%%%%%%%%%
% 59 supercell k
clf;
plot((freq-0.5)*dos_c(3)*UFREQ_IN_THZ,dos_c(4:3+dos_c(2))/C/dos_c(3)/UFREQ_IN_THZ);
xlabel ('THz');
ylabel ('1/THz');
hold on;
load dos_CTc.out;
plot(dos_CTc(:,1),filtfilt([0.25 0.25 0.25 0.25],[1],dos_CTc(:,2)),'o');
title ('<111> LDOS of C-interstitial at TC site (circles) compared to that in perfect crystal (line)');
axis([0 35 0 0.6]);
PRINT;

%%%%%%%%%
% CTsi  %
%%%%%%%%%
clf;
plot((freq-0.5)*dos_c(3)*UFREQ_IN_THZ,dos_c(4:3+dos_c(2))/C/dos_c(3)/UFREQ_IN_THZ);
xlabel ('THz');
ylabel ('1/THz');
hold on;
axis([0 35 0 0.7]);
load dos_CTsi.out;
plot(dos_CTsi(:,1),filtfilt([0.25 0.25 0.25 0.25],[1],dos_CTsi(:,2)),'o');
title ('<111> LDOS of C-interstitial at TSi site (circles) compared to that in perfect crystal (line)');
PRINT;

%%%%%%%%%%%%%%%
% C near Vsi  %
%%%%%%%%%%%%%%%
clf;
plot((freq-0.5)*dos_c(3)*UFREQ_IN_THZ,dos_c(4:3+dos_c(2))/C/dos_c(3)/UFREQ_IN_THZ);
xlabel ('THz');
ylabel ('1/THz');
hold on;
axis([0 35 0 0.8]);
load dos_VSi.out;
plot(dos_VSi(:,1),filtfilt([0.25 0.25 0.25 0.25],[1],dos_VSi(:,2)),'o');
title ('<111> LDOS of nearest C beside Si-vacancy (circles) compared to that in perfect crystal (line)');
PRINT;

%%%%%%%%%
% C_si  %
%%%%%%%%%
clf;
plot((freq-0.5)*dos_c(3)*UFREQ_IN_THZ,dos_c(4:3+dos_c(2))/C/dos_c(3)/UFREQ_IN_THZ);
xlabel ('THz');
ylabel ('1/THz');
hold on;
axis([0 40 0 0.8]);
load dos_C_si.out;
plot(dos_C_si(:,1),filtfilt([0.25 0.25 0.25 0.25],[1],dos_C_si(:,2)),'o');
title ('<111> LDOS of C-antisite (circles) compared to that in perfect crystal (line)');
PRINT;

