// *** フィッティング関数の定義 ***
// Morse関数
function E = Morse(r)
  E = a + b * exp(- lambda * r) + c * exp(-2 * lambda * r);
endfunction
// フィッティング関数
function err = fit(fitparam,m)
  a      = fitparam(1);
  b      = fitparam(2);
  c      = fitparam(3);
  lambda = fitparam(4);
  err    = Energy - Morse(rws);
endfunction

// *** ヘルムホルツの自由エネルギー ***
// デバイ温度の体積変化
function ThetaD = ThetaD(r)
  ThetaD = ThetaD0 * (r0_rigid ./ r) .^ (3 * gamma0)
endfunction
// Debye関数
function Debye = Debye(r, T)
//  Debye = 3 * (ThetaD(r) ./ T) .^ 3 .* integrate('(z .^ 4 * exp(z) ./ (exp(z) - 1) .^ 2)', 'z', 0, (ThetaD(r) ./ T))
//  Debye = 3 * (T ./ ThetaD(r)) .^ 3 .* integrate('(z .^ 4 * exp(z) ./ (exp(z) - 1) .^ 2)', 'z', 0, (ThetaD(r) ./ T))
  Debye = 3 * (T ./ ThetaD(r)) .^ 3 .* integrate('(z .^ 3 ./ (exp(z) - 1))', 'z', 0, (ThetaD(r) ./ T))
endfunction
// ゼロ点エネルギー
function Ezero = Ezero(r)
  Ezero = 9 / 8 * kBry * ThetaD(r)
endfunction
// 格子系のエネルギー
function Edebye = Edebye(r, T)
  Edebye = 3 * kBry .* T .* Debye(r, T) + Ezero(r)
endfunction
// 格子系のエントロピー
function Sdebye = Sdebye(r, T)
  Sdebye = 3 * kBry * (4 / 3 * Debye(r, T) - log(1 - exp(- ThetaD(r) ./ T)))
endfunction