ising.rb

J  = 1
Mu = 1

if ARGV.length != 2
  puts "Usage: #{File.basename(__FILE__)} <T> <H>"
  exit 1
else
  KbT = ARGV.shift.to_f
  H = ARGV.shift.to_f
end
  
sp = Array.new(50) { rand < 0.3 ? 1 : -1 }

def spins_display(spins)
  o = ""  
  spins.each do |s|
    if s == 1
      o << "⬆"
    else
      o << "\e[1;30m⬇\e[0m"
    end
  end
  o << " "
end

def flip_spins(spins)
  indices = (0...spins.length).to_a.sort_by { rand }
  indices.each do |i|
    s = spins[i]
    r = spins[i-1]
    t = ((i+1) == spins.length) ? spins[0] : spins[i+1]

    e_0 = (r *  s) + ( s * t) - (Mu * H *  s)
    e_1 = (r * -s) + (-s * t) - (Mu * H * -s)
    
    if (e_1 < e_0) || (Math.exp(-(e_1 - e_0)/KbT.to_f) > rand)
      spins[i] = -s
    end
  end
end

def energy(spins)
  interaction_energy = 0
  
  spins.each_cons(2) do |pair|
    interaction_energy += pair[0] * pair[1]
  end
  interaction_energy += spins[0] * spins[-1]
  
  interaction_energy *= J
  
  field_energy = - Mu * H * spins.inject(0) { |sum, x| sum += x }
    
  return interaction_energy + field_energy
end

def magnetization(spins)
  (Mu * spins.inject(0) { |sum, x| sum += x }) / spins.length.to_f
end

# all_up = Array.new(50) { 1 }
# all_down = Array.new(50) { -1 }
# even = Array.new(50) { |i| (i % 2) == 0 ? 1 : -1 }

# print spins_display(all_up)
# puts " #{energy(all_up)}"
# print spins_display(all_down)
# puts " #{energy(all_down)}"
# print spins_display(even)
# puts " #{energy(even)}"
# 
# puts "-" * 55

$stderr.print spins_display(sp)
puts "# E M" # energy, magnetization

loop do
  $stderr.print spins_display(sp)
  puts "%4i %5.2f" % [energy(sp), magnetization(sp)]
  
  flip_spins(sp)
end

puts