1 #!/usr/bin/python -OO 2 # The Computer Language Shootout Benchmarks 3 # http://shootout.alioth.debian.org/ 4 # 5 # contributed by Kevin Carson 6 7 import sys 8 9 pi = 3.14159265358979323 10 solar_mass = 4 * pi * pi 11 days_per_year = 365.24 12 13 class body : 14 pass 15 16 def advance(bodies, dt) : 17 for i in xrange(len(bodies)) : 18 b = bodies[i] 19 20 for j in xrange(i + 1, len(bodies)) : 21 b2 = bodies[j] 22 23 dx = b.x - b2.x 24 dy = b.y - b2.y 25 dz = b.z - b2.z 26 distance = (dx**2 + dy**2 + dz**2)**0.5 27 28 b_mass_x_mag = dt * b.mass / distance**3 29 b2_mass_x_mag = dt * b2.mass / distance**3 30 31 b.vx -= dx * b2_mass_x_mag 32 b.vy -= dy * b2_mass_x_mag 33 b.vz -= dz * b2_mass_x_mag 34 b2.vx += dx * b_mass_x_mag 35 b2.vy += dy * b_mass_x_mag 36 b2.vz += dz * b_mass_x_mag 37 38 for b in bodies : 39 b.x += dt * b.vx 40 b.y += dt * b.vy 41 b.z += dt * b.vz 42 43 def energy(bodies) : 44 e = 0.0 45 for i in xrange(len(bodies)) : 46 b = bodies[i] 47 e += 0.5 * b.mass * (b.vx**2 + b.vy**2 + b.vz**2) 48 49 for j in xrange(i + 1, len(bodies)) : 50 b2 = bodies[j] 51 52 dx = b.x - b2.x 53 dy = b.y - b2.y 54 dz = b.z - b2.z 55 distance = (dx**2 + dy**2 + dz**2)**0.5 56 57 e -= (b.mass * b2.mass) / distance 58 59 return e 60 61 def offset_momentum(bodies) : 62 global sun 63 px = py = pz = 0.0 64 65 for b in bodies : 66 px += b.vx * b.mass 67 py += b.vy * b.mass 68 pz += b.vz * b.mass 69 70 sun.vx = - px / solar_mass 71 sun.vy = - py / solar_mass 72 sun.vz = - pz / solar_mass 73 74 sun = body() 75 sun.x = sun.y = sun.z = sun.vx = sun.vy = sun.vz = 0.0 76 sun.mass = solar_mass 77 78 jupiter = body() 79 jupiter.x = 4.84143144246472090e+00 80 jupiter.y = -1.16032004402742839e+00 81 jupiter.z = -1.03622044471123109e-01 82 jupiter.vx = 1.66007664274403694e-03 * days_per_year 83 jupiter.vy = 7.69901118419740425e-03 * days_per_year 84 jupiter.vz = -6.90460016972063023e-05 * days_per_year 85 jupiter.mass = 9.54791938424326609e-04 * solar_mass 86 87 saturn = body() 88 saturn.x = 8.34336671824457987e+00 89 saturn.y = 4.12479856412430479e+00 90 saturn.z = -4.03523417114321381e-01 91 saturn.vx = -2.76742510726862411e-03 * days_per_year 92 saturn.vy = 4.99852801234917238e-03 * days_per_year 93 saturn.vz = 2.30417297573763929e-05 * days_per_year 94 saturn.mass = 2.85885980666130812e-04 * solar_mass 95 96 uranus = body() 97 uranus.x = 1.28943695621391310e+01 98 uranus.y = -1.51111514016986312e+01 99 uranus.z = -2.23307578892655734e-01 100 uranus.vx = 2.96460137564761618e-03 * days_per_year 101 uranus.vy = 2.37847173959480950e-03 * days_per_year 102 uranus.vz = -2.96589568540237556e-05 * days_per_year 103 uranus.mass = 4.36624404335156298e-05 * solar_mass 104 105 neptune = body() 106 neptune.x = 1.53796971148509165e+01 107 neptune.y = -2.59193146099879641e+01 108 neptune.z = 1.79258772950371181e-01 109 neptune.vx = 2.68067772490389322e-03 * days_per_year 110 neptune.vy = 1.62824170038242295e-03 * days_per_year 111 neptune.vz = -9.51592254519715870e-05 * days_per_year 112 neptune.mass = 5.15138902046611451e-05 * solar_mass 113 114 def main() : 115 try : 116 n = int(sys.argv[1]) 117 except : 118 print "Usage: %s <N>" % sys.argv[0] 119 120 bodies = [sun, jupiter, saturn, uranus, neptune] 121 122 offset_momentum(bodies) 123 124 print "%.9f" % energy(bodies) 125 126 for i in xrange(n) : 127 advance(bodies, 0.01) 128 129 print "%.9f" % energy(bodies) 130 131 main()