clear;

// 体心正方(bct)の格子体積を固定
v = 100;
// bctのc/aを変化
// c/a = 1:   bcc
// c/a = √2: fcc
//eta = linspace(1, sqrt(2));
eta = linspace(1, 1.5);

// 格子定数aの計算
a = (v ./ eta) .^ (1/3);

// ab面内で頂点同士が触れる時
rmt_ab = a ./ 2;
// // ac面内で頂点同士が触れる時
// rmt_ac = eta .* a ./ 2;
// 頂点と体心で触れる時
rmt_diag = sqrt(2 .* (a .^ 2) + (eta .^ 2) .* (a .^ 2)) ./ 4;

// *** グラフの描画 ***
// ab面内
plot(eta, rmt_ab, '-b');
// // ac面内
// plot(eta, rmt_ac, '-g');
// 頂点と体心
plot(eta, rmt_diag, '-r');
// グラフの装飾
xlabel("c/a");
ylabel("muffin-tin radius");