%%%%%%%%%
% amp.m %
%%%%%%%%%

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

a0 = 4.695;
TMAX = 4500;
T = 1 : 10 : TMAX;
amp2 = zeros(1+npa,1) * T;
totalmass = sum(mass);
batchk = 100;
batchw = 3 * npa * batchk;

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;
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,
    occ = .5 + 1 ./ ( exp(HBAR_IN_J_S__R * w(l) / ...
                          BOLZ_IN_J__K ./ T) - 1 );
    amp2 = amp2 + [sum(a(2:1+npa,l)./mass)/npa; ...
                   a(2:1+npa,l)./mass] * occ * HBAR_IN_J_S__R ./ w(l);
  end;
  kbatch = kbatch+1;
  clf;
  plot(297.13, 0.105, 'g^'); hold on;
  a = [25 832 1100 1430 1700; 0.087 0.161 0.185 0.206 0.225];
  plot(a(1,:)+297.13, a(2,:), 'o');
  legend('carbon (exp''t)', 'avg. (exp''t)', 4);
  for m = 1 : 1+npa;
    plot (T, sqrt(amp2(m,:)/kbatch/batchk)/1e-10, color{m});
  end;
  xlabel ('T [K]');
  sylabel ('<|\Delta{\b r}|^2>^{1/2} [A]');
  title (['r.m.s. vibrational amplitudes of Zr and C atoms in Rocksalt' ...
          ' ZrC']);
%   t = 0.083 * sqrt(3) * a0 / 2; %
  t = 0.12 * a0 / 2;
  axis ([0 TMAX 0 0.37]);
  stext(1800, 0.19,'\red<|\Delta{\b r}_{Zr}|^2>^{1/2}');
  stext(900, 0.21,'\green<|\Delta{\b r}_C|^2>^{1/2}');
  v = axis; plot ([v(1) v(2)], [t t], 'm');
  stext(1000, 0.35, ['\magenta Lindemann/Gilvarry law: '...
                     '<|\Delta{\b r}|^2>^{1/2} = 0.12(a_0/2)']);
  grid on;
  fprintf(1, 'Batch %d x %d k\n', kbatch, batchk);
  drawnow;
end;

print -depsc amp.eps;
