%%%%%%%%
% cv.m %
%%%%%%%%

loadf;
BOLZ_IN_J__K = 1.380658e-23;
HBAR_IN_J_S__R = 6.6260755e-34 / 2 / pi;

TMAX = 3000;
T = 1 : 10 : TMAX;
CV = 0 * T;
batchk = 100;
batchw = 3 * npa * batchk;

filename = ['/tmp/' fname '.acc'];
fprintf (1, 'Loading from accumulation file ''%s''...\n', filename);
fp = fopen(filename, 'r');

kbatch = 0;
more off;
while(1)
  [a,b] = fscanf(fp, '%e', [1+npa, batchw]);
  if (b ~= (1+npa)*batchw) break; end;
  w = sqrt(a(1,:));
  for l = 1 : batchw,
    x = HBAR_IN_J_S__R * w(l) / BOLZ_IN_J__K ./ T;
    y = exp(x);
    CV = CV + (y-1).^-2 .* y .* x.^2;
  end;
  kbatch = kbatch+1;
  clf;
  plot (T, CV/kbatch/batchw);
  xlabel ('T [K]');
  ylabel ('C_V / 3Nk_B');
  title (['Harmonic specific heat (constant-volume) of rocksalt ZrC']);
  grid on;
  fprintf(1, 'Batch %d x %d k\n', kbatch, batchk);
  drawnow;
end;

axis([0 1000 0 1]);
print -depsc cv.eps;
