%%%%%%%%%%%%%%%%%
% Plot out LDOS %
%%%%%%%%%%%%%%%%%

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

filehandle = fopen(fn_ldos, 'r');
[a b] = fscanf (filehandle, '%f');
fclose (filehandle);

freq = a(2+3*a(1):1+3*a(1):b); 

fprintf(1, 'There are %d atoms recorded: ', a(1));
for i=1:a(1), fprintf(1, '%d ', a(3*i-1)); end;
fprintf(1, '\n');
ip = input (['Please choose one to plot out (default=' num2str(a(2)) '): '], 's');
if ~isempty(ip), ip = str2num(ip); else ip = a(2); end;

for i=1:a(1), if (ip==a(3*i-1)) break; end; end;
ldos_x = a(3*i+3*a(1):1+3*a(1):b); 
ldos_y = a(3*i+3*a(1)+1:1+3*a(1):b); 
ldos_z = a(3*i+3*a(1)+2:1+3*a(1):b); 

clf;
plot(freq, ldos_x/(freq(2)-freq(1)), 'r');
hold on;
plot(freq, ldos_y/(freq(2)-freq(1)), 'g');
plot(freq, ldos_z/(freq(2)-freq(1)), 'b');
% stext(0.2, 1.4, '\red x');
% stext(0.2, 1.2, '\green y');
% stext(0.2, 1.0, '\blue z');

xlabel('\nu [THz]');
ylabel(['\rho_{' num2str(ip) '\alpha} (\nu) [1/THz]']);
ans = ['Vibrational LDOS ({\red x}{\green y}{\blue z}) of atom ' num2str(ip) ' in Argon perfect crystal at P=0 (LJ6-12, R_{cut}=2.5)'];
stitle (ans);

axis([0 2.1 0 1.5]);

% print -dpsc ldos.ps 
