clear;

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

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

// *** プロットと画像出力 ***
driver('GIF');				// GIF形式で出力
xinit('ising2d0000.GIF');		// 画像ファイルの指定
isoview(0,n,0,n)			// 画像を正方形にする
Matplot((spin + 1) .* 15);		// スピン状態をプロット
xend();					// 画像ファイルを閉じる

// *** エネルギーの計算関数 ***
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
  // ファイル番号の先頭にゼロを補う
  filenum = string(t);
  if length(filenum) == 1 then
    filenum = strcat(["000",filenum]);
  elseif length(filenum) == 2 then
    filenum = strcat(["00",filenum]);
  elseif length(filenum) == 3 then
    filenum = strcat(["0",filenum]);
  end
  // 画像ファイルの指定
  xinit(strcat(['ising2d',filenum,'.GIF']));

  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);

  // スピン状態をプロット
  isoview(0,n,0,n)			// 画像を正方形にする
  Matplot((spin + 1) .* 15);
  xend();
end

// グラフ出力先をコンピュータ画面出力へもどす
driver('X11');