clear;

// *** 物理定数 ***
// アボガドロ数 (/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)
// 電気素量
chage = 1.60217662E-19;

// *** 状態密度の読み出し ***
X = fscanfMat("Kelvin-Pt-calc.txt");
Edat = X(:,1);
Ddat = 2 * X(:,2);

// *** 計算用 ***
// エネルギー
E = linspace(min(Edat), max(Edat), 10000);
// 温度
tstart = 10; tstep = 10; tend = 2500;
T = [tstart:tstep:tend];
// 状態密度
D = interp1(Edat, Ddat, E, "linear");

// *** フェルミ分布関数 ***
function fermi = fermi(mu, energy, temp)
  fermi = 1 ./ (exp((energy - mu) ./ (kB * temp)) + 1)
endfunction

// *** フェルミ分布関数 ***
n = intsplin(E, (D .* fermi(0, E, 0)));
function y = f1(x, temp)
  y = intsplin(E, D .* fermi(x, E, temp)) - n
endfunction

// *** 化学ポテンシャル ***
Snum = ones(T);
for i = 1:length(T) do
  temp = 1.01 * T(i);
  mu1 = fsolve(0, f1);
  temp = 0.99 * T(i);
  mu2 = fsolve(0, f1);
  // 数値計算による電子比熱
  Snum(i) =  - eRy * (mu1 - mu2) / (0.02 * T(i)) / chage;
end

// 数値計算による電子比熱
plot(T, 1E6 * Snum, "-b");
Y = fscanfMat("Kelvin-Pt-lit.txt");
plot(Y(:,1), Y(:,2), ".b");

// *** グラフの装飾 ***
xlabel("Temperature (K)");
ylabel("Seebeck coefficient (uV/K)");
xgrid(color("gray"));