1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/tests/shootout/nbody.py Sun Feb 18 01:39:11 2007 +0100
1.3 @@ -0,0 +1,131 @@
1.4 +#!/usr/bin/python -OO
1.5 +# The Computer Language Shootout Benchmarks
1.6 +# http://shootout.alioth.debian.org/
1.7 +#
1.8 +# contributed by Kevin Carson
1.9 +
1.10 +import sys
1.11 +
1.12 +pi = 3.14159265358979323
1.13 +solar_mass = 4 * pi * pi
1.14 +days_per_year = 365.24
1.15 +
1.16 +class body :
1.17 + pass
1.18 +
1.19 +def advance(bodies, dt) :
1.20 + for i in xrange(len(bodies)) :
1.21 + b = bodies[i]
1.22 +
1.23 + for j in xrange(i + 1, len(bodies)) :
1.24 + b2 = bodies[j]
1.25 +
1.26 + dx = b.x - b2.x
1.27 + dy = b.y - b2.y
1.28 + dz = b.z - b2.z
1.29 + distance = (dx**2 + dy**2 + dz**2)**0.5
1.30 +
1.31 + b_mass_x_mag = dt * b.mass / distance**3
1.32 + b2_mass_x_mag = dt * b2.mass / distance**3
1.33 +
1.34 + b.vx -= dx * b2_mass_x_mag
1.35 + b.vy -= dy * b2_mass_x_mag
1.36 + b.vz -= dz * b2_mass_x_mag
1.37 + b2.vx += dx * b_mass_x_mag
1.38 + b2.vy += dy * b_mass_x_mag
1.39 + b2.vz += dz * b_mass_x_mag
1.40 +
1.41 + for b in bodies :
1.42 + b.x += dt * b.vx
1.43 + b.y += dt * b.vy
1.44 + b.z += dt * b.vz
1.45 +
1.46 +def energy(bodies) :
1.47 + e = 0.0
1.48 + for i in xrange(len(bodies)) :
1.49 + b = bodies[i]
1.50 + e += 0.5 * b.mass * (b.vx**2 + b.vy**2 + b.vz**2)
1.51 +
1.52 + for j in xrange(i + 1, len(bodies)) :
1.53 + b2 = bodies[j]
1.54 +
1.55 + dx = b.x - b2.x
1.56 + dy = b.y - b2.y
1.57 + dz = b.z - b2.z
1.58 + distance = (dx**2 + dy**2 + dz**2)**0.5
1.59 +
1.60 + e -= (b.mass * b2.mass) / distance
1.61 +
1.62 + return e
1.63 +
1.64 +def offset_momentum(bodies) :
1.65 + global sun
1.66 + px = py = pz = 0.0
1.67 +
1.68 + for b in bodies :
1.69 + px += b.vx * b.mass
1.70 + py += b.vy * b.mass
1.71 + pz += b.vz * b.mass
1.72 +
1.73 + sun.vx = - px / solar_mass
1.74 + sun.vy = - py / solar_mass
1.75 + sun.vz = - pz / solar_mass
1.76 +
1.77 +sun = body()
1.78 +sun.x = sun.y = sun.z = sun.vx = sun.vy = sun.vz = 0.0
1.79 +sun.mass = solar_mass
1.80 +
1.81 +jupiter = body()
1.82 +jupiter.x = 4.84143144246472090e+00
1.83 +jupiter.y = -1.16032004402742839e+00
1.84 +jupiter.z = -1.03622044471123109e-01
1.85 +jupiter.vx = 1.66007664274403694e-03 * days_per_year
1.86 +jupiter.vy = 7.69901118419740425e-03 * days_per_year
1.87 +jupiter.vz = -6.90460016972063023e-05 * days_per_year
1.88 +jupiter.mass = 9.54791938424326609e-04 * solar_mass
1.89 +
1.90 +saturn = body()
1.91 +saturn.x = 8.34336671824457987e+00
1.92 +saturn.y = 4.12479856412430479e+00
1.93 +saturn.z = -4.03523417114321381e-01
1.94 +saturn.vx = -2.76742510726862411e-03 * days_per_year
1.95 +saturn.vy = 4.99852801234917238e-03 * days_per_year
1.96 +saturn.vz = 2.30417297573763929e-05 * days_per_year
1.97 +saturn.mass = 2.85885980666130812e-04 * solar_mass
1.98 +
1.99 +uranus = body()
1.100 +uranus.x = 1.28943695621391310e+01
1.101 +uranus.y = -1.51111514016986312e+01
1.102 +uranus.z = -2.23307578892655734e-01
1.103 +uranus.vx = 2.96460137564761618e-03 * days_per_year
1.104 +uranus.vy = 2.37847173959480950e-03 * days_per_year
1.105 +uranus.vz = -2.96589568540237556e-05 * days_per_year
1.106 +uranus.mass = 4.36624404335156298e-05 * solar_mass
1.107 +
1.108 +neptune = body()
1.109 +neptune.x = 1.53796971148509165e+01
1.110 +neptune.y = -2.59193146099879641e+01
1.111 +neptune.z = 1.79258772950371181e-01
1.112 +neptune.vx = 2.68067772490389322e-03 * days_per_year
1.113 +neptune.vy = 1.62824170038242295e-03 * days_per_year
1.114 +neptune.vz = -9.51592254519715870e-05 * days_per_year
1.115 +neptune.mass = 5.15138902046611451e-05 * solar_mass
1.116 +
1.117 +def main() :
1.118 + try :
1.119 + n = int(sys.argv[1])
1.120 + except :
1.121 + print "Usage: %s <N>" % sys.argv[0]
1.122 +
1.123 + bodies = [sun, jupiter, saturn, uranus, neptune]
1.124 +
1.125 + offset_momentum(bodies)
1.126 +
1.127 + print "%.9f" % energy(bodies)
1.128 +
1.129 + for i in xrange(n) :
1.130 + advance(bodies, 0.01)
1.131 +
1.132 + print "%.9f" % energy(bodies)
1.133 +
1.134 +main()