%%%%%%%%%%%%%%%
% gruneisen.m %
%%%%%%%%%%%%%%%

BOLZ_IN_J__K = 1.380658e-23;
HBAR_IN_J_S__R = 6.6260755e-34 / 2 / pi;
AVO = 6.02214199e+23;
CAL_IN_J = 4.1868;

TMAX = 3000;
T = 1 : 10 : TMAX;
batchk = 5000;
ks = rand(batchk,3);

fname = 'fconst';
loadf;
batchw = 3 * npa * batchk;
w  = zeros(batchw,1);
wg = zeros(batchw,1);
for l = 1 : batchk,
  w(3*npa*(l-1)+1:3*npa*l) = sqrt(sort(dk(ks(l,:))));
end;
fname = 'fconst.g';
loadf;
for l = 1 : batchk,
  wg(3*npa*(l-1)+1:3*npa*l) = sqrt(sort(dk(ks(l,:))));
end;

dlnV = 3 * log(1.001);
g = 0 * T;
CV = 0 * T;
for l = 1 : batchw,
  x = HBAR_IN_J_S__R * w(l) / BOLZ_IN_J__K ./ T;
  y = exp(x);
  cv = (y-1).^-2 .* y .* x.^2;
  CV = CV + cv;
  g = g - cv * log(wg(l)/w(l)) / dlnV;
end;
g = g ./ CV;
CV = CV / batchw;

clf;
plot (T, g);
xlabel ('T [K]');
ylabel ('Gruneisen Parameter');
title (['Gruneisen parameter of rocksalt ZrC']);
v = axis; axis([0 1000 v(3) v(4)]);
grid on;
print -depsc gruneisen.eps;

input('Press a key ... ');
CV = npa / det(H) * CV * 3 * BOLZ_IN_J__K;
alpha = CV .* g / 3 / B;
clf;
plot (T, alpha);
xlabel ('T [K]');
ylabel ('\alpha');
title (['Thermal expansion coefficient of rocksalt ZrC']);
axis([0 1000 0 7e-6]);
print -depsc gruneisen_expansion.eps;

input('Press a key ... ');
ALPHA = 3 * alpha;
CP = CV + ALPHA.^2 * B .* T;
clf;
plot (T, CV * det(H) / npa * AVO, '--'); hold on;
plot (T, CP * det(H) / npa * AVO );
xlabel ('T [K]');
ylabel ('[J mol^{-1}K^{-1}]');
title (['Specific heats of rocksalt ZrC']);
a = [
    298.15  9.016
    300     9.066
    400     10.76
    500     11.45
    600     11.77
    700     11.93
    800     12.01
    900     12.07
    1000    12.11
    1100    12.16
    1200    12.22
    1300    12.30
    1400    12.39
    1500    12.51
    1600    12.65
    1700    12.81
    1800    13.00
    1900    13.22
    2000    13.46
    2100    13.73
    2200    14.02
    2300    14.34
    2400    14.69
    2500    15.07
    2600    15.47
    2700    15.90
    2800    16.36
    2900    16.84
    3000    17.35 ];
plot(a(:,1), a(:,2)*CAL_IN_J/2, 'o');
legend('Phonon C_V', 'Phonon C_P', 'ZrC_{0.96} exp''t', 4);
axis([0 3000 0 35]);
print -depsc heat_capacities.eps;
