%%%%%%%%%
% tot.m %
%%%%%%%%%

loadf;
more off;
batchk = 1000;
batchw = 3 * npa * batchk;
numin = 0;
numax = 20;
numesh = 200;
nudel = (numax-numin) / numesh;
nu = numin + nudel/2 : nudel : numax - nudel/2;
DOS = zeros(1+npa, numesh);
filename = ['/tmp/' fname '.acc'];
fprintf (1, 'Loading from accumulation file ''%s''...\n', filename);
fp = fopen(filename, 'r');
color = {'', 'r:', 'g:'};
for m = 4 : 1+npa,
  color{m} = '';
end;
kbatch = 0;
while(1)
  [a,b] = fscanf(fp, '%e', [1+npa, batchw]);
  if (b ~= (1+npa)*batchw) break; end;
  n = ceil( sqrt(a(1,:)) / 2 / pi / 1e12 / nudel );
  for l = 1 : batchw,
    if ((n(l) >= 1) & (n(l) <= numesh))
      DOS(1,n(l)) = DOS(1,n(l)) + 1;
      for m = 2 : 1+npa,
        DOS(m,n(l)) = DOS(m,n(l)) + a(m,l);
      end;
    end;
  end;
  kbatch = kbatch + 1;
  clf;
  plot (nu, DOS(1,:)/kbatch/batchw/nudel, color{1});
  hold on;
  for m = 2 : 1+npa;
    plot (nu, DOS(m,:)/kbatch/batchw/nudel, color{m});
  end;
  xlabel ('\nu [THz]');
  ylabel ('DOS [THz^{-1}]');
  v=axis; axis([numin numax 0 .3]);
  title ('Phonon DOS of Rocksalt ZrC');
  fprintf(1, 'Batch %d x %d k\n', kbatch, batchk);
  drawnow;
end;

print -depsc tot.eps;
