// アルミニウムの電子比熱係数(J/mol/K^2)
egamma = 1.35e-3;
// アルミニウムのデバイ温度
dt = 428;
// 気体定数 (J/K/mol)
r = 8.314

// 格子比熱 (Debye model)
function Cl = Cl(T)
Cl = 9 * r * ((T ./ dt) .^ 3) .* integrate('(x .^ 4) .* exp(x) ./ ((exp(x) - 1) .^ 2)','x',0,dt ./ T);
endfunction
// 電子比熱
function Ce = Ce(T)
Ce = egamma .* T;
endfunction

// 高温までのプロット
// 温度ベクトル
T = [1:1:500];
// 格子比熱のプロット
//     絶対零度の計算は出来ないので後から補う
plot([0,T],[0,Cl(T)],'-b');
// 電子比熱のプロット
plot([0,T],[0,Ce(T)],'--r');
legend(['Lattice specific heat';'Electronic specific heat'],2);
xlabel("Temperature (K)");
ylabel("Specific heat (J/K/mol)");

// 低温部分のプロット
xsetech([0.4,0.32,0.5,0.5]);
// 温度ベクトル
T = [1:0.1:10];
// 格子比熱のプロット
//     絶対零度の計算は出来ないので後から補う
plot([0,T],[0,Cl(T)],'-b');
// 電子比熱のプロット
plot([0,T],[0,Ce(T)],'--r');
xlabel("Temperature (K)");
ylabel("Specific heat (J/K/mol)");