%%%%%%%%%%%%%%%%
% Plot out DOS %
%%%%%%%%%%%%%%%%

if exist('fn_dos'), 
 a = input(['fn_dos = ' fn_dos '? '], 's');
 if ~isempty(a), fn_dos = a; end;
else
 fn_dos = input('fn_dos = ', 's');
end;

filehandle = fopen(fn_dos, 'r');
[a b] = fscanf (filehandle, '%f');
freq = a(1:2:b); dos = a(2:2:b); 
fclose (filehandle);

clf;
plot(freq, dos/(freq(2)-freq(1)));

xlabel('\nu [THz]');
ylabel('\rho(\nu) [1/THz]');
title ('Vibrational DOS of Argon perfect crystal at P=0 (LJ6-12, R_{cut}=2.5)');

axis([0 2.1 0 1.5]);

% print -dpsc dos.ps 
