clear;

// *** 定数の設定 ***
n = 100;		// 粒子の数
m = 200;		// 熱力学的な平均を取る回数
J = 1;			// 交換エネルギー
rand("uniform");	// 乱数は一様乱数とする
tmax = 5 * n;		// 時間の最大ステップ
h = 0.0;		// 外部磁場
// 温度
ktmin = 0.5;		// 最低温度
ktmax = 5.0;		// 最高温度
nkt = 19;		// 温度の分割数
//T = linspace(ktmin, ktmax, nkt);	// 低温から開始
//spin = ones(1,n);			// 各粒子におけるスピン(コールドスタート)
T = linspace(ktmax, ktmin, nkt);	// 高温から開始
spin = 1 - 2 * round(rand(1,n));	// 各粒子におけるスピン(ランダム)

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

// *** 行列の初期化 ***
E  = [];	// エネルギーの和
E2 = [];	// エネルギーの二乗の和
M  = [];	// 磁化の和
M2 = [];	// 磁化の二乗の和

// *** 温度のループ ***
for kt = 1:nkt do
  // エネルギーの初期化
  ene1 = 0;			// エネルギーの和
  ene2 = 0;			// エネルギーの二乗和
  // 磁化の初期化
  mag1 = 0;			// 磁化の和
  mag2 = 0;			// 磁化の二乗和f
  // *** 熱力学平均のループ ***
  for samp = 1:m do
    // *** 時間発展のループ ***
    for t = 1:tmax do
      oldenergy = energy(spin);
      element = ceil(n * rand());		// 粒子を一つ選ぶ
      spin(element) = -1 * spin(element);	// スピンを反転
      newenergy = energy(spin);
      if (newenergy > oldenergy) & (exp((- newenergy + oldenergy) / T(kt)) < rand()) then
        spin(element) = -1 * spin(element);	// 棄却
      end
    end
    ene1 = ene1 + energy(spin);		// エネルギーの和
    ene2 = ene2 + energy(spin)^2;	// エネルギーの二乗の和
    mag1 = mag1 + sum(sum(spin));	// 磁化の和
    mag2 = mag2 + sum(sum(spin))^2;	// 磁化の二乗和
  end
  E  = [E,  ene1 / m];			// エネルギーの和
  E2 = [E2, ene2 / m];			// エネルギーの二乗の和
  M  = [M,  mag1 / m];			// 磁化の和
  M2 = [M2, mag2 / m];			// 磁化の二乗和
end

// *** エネルギーと磁化の揺らぎ ***
C = (E2 - E .^ 2) ./ (n * T .^ 2);	// 比熱
X = (M2 - M .^ 2) ./ (n * T);		// 磁化率

// *** 厳密解の計算 ***
// 温度ベクトル
Ta = linspace(0.1,5,50);
// 粒子1個あたりの平均エネルギー
Ea = - tanh(J ./ Ta);
// 比熱
Ca = (J ./ Ta) .^ 2 ./ cosh(J ./ Ta) .^ 2;
// 磁化
Ma = sinh(h ./ Ta) ./ sqrt(sinh(h ./ Ta) .^ 2 + exp(-4 * J ./ Ta));
// 磁化率
Xa = exp(2 * J ./ Ta) ./ Ta;
RXa = 1 ./ Xa;

// *** グラフのプロット ***
// エネルギー
subplot(2,2,1);
plot(T, E ./ n, 'or');
plot(Ta, Ea, '--g');
xlabel("kT/J");
ylabel("E/NJ");
// 比熱
subplot(2,2,2);
plot(T, C, 'or');
plot(Ta,Ca,'--g');
xlabel("kT/J");
ylabel("C/Nk");
// 磁化
subplot(2,2,3);
plot(T, M ./ n, 'or');
plot(Ta,Ma ./ n,'--g');
xlabel("kT/J");
ylabel("M/N");
// 磁化率
subplot(2,2,4);
plot(T, 1 ./ X, 'or');
plot(Ta,RXa,'--g');
xlabel("kT/J");
ylabel("N/JX");