clear;

// *** 定数の設定 ***
n = 100;		// 粒子の数
kt = 1.5;		// 温度
J = 1;			// 交換エネルギー (1: 強磁性, -1:反強磁性)
rand("uniform");	// 乱数は一様乱数とする
tmax = 1000;		// 時間の最大ステップ

// *** 初期化 ***
// 各粒子におけるスピン
spin = ones(1,n);			// コールドスタート
//spin = 1 - 2 * round(rand(1,n));	// 乱数スタート

// *** エネルギーの計算関数 ***
function e = energy(spin)
  e = - J * sum(spin .* [spin(2:n), spin(1)]);
endfunction

SPIN = [];

// *** 時間発展 ***
for t = 1:tmax do
  oldenergy = energy(spin);
  element = ceil(n * rand());		// 粒子を一つ選ぶ
  spin(element) = -1 * spin(element);	// スピンを反転
  newenergy = energy(spin);
  spin(element) = (- 2 * ((newenergy > oldenergy) & (exp((- newenergy + oldenergy) / kt) <= rand())) + 1) * spin(element);
  // 時刻tにおけるスピンを保存
  SPIN = [SPIN;spin];
end

// *** スピンの時間変化をプロット ***
Matplot((SPIN' + 1) .* 10);
zoom_rect([0,0,tmax,n]);
xlabel("Time");
ylabel("Position");