clear;

// *** 外部ファイルの読み込み ***
exec('physics.sci',0);		// 物理定数の読み込み
exec('Moruzzi1988.sci',0);	// 関数の読み込み

// *** 計算設定 ***
// 計算する温度範囲
tstep = 5;
tend = 300 + tstep;
temperature = [0:tstep:tend]';

// *** データの読み出し ***
X = fscanfMat("K-lattice.info");
abcc = X(:,1);	// 格子定数 a (Bohr)
E    = X(:,2);	// エネルギー (Ry)

// *** 物理定数などの読み込み ***
// 平均原子量
M = Mk;
// ボルツマン定数
kBry = kB / ry;
// Wigner-Seitz半径
rws = abcc2rws(abcc);
// プロット用
rplot = linspace(floor(10 * min(rws)) / 10, ceil(10 * max(rws)) / 10);

// *** Rigid Latticeのフィッティング ***
Energy = E;
[fitparam, v] = lsqrsolve([-100;-100;25000;1], fit, size(rws,1));
a_rigid = fitparam(1); b_rigid = fitparam(2); c_rigid = fitparam(3); lambda_rigid  = fitparam(4);
// フィッティング結果
r0_rigid = -1 / lambda_rigid * log(- b_rigid / 2 / c_rigid);
lambda_rigid;
D_rigid = b_rigid ^ 2 / 4 / c_rigid;
A_rigid = a_rigid;
x0_rigid = - b_rigid / 2 / c_rigid;
// デバイモデル
gamma0 = lambda_rigid * r0_rigid / 2;
B0 = - c_rigid * x0_rigid ^ 2 * lambda_rigid ^ 3 / (6 * %pi * log(x0_rigid)) * (ry / bohr ^ 3 * 1E-8);
ThetaD0 = 41.63 * sqrt(r0_rigid * B0 / M);

// *** 絶対零度のヘルムホルツエネルギー ***
// 電子系 + ゼロ点エネルギー
F_0K = E + Ezero(rws);

// *** ゼロ点振動 ***
Energy = F_0K;
[fitparam, v] = lsqrsolve([a_rigid;b_rigid;c_rigid;lambda_rigid], fit, size(rws,1));
a_vec(1) = fitparam(1); b_vec(1) = fitparam(2); c_vec(1) = fitparam(3); lambda_vec(1) = fitparam(4);
// フィッティング結果
r0_vec(1) = -1 / lambda_vec(1) * log(- b_vec(1) / 2 / c_vec(1));
lambda_vec(1);
D_vec(1) = b_vec(1) ^ 2 / 4 / c_vec(1);
A_vec(1) = a_vec(1);
x0_vec(1) = - b_vec(1) / 2 / c_vec(1);
// デバイモデル
gamma0_vec(1) = lambda_vec(1) * r0_vec(1) / 2;
B0_vec(1) = - c_vec(1) * x0_vec(1) ^ 2 * lambda_vec(1) ^ 3 / (6 * %pi * log(x0_vec(1))) * (ry / bohr ^ 3 * 1E-8);
ThetaD0_vec(1) = 41.63 * sqrt(r0_vec(1) * B0 / M);

// *** 有限温度の計算 ***
// 各温度について計算
for i = 2:1:length(temperature) do
  T = temperature(i);
  Energy = E - kBry .* T .* (Debye(rws, T) - 3 * log(1 - exp(- ThetaD(rws) ./ T))) + Ezero(rws);
  exec('finit.sce');
  // フィッティング結果
  a_vec(i) = fitparam(1); b_vec(i) = fitparam(2); c_vec(i) = fitparam(3); lambda_vec(i) = fitparam(4);
  r0_vec(i) = -1 / lambda_vec(i) * log(- b_vec(i) / 2 / c_vec(i));
  lambda_vec(i);
  D_vec(i) = b_vec(i) ^ 2 / 4 / c_vec(i);
  A_vec(i) = a_vec(i);
  x0_vec(i) = - b_vec(i) / 2 / c_vec(i);
  // デバイモデル
  gamma0_vec(i) = lambda_vec(i) * r0_vec(i) / 2;
  B0_vec(i) = - c_vec(i) * x0_vec(i) ^ 2 * lambda_vec(i) ^ 3 / (6 * %pi * log(x0_vec(i))) * (ry / bohr ^ 3 * 1E-8);
  ThetaD0_vec(i) = 41.63 * sqrt(r0_vec(i) * B0 / M);
