clear;

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

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

// *** 分散関係の計算 ***
// 占有されている準位の計算
Kg = [];
Eg = [];
for g=1:int(n_electron / 2)
  // 波数ベクトル (式 2.8)
  kg = %pi / l * g;
  Kg = [Kg, kg];
  // エネルギー (式 2.18)
  eg = h * h / (8 * m * l * l) * g * g;
  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([0, 0, kf, ef]);