clf;

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;

load dat.out;
freq = 1:dat(2);

% zero-phase filter
da = filtfilt([0.6, 0.4],[1],dat(4:3+dat(2)));

if (dat(1)<=6000)
 % must be the DAT for disordered configuration. # of atoms < 1000
 % So we plot them in bar plot.
 bar((freq-0.5)*dat(3)*UFREQ_IN_THZ,dat(4:3+dat(2)));
 title ('Vibrational Eigenfrequency Distribution of Perfect \beta-SiC (216-atom cell) ')
else
 % The perfect crystal case.
% plot((freq-0.5)*dat(3)*UFREQ_IN_THZ,dat(4:3+dat(2)));
 plot((freq-0.5)*dat(3)*UFREQ_IN_THZ,da);
 title ('Total Vibrational Density of States of SiC Crystal')
end;
xlabel ('THz');

PRINT;

% Zero-pt temperature of perfect crystal = 512.564699

