clear;
format('v',16)
bohr = 0.5291772108 // 1 bohr = 0.529 A

// *** データの読み込み ***
e = fscanfMat("Energy.txt")';
// 第一原理計算の格子定数
a = [4.70:0.005:4.79];
c = [1.60:0.005:1.69];

// *** スプライン補間 ***
E = splin2d(a,c,e);
aa = [4.70:0.001:4.79];
cc = [1.60:0.001:1.69];
[AA,CC] = ndgrid(aa,cc);
ee = interp2d(AA, CC, a, c, E);

// *** カラーコンターのプロット ***
// グラフを正方形に
isoview(4.70,4.79,1.60,1.69);
// 色の設定
xset("colormap",jetcolormap(64))
// カラコンターのプロット
Sgrayplot(aa, cc, ee - min(ee));
// カラーバーの表示
colorbar(0, max(ee) - min(ee));
// グラフの軸ラベル
xlabel("a (bohr)");
ylabel("c/a");

// *** 平衡格子定数の計算 ***
// 補間データから最小値を探す
[emin,ind] = min(ee);
// 平衡格子定数のプロット(白丸)
plot(aa(ind(1)),cc(ind(2)),'ow');
// 文献値(実験データ)のプロット
// http://www.tuat.ac.jp/~sameken/lecture/2012raremetal.pdf
plot(2.514 / bohr,  4.105 / 2.514,'xr');	// 赤クロス
// https://github.com/cryos/avogadro/blob/master/crystals/elements/Co-Cobalt.cif
plot(2.5071 / bohr, 4.0686 / 2.5071,'xy');	// 黄クロス