|
| 1 | +/* |
| 2 | + * The Computer Language Benchmarks Game |
| 3 | + * http://shootout.alioth.debian.org/ |
| 4 | + * |
| 5 | + * Translated from Christoph Bauer's nbody.c by Ian Martins |
| 6 | + */ |
| 7 | +typedef Planet = { |
| 8 | + var x:Float; |
| 9 | + var y:Float; |
| 10 | + var z:Float; |
| 11 | + var vx:Float; |
| 12 | + var vy:Float; |
| 13 | + var vz:Float; |
| 14 | + var mass:Float; |
| 15 | +} |
| 16 | + |
| 17 | +class App { |
| 18 | + private static var SOLAR_MASS = 4 * Math.PI * Math.PI; |
| 19 | + private static var DAYS_PER_YEAR = 365.24; |
| 20 | + private static var DT = 1e-2; |
| 21 | + private static var RECIP_DT = (1.0 / DT); |
| 22 | + |
| 23 | + private static function advance(bodies:Array<Planet>) { |
| 24 | + for (i in 0...bodies.length) { |
| 25 | + var b = bodies[i]; |
| 26 | + for (j in (i + 1)...bodies.length) { |
| 27 | + var b2 = bodies[j]; |
| 28 | + var dx = b.x - b2.x; |
| 29 | + var dy = b.y - b2.y; |
| 30 | + var dz = b.z - b2.z; |
| 31 | + var invDist = 1.0 / Math.sqrt(dx * dx + dy * dy + dz * dz); |
| 32 | + var mag = invDist * invDist * invDist; |
| 33 | + b.vx -= dx * b2.mass * mag; |
| 34 | + b.vy -= dy * b2.mass * mag; |
| 35 | + b.vz -= dz * b2.mass * mag; |
| 36 | + b2.vx += dx * b.mass * mag; |
| 37 | + b2.vy += dy * b.mass * mag; |
| 38 | + b2.vz += dz * b.mass * mag; |
| 39 | + } |
| 40 | + } |
| 41 | + for (b in bodies) { |
| 42 | + b.x += b.vx; |
| 43 | + b.y += b.vy; |
| 44 | + b.z += b.vz; |
| 45 | + } |
| 46 | + } |
| 47 | + |
| 48 | + private static function energy(bodies:Array<Planet>) { |
| 49 | + var e = 0.0; |
| 50 | + for (i in 0...bodies.length) { |
| 51 | + var b = bodies[i]; |
| 52 | + e += 0.5 * b.mass * (b.vx * b.vx + b.vy * b.vy + b.vz * b.vz); |
| 53 | + for (j in (i + 1)...bodies.length) { |
| 54 | + var b2 = bodies[j]; |
| 55 | + var dx = b.x - b2.x; |
| 56 | + var dy = b.y - b2.y; |
| 57 | + var dz = b.z - b2.z; |
| 58 | + var distance = Math.sqrt(dx * dx + dy * dy + dz * dz); |
| 59 | + e -= (b.mass * b2.mass) / distance; |
| 60 | + } |
| 61 | + } |
| 62 | + return e; |
| 63 | + } |
| 64 | + |
| 65 | + private static function offsetMomentum(bodies:Array<Planet>) { |
| 66 | + var px = 0.0, py = 0.0, pz = 0.0; |
| 67 | + for (b in bodies) { |
| 68 | + px += b.vx * b.mass; |
| 69 | + py += b.vy * b.mass; |
| 70 | + pz += b.vz * b.mass; |
| 71 | + } |
| 72 | + bodies[0].vx = -px / SOLAR_MASS; |
| 73 | + bodies[0].vy = -py / SOLAR_MASS; |
| 74 | + bodies[0].vz = -pz / SOLAR_MASS; |
| 75 | + } |
| 76 | + |
| 77 | + /* |
| 78 | + * Rescale certain properties of bodies. That allows doing |
| 79 | + * consequential advance()'s as if dt were equal to 1.0. |
| 80 | + * |
| 81 | + * When all advances done, rescale bodies back to obtain correct energy. |
| 82 | + */ |
| 83 | + private static function scaleBodies(bodies:Array<Planet>, scale:Float) { |
| 84 | + for (b in bodies) { |
| 85 | + b.mass *= scale * scale; |
| 86 | + b.vx *= scale; |
| 87 | + b.vy *= scale; |
| 88 | + b.vz *= scale; |
| 89 | + } |
| 90 | + } |
| 91 | + |
| 92 | + inline private static function round(val:Float) { |
| 93 | + return Math.round(val * 1e9) / 1e9; |
| 94 | + } |
| 95 | + |
| 96 | + public static function main() { |
| 97 | + var bodies = new Array<Planet>(); |
| 98 | + bodies.push({ |
| 99 | + x: 0, // sun |
| 100 | + y: 0, |
| 101 | + z: 0, |
| 102 | + vx: 0, |
| 103 | + vy: 0, |
| 104 | + vz: 0, |
| 105 | + mass: SOLAR_MASS |
| 106 | + }); |
| 107 | + bodies.push({ |
| 108 | + x: 4.84143144246472090e+00, // jupiter |
| 109 | + y: -1.16032004402742839e+00, |
| 110 | + z: -1.03622044471123109e-01, |
| 111 | + vx: 1.66007664274403694e-03 * DAYS_PER_YEAR, |
| 112 | + vy: 7.69901118419740425e-03 * DAYS_PER_YEAR, |
| 113 | + vz: -6.90460016972063023e-05 * DAYS_PER_YEAR, |
| 114 | + mass: 9.54791938424326609e-04 * SOLAR_MASS |
| 115 | + }); |
| 116 | + bodies.push({ |
| 117 | + x: 8.34336671824457987e+00, // saturn |
| 118 | + y: 4.12479856412430479e+00, |
| 119 | + z: -4.03523417114321381e-01, |
| 120 | + vx: -2.76742510726862411e-03 * DAYS_PER_YEAR, |
| 121 | + vy: 4.99852801234917238e-03 * DAYS_PER_YEAR, |
| 122 | + vz: 2.30417297573763929e-05 * DAYS_PER_YEAR, |
| 123 | + mass: 2.85885980666130812e-04 * SOLAR_MASS |
| 124 | + }); |
| 125 | + bodies.push({ |
| 126 | + x: 1.28943695621391310e+01, // uranus |
| 127 | + y: -1.51111514016986312e+01, |
| 128 | + z: -2.23307578892655734e-01, |
| 129 | + vx: 2.96460137564761618e-03 * DAYS_PER_YEAR, |
| 130 | + vy: 2.37847173959480950e-03 * DAYS_PER_YEAR, |
| 131 | + vz: -2.96589568540237556e-05 * DAYS_PER_YEAR, |
| 132 | + mass: 4.36624404335156298e-05 * SOLAR_MASS |
| 133 | + }); |
| 134 | + bodies.push({ |
| 135 | + x: 1.53796971148509165e+01, // neptune |
| 136 | + y: -2.59193146099879641e+01, |
| 137 | + z: 1.79258772950371181e-01, |
| 138 | + vx: 2.68067772490389322e-03 * DAYS_PER_YEAR, |
| 139 | + vy: 1.62824170038242295e-03 * DAYS_PER_YEAR, |
| 140 | + vz: -9.51592254519715870e-05 * DAYS_PER_YEAR, |
| 141 | + mass: 5.15138902046611451e-05 * SOLAR_MASS |
| 142 | + }); |
| 143 | + |
| 144 | + var n = Std.parseInt(Sys.args()[0]); |
| 145 | + offsetMomentum(bodies); |
| 146 | + Sys.println(round(energy(bodies))); |
| 147 | + scaleBodies(bodies, DT); |
| 148 | + for (i in 0...n) |
| 149 | + advance(bodies); |
| 150 | + scaleBodies(bodies, RECIP_DT); |
| 151 | + Sys.println(round(energy(bodies))); |
| 152 | + } |
| 153 | +} |
0 commit comments