end

// *** グラフのプロット ***
scf(0);
// rigid latticeのプロット
a = a_rigid; b = b_rigid; c = c_rigid; lambda = lambda_rigid;
Eref = Morse(r0_rigid);
plot(rplot, Morse(rplot) - Eref, '--k');
plot(r0_rigid, Morse(r0_rigid) - Eref, '.k');
// ゼロ点振動を含むプロット
a = a_vec(1); b = b_vec(1); c = c_vec(1); lambda = lambda_vec(1);
plot(rplot, Morse(rplot) - Eref, '-b');
plot(r0_vec(1), Morse(r0_vec(1)) - Eref, '.b');
// 100Kのプロット
i100K = find(temperature==100);
if length(i100K) == 1 then
  a = a_vec(i100K); b = b_vec(i100K); c = c_vec(i100K); lambda = lambda_vec(i100K);
  plot(rplot, Morse(rplot) - Eref, '-c');
  plot(r0_vec(i100K), Morse(r0_vec(i100K)) - Eref, '.c');
end
// 200Kのプロット
i200K = find(temperature==200);
if length(i200K) == 1 then
  a = a_vec(i200K); b = b_vec(i200K); c = c_vec(i200K); lambda = lambda_vec(i200K);
  plot(rplot, Morse(rplot) - Eref, '-g');
  plot(r0_vec(i200K), Morse(r0_vec(i200K)) - Eref, '.g');
end
// 300Kのプロット
i300K = find(temperature==300);
if length(i300K) == 1 then
  a = a_vec(i300K); b = b_vec(i300K); c = c_vec(i300K); lambda = lambda_vec(i300K);
  plot(rplot, Morse(rplot) - Eref, '-r');
  plot(r0_vec(i300K), Morse(r0_vec(i300K)) - Eref, '.r');
end
// グラフの装飾
xlabel("Wigner-Seitz Radius (Bohr)");
ylabel("Free Energy (Ry)");
// 描画領域の制限
xmin=4.6; xmax=5.1; ymin=-10E-3; ymax=3E-3;
zoom_rect([xmin,ymin,xmax,ymax]);
xgrid(color("gray"));

// *** 温度依存性のプロット ***
scf(1);
// 熱膨張率
subplot(3,1,1);
// 前進差分
// alpha_f = (r0_vec(2:$) ./ r0_vec(1:$-1) - 1) ./ tstep;
// plot(temperature(1:$-1), 1E6 * alpha_f, '-b');
// 中心差分
alpha_c = (r0_vec(3:$) - r0_vec(1:$-2)) ./ (2 * tstep) ./ r0_vec(2:$-1);
plot(temperature(2:$-1), 1E6 * alpha_c, '-b');
// 中心差分のRomberg1段
// alpha_cr1 = 2 * ((r0_vec(4:$-1) - r0_vec(2:$-3)) ./ (2 * tstep) - 0.5 * (r0_vec(5:$) - r0_vec(1:$-4)) ./ (4 * tstep)) ./ r0_vec(3:$-2);
// plot(temperature(3:$-2), 1E6 * alpha_cr1, '-b');
// グラフの装飾
xgrid(color("gray"));
xlabel("Temperature (K)");
ylabel("alpha (10^-6/K)");

// Wigner-Seitz半径
subplot(3,1,2);
plot(temperature(1:$-1), r0_vec(1:$-1) ./ r0_vec(1));
// グラフの装飾
xgrid(color("gray"));
xlabel("Temperature (K)");
ylabel("r0/r0(T=0) (P=0)");

// 体積弾性率
subplot(3,1,3);
plot(temperature(1:$-1), B0_vec(1:$-1) ./ B0_vec(1));
// グラフの装飾
xgrid(color("gray"));
xlabel("Temperature (K)");
ylabel("B0/B0(T=0) (P=0)");