clear;

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

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

isoview(0,n,0,n)

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

// *** 時間発展 ***
for t = 1:tmax do
  oldenergy = energy(spin);
  // 粒子を一つ選ぶ
  elementx = ceil(n * rand());
  elementy = ceil(n * rand());
  // スピンを反転
  spin(elementx,elementy) = -1 * spin(elementx,elementy);
  newenergy = energy(spin);
  spin(elementx, elementy) = (- 2 * ((newenergy > oldenergy) & (exp((- newenergy + oldenergy) / kt) <= rand())) + 1) * spin(elementx, elementy);
  // スピン状態をプロット
  Matplot((spin + 1) .* 15);
end