clear;
format('e',12);

// ************************************
//
//         アルミニウムの物性
//
// ************************************
// アルミニウムの格子定数: a (bohr)
a = 7.6517;
// アルミニウムの価電子数
n = 3;
// アルミニウムの原子体積(bohr3/atom)
// fccは単位格子に原子4個
vat = (a .^ 3) / 4;

// ************************************
//
//              物理定数
//
// ************************************
// アボガドロ数 (/mol)
na = 6.0221413E23;
// 1 (Ry) =  2.179872E-18 (J)
eRy = 2.179872E-18; //(J)
// リュードベリ原子単位系でのボルツマン定数
// Boltzmann constant kB = 1.3806488E-23 (J/K)
kB = 1.3806488E-23 / eRy;  // (Ry/K)

// ************************************
//
//              行列定義
//
// ************************************
// エネルギーベクトルの生成
E = [0:0.001:2];
// 温度ベクトルの生成
tstart = 10;
tstep = 10;
tend = 20000;
T = [tstart:tstep:tend];

// ************************************
//
//             自由電子近似
//
// ************************************
// 自由電子近似でのフェルミエネルギー(Ry)
ef = (3 * %pi * %pi * n / vat) ^ (2/3);
// 自由電子近似での状態密度 (state/Ry/atom)
function Dfree = Dfree(energy);
 Dfree = vat / 2 / %pi / %pi .* (energy ^ (1/2))
endfunction

// ************************************
//
//         ゾンマーフェルト展開
//
// ************************************
// ゾンマーフェルト展開と自由電子近似による化学ポテンシャル
Musomm = ef - %pi * %pi * kB * kB / 12 / ef .* T .* T;
// ゾンマーフェルト展開と自由電子近似による電子比熱
Cesomm = %pi * %pi * kB * kB * n / 2 / ef .* T;

// ************************************
//
//              数値計算
//
// ************************************
// フェルミ分布関数の定義
function fermi = fermi(mu, energy, temp);
fermi = 1 ./ (exp((energy - mu) ./ (kB * temp)) + 1)
endfunction
// 数値計算による化学ポテンシャル
deff('[y] = f1(x,temp)','y = intsplin(E,Dfree(E) .* fermi(x,E,temp)) - 3');
Munum = ones(T);
Uenum = ones(T);
for i = 1:length(T) do
  temp = T(i);
  Munum(i) = fsolve(0,f1);
  Uenum(i) = inttrap(E,Dfree(E) .* fermi(Munum(i),E,T(i)) .* E);	// 内部エネルギー
end
// 数値計算による電子比熱
dUdT = diff(Uenum) / tstep;
Cenum = [0,dUdT];

// ************************************
//
//             グラフ描画
//
// ************************************
// 化学ポテンシャルのグラフ描画
scf(0)
xsetech([0,0,0.95,0.95]);
// Sommerfeld展開による化学ポテンシャル
plot(T,Musomm - ef,'-b');
// 数値積分による化学ポテンシャル
plot(T,Munum - ef,'--r');
legend(['Sommerfeld expansion';'Numerical calculation'],1);
xlabel('Temperature (K)');
ylabel('Chemical potential (Ry)');
// インセットで描画
xsetech([0.15,0.38,0.4,0.4]);
// 状態密度
plot(E-ef,Dfree(E));
// 最低温度()
plot(E-ef,Dfree(E).*fermi(Munum(1),E,T(1)),'-g');
// 最高温度(20,000K)での電子の分布
plot(E-ef,Dfree(E).*fermi(Munum(length(T)),E,T(length(T))),'-r');
legend(['DOS';'T = 10 K';'T = 20,000 K'],2);
xlabel('Energy (Ry)');
ylabel('Densitiy of states (state/Ry/atom)');
// 電子比熱のグラフ描画
scf(1)
xsetech([0,0,0.95,0.95]);
// ゾンマーフェルト展開による電子比熱
plot(T,eRy * na * Cesomm,'-b');
// 数値計算による電子比熱
plot(T,eRy * na * Cenum,'--r');
legend(['Sommerfeld expansion';'Numerical calculation'],2);
xlabel('Temperature (K)');
ylabel('Electronic specific heat (J/mol/K)');