Float pi = 3.141592653589793; Float solar_mass = 4.0 * pi * pi; Float days_per_year = 365.24; type Vector = vector(x: Float, y: Float, z: Float); Vector new_velocity(Nat body_idx, Vector velocity, Vector* positions, Float* masses, Float dt) { position = positions(body_idx); x = position.x; y = position.y; z = position.z; dvx = 0.0; dvy = 0.0; dvz = 0.0; for p @ i <- positions: if i != body_idx: dx = x - p.x; dy = y - p.y; dz = z - p.z; distance = sqrt(dx * dx + dy * dy + dz * dz); a_quantity = (masses(i) * dt) / (distance * distance * distance); dvx = dvx - dx * a_quantity; dvy = dvy - dy * a_quantity; dvz = dvz - dz * a_quantity; ; ; return vector(x: velocity.x + dvx, y: velocity.y + dvy, z: velocity.z + dvz); } (Vector*, Vector*) advance(Nat steps, Vector* positions, Vector* velocities, Float* masses, Float dt) { curr_positions = positions; curr_velocities = velocities; for steps: curr_velocities = (new_velocity(i, v, curr_positions, masses, dt) : v @ i <- curr_velocities); curr_positions = ({ velocity = curr_velocities(i); return vector( x: p.x + dt * velocity.x, y: p.y + dt * velocity.y, z: p.z + dt * velocity.z ); } : p @ i <- curr_positions ); ; return (curr_positions, curr_velocities); } Float energy(Vector* positions, Vector* velocities, Float* masses) { n = |positions|; e = 0.0; for i1 < n: p1 = positions(i1); v1 = velocities(i1); m1 = masses(i1); e = e + 0.5 * m1 * (v1.x * v1.x + v1.y * v1.y + v1.z * v1.z); for i2 = i1+1 .. n: p2 = positions(i2); dx = p1.x - p2.x; dy = p1.y - p2.y; dz = p1.z - p2.z; distance = sqrt(dx * dx + dy * dy + dz * dz); e = e - (masses(i1) * masses(i2)) / distance; ; ; return e; } Vector* offset_momentum(Vector* velocities, Float* masses) { px = 0.0; py = 0.0; pz = 0.0; for v @ i <- velocities: m = masses(i); px = px + v.x * m; py = py + v.y * m; pz = pz + v.z * m; ; v0 = velocities(0); v0 = vector( x: v0.x - px / solar_mass, y: v0.y - py / solar_mass, z: v0.z - pz / solar_mass ); return (if i == 0 then v0 else v : v @ i <- velocities); } (Float, Float, Float, Float, Float, Float, Float)* bodies = ( // Sun (0.0, 0.0, 0.0, 0.0, 0.0, 0.0, solar_mass), // Jupiter ( 4.84143144246472090e00, -1.16032004402742839e00, -1.03622044471123109e-01, 1.66007664274403694e-03 * days_per_year, 7.69901118419740425e-03 * days_per_year, -6.90460016972063023e-05 * days_per_year, 9.54791938424326609e-04 * solar_mass ), // Saturn ( 8.34336671824457987e00, 4.12479856412430479e00, -4.03523417114321381e-01, -2.76742510726862411e-03 * days_per_year, 4.99852801234917238e-03 * days_per_year, 2.30417297573763929e-05 * days_per_year, 2.85885980666130812e-04 * solar_mass ), // Uranus ( 1.28943695621391310e01, -1.51111514016986312e01, -2.23307578892655734e-01, 2.96460137564761618e-03 * days_per_year, 2.37847173959480950e-03 * days_per_year, -2.96589568540237556e-05 * days_per_year, 4.36624404335156298e-05 * solar_mass ), // Neptune ( 1.53796971148509165e01, -2.59193146099879641e01, 1.79258772950371181e-01, 2.68067772490389322e-03 * days_per_year, 1.62824170038242295e-03 * days_per_year, -9.51592254519715870e-05 * days_per_year, 5.15138902046611451e-05 * solar_mass ) ); Main(String* args) { if |args| != 1: Print("Usage: nbody \n"); return; ; res = _parse_(args(0)); if not res :: : Print("Usage: nbody \n"); return; ; iteractions = _untag_(res); positions = (vector(x: x, y: y, z: z) : x, y, z, vx, vy, vz, m <- bodies); velocities = (vector(x: vx, y: vy, z: vz) : x, y, z, vx, vy, vz, m <- bodies); masses = (m : x, y, z, vx, vy, vz, m <- bodies); velocities = offset_momentum(velocities, masses); e = energy(positions, velocities, masses); Print(_print_(e)); Print("\n"); positions, velocities = advance(iteractions, positions, velocities, masses, 0.01); e = energy(positions, velocities, masses); Print(_print_(e)); Print("\n"); } //////////////////////////////////////////////////////////////////////////////// type Int = <*..*>; type Float = ; type Nat = <0..*>; type String = string(Nat*); Int (_+_) (Int a, Int b) = _add_(a, b); Float (-_) (Float x) = _fneg_(x); Float (_+_) (Float x, Float y) = _fadd_(x, y); Float (_-_) (Float x, Float y) = _fsub_(x, y); Float (_*_) (Float x, Float y) = _fmult_(x, y); Float (_/_) (Float x, Float y) = _fdiv_(x, y); Float sqrt(Float x) = _fsqrt_(x);