gistfile1.rs

import std.math.{PI, sqrt};

struct Body {
    x: f64;
    y: f64;
    z: f64;
    vx: f64;
    vy: f64;
    vz: f64;
    mass: f64;
};

const SOLAR_MASS: f64 = 4 * PI * PI;
const DAYS_PER_YEAR: f64 = 365.24;

fn(*Body) offset_momentum(px: f64, py: f64, pz: f64) {
    b.vx = -px / SOLAR_MASS;
    b.vy = -py / SOLAR_MASS;
    b.vz = -pz / SOLAR_MASS;
}

struct System(*[Body]);

fn new_system(bodies: *[Body]) -> System {
    let system = System(bodies);
    let mut (px, py, pz) = (0.0, 0.0, 0.0);
    for bodies.each |body| {
        px += body.vx * body.mass;
        py += body.vy * body.mass;
        pz += body.vz * body.mass;
    }

    bodies[0].offset_momentum(px, py, pz);
    system;
}

fn(*System) energy() -> f64 {
    let mut e = 0.0;
    for self.eachi |body_a, i| {
        e += 0.5 * body_a.mass * (body_a.vx * body_a.vx +
                                  body_a.vy * body_a.vy +
                                  body_a.vz * body_a.vz);

        for self.slice_from(i + 1).each |body_b| {
            let dx = body_a.x - body_b.x;
            let dy = body_a.y - body_b.y;
            let dz = body_a.z - body_b.z;
            let distance = sqrt(dx*dx + dy*dy + dz*dz);
            e -= (body_a.mass * body_b.mass) / distance;
        }
    }

    e;
}

fn(*System) advance(dt: f64) {
    for self.eachi |body_a, i| {
        for self.slice_from(i + 1).each |body_b| {
            let dx = body_a.x - body_b.x;
            let dy = body_a.y - body_b.y;
            let dz = body_a.z - body_b.z;

            let d_squared = dx*dx + dy*dy + dz*dz;
            let distance = sqrt(d_squared);
            let magnitude = dt / (d_squared * distance);

            body_a.vx -= dx * body_b.mass * magnitude;
            body_a.vy -= dy * body_b.mass * magnitude;
            body_a.vz -= dz * body_b.mass * magnitude;

            body_b.vx += dx * body_a.mass * magnitude;
            body_b.vy += dy * body_a.mass * magnitude;
            body_b.vz += dz * body_a.mass * magnitude;
        }
    }

    for self.each |body| {
        body.x += dt * body.vx;
        body.y += dt * body.vy;
        body.z += dt * body.vz;
    }
}

const BODIES: [Body/5] = [
   JUPITER = Body {
      x: 4.84143144246472090e+00,
      y: -1.16032004402742839e+00,
      z: -1.03622044471123109e-01,
      vx: 1.66007664274403694e-03 * DAYS_PER_YEAR,
      vy: 7.69901118419740425e-03 * DAYS_PER_YEAR,
      vz: -6.90460016972063023e-05 * DAYS_PER_YEAR,
      mass: 9.54791938424326609e-04 * SOLAR_MASS,
   },
   SATURN = Body {
      x: 8.34336671824457987e+00,
      y: 4.12479856412430479e+00,
      z: -4.03523417114321381e-01,
      vx: -2.76742510726862411e-03 * DAYS_PER_YEAR,
      vy: 4.99852801234917238e-03 * DAYS_PER_YEAR,
      vz: 2.30417297573763929e-05 * DAYS_PER_YEAR,
      mass: 2.85885980666130812e-04 * SOLAR_MASS,
   },
   URANUS = Body {
      x: 1.28943695621391310e+01,
      y: -1.51111514016986312e+01,
      z: -2.23307578892655734e-01,
      vx: 2.96460137564761618e-03 * DAYS_PER_YEAR,
      vy: 2.37847173959480950e-03 * DAYS_PER_YEAR,
      vz: -2.96589568540237556e-05 * DAYS_PER_YEAR,
      mass: 4.36624404335156298e-05 * SOLAR_MASS,
   },
   NEPTUNE = Body {
      x: 1.53796971148509165e+01,
      y: -2.59193146099879641e+01,
      z: 1.79258772950371181e-01,
      vx: 2.68067772490389322e-03 * DAYS_PER_YEAR,
      vy: 1.62824170038242295e-03 * DAYS_PER_YEAR,
      vz: -9.51592254519715870e-05 * DAYS_PER_YEAR,
      mass: 5.15138902046611451e-05 * SOLAR_MASS,
   },
   SUN = Body {
      x: 0.0,
      y: 0.0,
      z: 0.0,
      vx: 0.0,
      vy: 0.0,
      vz: 0.0,
      mass: SOLAR_MASS,
   }
];

fn main(args: [str]) {
    let n = int::from_str(args[0]);
    let system = new_system([ JUPITER, SATURN, URANUS, NEPTUNE, SUN ]);

    printf!("%.9f", system.energy());

    for n.times {
        system.advance(0.01);
    }

    printf!("%.9f", system.energy());
}