Skip to content

Instantly share code, notes, and snippets.

@dherman
Forked from pcwalton/gist:2927607
Created June 14, 2012 05:30
Show Gist options
  • Save dherman/2928167 to your computer and use it in GitHub Desktop.
Save dherman/2928167 to your computer and use it in GitHub Desktop.
Rust of Patrick's dreams, with Dave's method syntax :)
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());
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment