bench/yarv/bm_so_nbody.rb

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