benchmark/bm_so_nbody.rb

   1 # The Computer Language Shootout
   2 # http://shootout.alioth.debian.org
   3 #
   4 # Optimized for Ruby by Jesse Millikan
   5 # From version ported by Michael Neumann from the C gcc version,
   6 # which was written by Christoph Bauer.
   7
   8 SOLAR_MASS = 4 * Math::PI**2
   9 DAYS_PER_YEAR = 365.24
  10
  11 def _puts *args
  12 end
  13
  14 class Planet
  15  attr_accessor :x, :y, :z, :vx, :vy, :vz, :mass
  16
  17  def initialize(x, y, z, vx, vy, vz, mass)
  18   @x, @y, @z = x, y, z
  19   @vx, @vy, @vz = vx * DAYS_PER_YEAR, vy * DAYS_PER_YEAR, vz * DAYS_PER_YEAR
  20   @mass = mass * SOLAR_MASS
  21  end
  22
  23  def move_from_i(bodies, nbodies, dt, i)
  24   while i < nbodies
  25    b2 = bodies[i]
  26    dx = @x - b2.x
  27    dy = @y - b2.y
  28    dz = @z - b2.z
  29
  30    distance = Math.sqrt(dx * dx + dy * dy + dz * dz)
  31    mag = dt / (distance * distance * distance)
  32    b_mass_mag, b2_mass_mag = @mass * mag, b2.mass * mag
  33
  34    @vx -= dx * b2_mass_mag
  35    @vy -= dy * b2_mass_mag
  36    @vz -= dz * b2_mass_mag
  37    b2.vx += dx * b_mass_mag
  38    b2.vy += dy * b_mass_mag
  39    b2.vz += dz * b_mass_mag
  40    i += 1
  41   end
  42
  43   @x += dt * @vx
  44   @y += dt * @vy
  45   @z += dt * @vz
  46  end
  47 end
  48
  49 def energy(bodies)
  50   e = 0.0
  51   nbodies = bodies.size
  52
  53   for i in 0 ... nbodies
  54     b = bodies[i]
  55     e += 0.5 * b.mass * (b.vx * b.vx + b.vy * b.vy + b.vz * b.vz)
  56     for j in (i + 1) ... nbodies
  57       b2 = bodies[j]
  58       dx = b.x - b2.x
  59       dy = b.y - b2.y
  60       dz = b.z - b2.z
  61       distance = Math.sqrt(dx * dx + dy * dy + dz * dz)
  62       e -= (b.mass * b2.mass) / distance
  63     end
  64   end
  65   e
  66 end
  67
  68 def offset_momentum(bodies)
  69   px, py, pz = 0.0, 0.0, 0.0
  70
  71   for b in bodies
  72     m = b.mass
  73     px += b.vx * m
  74     py += b.vy * m
  75     pz += b.vz * m
  76   end
  77
  78   b = bodies[0]
  79   b.vx = - px / SOLAR_MASS
  80   b.vy = - py / SOLAR_MASS
  81   b.vz = - pz / SOLAR_MASS
  82 end
  83
  84 BODIES = [
  85   # sun
  86   Planet.new(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0),
  87
  88   # jupiter
  89   Planet.new(
  90     4.84143144246472090e+00,
  91     -1.16032004402742839e+00,
  92     -1.03622044471123109e-01,
  93     1.66007664274403694e-03,
  94     7.69901118419740425e-03,
  95     -6.90460016972063023e-05,
  96     9.54791938424326609e-04),
  97
  98   # saturn
  99   Planet.new(
 100     8.34336671824457987e+00,
 101     4.12479856412430479e+00,
 102     -4.03523417114321381e-01,
 103     -2.76742510726862411e-03,
 104     4.99852801234917238e-03,
 105     2.30417297573763929e-05,
 106     2.85885980666130812e-04),
 107
 108   # uranus
 109   Planet.new(
 110     1.28943695621391310e+01,
 111     -1.51111514016986312e+01,
 112     -2.23307578892655734e-01,
 113     2.96460137564761618e-03,
 114     2.37847173959480950e-03,
 115     -2.96589568540237556e-05,
 116     4.36624404335156298e-05),
 117
 118   # neptune
 119   Planet.new(
 120     1.53796971148509165e+01,
 121     -2.59193146099879641e+01,
 122     1.79258772950371181e-01,
 123     2.68067772490389322e-03,
 124     1.62824170038242295e-03,
 125     -9.51592254519715870e-05,
 126     5.15138902046611451e-05)
 127 ]
 128
 129 init = 200_000 # ARGV[0]
 130 n = Integer(init)
 131
 132 offset_momentum(BODIES)
 133
 134 puts "%.9f" % energy(BODIES)
 135
 136 nbodies = BODIES.size
 137 dt = 0.01
 138
 139 n.times do
 140   i = 0
 141   while i < nbodies
 142     b = BODIES[i]
 143     b.move_from_i(BODIES, nbodies, dt, i + 1)
 144     i += 1
 145   end
 146 end
 147
 148 puts "%.9f" % energy(BODIES)