// *** リュードベリ(Rydberg)原子単位系 ***
hbar = 1;		// 換算プランク定数
h = 2 * %pi * hbar;	// プランク定数
m = 0.5;		// 電子の質量
a0 = 1;			// ボーア半径

// *** 系の設定 ***
// 原子の数
//n_atom = 20;
n_atom = 4 * gmax + 2;
// 原子間距離
r0 = 1.0 * a0;
// 結晶の長さ
l = n_atom * r0;
// 価数と電子の総数
n_valence = 1;
n_electron = n_valence *  n_atom;

// *** 分散関係の計算 ***
g = 0;
// 波数ベクトル (式 3.11)
kg = 2 * %pi / l * g;
Kg = [kg];
// エネルギー (式 3.13)
eg = hbar * hbar / (2 * m) * (2 * %pi / l * g) * (2 * %pi / l * g);
Eg = [eg];
for g=1:int((n_electron - 2) / 4)
  // 波数ベクトル (式 3.11)
  kg = 2 * %pi / l * g;
  Kg = [-kg, Kg, kg];
  // エネルギー (式 3.13)
  eg = hbar * hbar / (2 * m) * (2 * %pi / l * g) * (2 * %pi / l * g);
  Eg = [eg, Eg, eg];
end
// フェルミ波数ベクトル (1/Bohr)
kf = kg;
// フェルミエネルギー (Ry)
ef = eg;

// // *** グラフの描画 ***
// plot(0, 0); // グラフ描画範囲のためのダミー
// plot(Kg, Eg, ".b");
//
// // *** グラフの装飾 ***
// xgrid(color("gray"));
// xlabel("Wave vector (1/Bohr)");
// ylabel("Energy (Ry)");
// //zoom_rect([-kf, 0, kf, ef]);