Ising model

library(lattice)
range_f = function(X, loc) c(X[loc[1], c(loc[2]-1, loc[2]+1)], X[c(loc[1]-1, loc[1]+1), loc[2]])
ising_model = function(N, beta, n, range_f){
  X = replicate(N+2, rbinom(N+2, 1, .5))
  for (i in 1:n)
  {
      X[c(1, N+2),] = X[c(N+1, 2),]
      X[,c(1, N+2)] = X[,c(N+1, 2)]
      loc = floor(N*runif(2)) + 2
      if (runif(1) < exp(2*beta*(2-sum(range_f(X, loc) == X[loc[1], loc[2]]))))
        X[loc[1],loc[2]] = 1 - X[loc[1],loc[2]]
      if(i %% 50 == 0)
      {
        print(levelplot(X, col.regions = rgb(0,0,0,seq(0,1,length=16))))
      }
  }
}
library(animation)
ani.options(ffmpeg = "C:/ffmpeg/bin/ffmpeg.exe", interval = 1/25)
saveVideo(ising_model(100, -1, 100000, range_f), "beta_is_-1.mp4")
saveVideo(ising_model(100, 2, 100000, range_f), "beta_is_2.mp4")
saveVideo(ising_model(100, 6, 100000, range_f), "beta_is_6.mp4")
saveVideo(ising_model(100, 10, 100000, range_f), "beta_is_10.mp4")

ani.options(convert ="C:/Program Files/ImageMagick-6.9.1-Q16/convert.exe")
saveGIF(ising_model(100, -1, 100000, range_f), "beta_is_-1.gif")
saveGIF(ising_model(100, 2, 100000, range_f), "beta_is_2.gif")
saveGIF(ising_model(100, 6, 100000, range_f), "beta_is_6.gif")
saveGIF(ising_model(100, 10, 100000, range_f), "beta_is_10.gif")