Got the nbody demo working, and hopefully faster, too!
This commit is contained in:
Lindsey Kuper
Graydon Hoare
parent
1c1c67229a
commit
2650487120
1 changed files with 111 additions and 105 deletions
|
|
@ -4,8 +4,6 @@
|
|||
|
||||
fn main() {
|
||||
|
||||
log "foo";
|
||||
|
||||
let vec[int] inputs = vec(
|
||||
50000
|
||||
//these segfault :(
|
||||
|
|
@ -17,13 +15,12 @@ fn main() {
|
|||
let vec[Body.props] bodies = NBodySystem.MakeNBodySystem();
|
||||
|
||||
for (int n in inputs) {
|
||||
// TODO: make #fmt handle floats?
|
||||
log NBodySystem.energy(bodies);
|
||||
|
||||
let int i = 0;
|
||||
while (i < n) {
|
||||
bodies = NBodySystem.advance(bodies, 0.01);
|
||||
i = i+1;
|
||||
NBodySystem.advance(bodies, 0.01);
|
||||
i += 1;
|
||||
}
|
||||
log NBodySystem.energy(bodies);
|
||||
}
|
||||
|
|
@ -40,7 +37,6 @@ native "rust" mod rustrt {
|
|||
mod NBodySystem {
|
||||
|
||||
fn MakeNBodySystem() -> vec[Body.props] {
|
||||
// can't iterate over a record? how about a vector, then?
|
||||
let vec[Body.props] bodies = vec(
|
||||
// these each return a Body.props
|
||||
Body.sun(),
|
||||
|
|
@ -53,10 +49,13 @@ mod NBodySystem {
|
|||
let float py = 0.0;
|
||||
let float pz = 0.0;
|
||||
|
||||
for (Body.props body in bodies) {
|
||||
px += body.vx * body.mass;
|
||||
py += body.vy * body.mass;
|
||||
pz += body.vz * body.mass;
|
||||
let int i = 0;
|
||||
while (i < 5) {
|
||||
px += bodies.(i).vx * bodies.(i).mass;
|
||||
py += bodies.(i).vy * bodies.(i).mass;
|
||||
pz += bodies.(i).vz * bodies.(i).mass;
|
||||
|
||||
i += 1;
|
||||
}
|
||||
|
||||
// side-effecting
|
||||
|
|
@ -65,15 +64,15 @@ mod NBodySystem {
|
|||
ret bodies;
|
||||
}
|
||||
|
||||
fn advance(vec[Body.props] bodies, float dt) -> vec[Body.props] {
|
||||
for (Body.props ibody in bodies) {
|
||||
fn advance(vec[Body.props] bodies, float dt) -> () {
|
||||
|
||||
let Body.props iBody = ibody;
|
||||
|
||||
for (Body.props jbody in bodies) {
|
||||
let float dx = iBody.x - jbody.x;
|
||||
let float dy = iBody.y - jbody.y;
|
||||
let float dz = iBody.z - jbody.z;
|
||||
let int i = 0;
|
||||
while (i < 5) {
|
||||
let int j = i+1;
|
||||
while (j < 5) {
|
||||
let float dx = bodies.(i).x - bodies.(j).x;
|
||||
let float dy = bodies.(i).y - bodies.(j).y;
|
||||
let float dz = bodies.(i).z - bodies.(j).z;
|
||||
|
||||
let float dSquared = dx * dx + dy * dy + dz * dz;
|
||||
|
||||
|
|
@ -81,23 +80,30 @@ mod NBodySystem {
|
|||
rustrt.squareroot(dSquared, distance);
|
||||
let float mag = dt / (dSquared * distance);
|
||||
|
||||
iBody.vx -= dx * jbody.mass * mag;
|
||||
iBody.vy -= dy * jbody.mass * mag;
|
||||
iBody.vz -= dz * jbody.mass * mag;
|
||||
bodies.(i).vx -= dx * bodies.(j).mass * mag;
|
||||
bodies.(i).vy -= dy * bodies.(j).mass * mag;
|
||||
bodies.(i).vz -= dz * bodies.(j).mass * mag;
|
||||
|
||||
bodies.(j).vx += dx * bodies.(i).mass * mag;
|
||||
bodies.(j).vy += dy * bodies.(i).mass * mag;
|
||||
bodies.(j).vz += dz * bodies.(i).mass * mag;
|
||||
|
||||
j += 1;
|
||||
|
||||
jbody.vx += dx * iBody.mass * mag;
|
||||
jbody.vy += dy * iBody.mass * mag;
|
||||
jbody.vz += dz * iBody.mass * mag;
|
||||
}
|
||||
}
|
||||
|
||||
for (Body.props body in bodies) {
|
||||
body.x += dt * body.vx;
|
||||
body.y += dt * body.vy;
|
||||
body.z += dt * body.vz;
|
||||
i += 1;
|
||||
}
|
||||
|
||||
ret bodies;
|
||||
i = 0;
|
||||
while (i < 5) {
|
||||
|
||||
bodies.(i).x += dt * bodies.(i).vx;
|
||||
bodies.(i).y += dt * bodies.(i).vy;
|
||||
bodies.(i).z += dt * bodies.(i).vz;
|
||||
|
||||
i += 1;
|
||||
}
|
||||
}
|
||||
|
||||
fn energy(vec[Body.props] bodies) -> float {
|
||||
|
|
@ -107,110 +113,110 @@ mod NBodySystem {
|
|||
let float distance;
|
||||
let float e = 0.0;
|
||||
|
||||
for (Body.props ibody in bodies) {
|
||||
let int i = 0;
|
||||
while (i < 5) {
|
||||
e += 0.5 * bodies.(i).mass *
|
||||
( bodies.(i).vx * bodies.(i).vx
|
||||
+ bodies.(i).vy * bodies.(i).vy
|
||||
+ bodies.(i).vz * bodies.(i).vz );
|
||||
|
||||
// do we need this?
|
||||
let Body.props iBody = ibody;
|
||||
|
||||
e += 0.5 * iBody.mass *
|
||||
( iBody.vx * iBody.vx
|
||||
+ iBody.vy * iBody.vy
|
||||
+ iBody.vz * iBody.vz );
|
||||
|
||||
for (Body.props jbody in bodies) {
|
||||
|
||||
// do we need this?
|
||||
let Body.props jBody = jbody;
|
||||
|
||||
dx = iBody.x - jBody.x;
|
||||
dy = iBody.y - jBody.y;
|
||||
dz = iBody.z - jBody.z;
|
||||
let int j = i+1;
|
||||
while (j < 5) {
|
||||
dx = bodies.(i).x - bodies.(j).x;
|
||||
dy = bodies.(i).y - bodies.(j).y;
|
||||
dz = bodies.(i).z - bodies.(j).z;
|
||||
|
||||
rustrt.squareroot(dx*dx + dy*dy + dz*dz, distance);
|
||||
e -= (iBody.mass * jBody.mass) / distance;
|
||||
e -= (bodies.(i).mass * bodies.(j).mass) / distance;
|
||||
|
||||
j += 1;
|
||||
}
|
||||
|
||||
i += 1;
|
||||
}
|
||||
ret e;
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
mod Body {
|
||||
|
||||
const float PI = 3.14;
|
||||
const float SOLAR_MASS = 39.47; // was 4 * PI * PI originally
|
||||
|
||||
const float PI = 3.141592653589793;
|
||||
const float SOLAR_MASS = 39.478417604357432; // was 4 * PI * PI originally
|
||||
const float DAYS_PER_YEAR = 365.24;
|
||||
|
||||
type props = rec(float x,
|
||||
float y,
|
||||
float z,
|
||||
float vx,
|
||||
float vy,
|
||||
float vz,
|
||||
type props = rec(mutable float x,
|
||||
mutable float y,
|
||||
mutable float z,
|
||||
mutable float vx,
|
||||
mutable float vy,
|
||||
mutable float vz,
|
||||
float mass);
|
||||
|
||||
fn jupiter() -> Body.props {
|
||||
// current limitation of the float lexer: decimal part has to
|
||||
// fit into a 32-bit int.
|
||||
|
||||
let Body.props p;
|
||||
p.x = 4.84e+00;
|
||||
p.y = -1.16e+00;
|
||||
p.z = -1.03e-01;
|
||||
p.vx = 1.66e-03 * DAYS_PER_YEAR;
|
||||
p.vy = 7.69e-03 * DAYS_PER_YEAR;
|
||||
p.vz = -6.90e-05 * DAYS_PER_YEAR;
|
||||
p.mass = 9.54e-04 * SOLAR_MASS;
|
||||
ret p;
|
||||
ret rec(
|
||||
mutable x = 4.84143144246472090e+00,
|
||||
mutable y = -1.16032004402742839e+00,
|
||||
mutable z = -1.03622044471123109e-01,
|
||||
mutable vx = 1.66007664274403694e-03 * DAYS_PER_YEAR,
|
||||
mutable vy = 7.69901118419740425e-03 * DAYS_PER_YEAR,
|
||||
mutable vz = -6.90460016972063023e-05 * DAYS_PER_YEAR,
|
||||
mass = 9.54791938424326609e-04 * SOLAR_MASS
|
||||
);
|
||||
}
|
||||
|
||||
fn saturn() -> Body.props {
|
||||
let Body.props p;
|
||||
p.x = 8.34e+00;
|
||||
p.y = 4.12e+00;
|
||||
p.z = -4.03e-01;
|
||||
p.vx = -2.76e-03 * DAYS_PER_YEAR;
|
||||
p.vy = 4.99e-03 * DAYS_PER_YEAR;
|
||||
p.vz = 2.30e-05 * DAYS_PER_YEAR;
|
||||
p.mass = 2.85e-04 * SOLAR_MASS;
|
||||
ret p;
|
||||
ret rec(
|
||||
mutable x = 8.34336671824457987e+00,
|
||||
mutable y = 4.12479856412430479e+00,
|
||||
mutable z = -4.03523417114321381e-01,
|
||||
mutable vx = -2.76742510726862411e-03 * DAYS_PER_YEAR,
|
||||
mutable vy = 4.99852801234917238e-03 * DAYS_PER_YEAR,
|
||||
mutable vz = 2.30417297573763929e-05 * DAYS_PER_YEAR,
|
||||
mass = 2.85885980666130812e-04 * SOLAR_MASS
|
||||
);
|
||||
}
|
||||
|
||||
fn uranus() -> Body.props {
|
||||
let Body.props p;
|
||||
p.x = 1.28e+01;
|
||||
p.y = -1.51e+01;
|
||||
p.z = -2.23e-01;
|
||||
p.vx = 2.96e-03 * DAYS_PER_YEAR;
|
||||
p.vy = 2.37e-03 * DAYS_PER_YEAR;
|
||||
p.vz = -2.96e-05 * DAYS_PER_YEAR;
|
||||
p.mass = 4.36e-05 * SOLAR_MASS;
|
||||
ret p;
|
||||
ret rec(
|
||||
mutable x = 1.28943695621391310e+01,
|
||||
mutable y = -1.51111514016986312e+01,
|
||||
mutable z = -2.23307578892655734e-01,
|
||||
mutable vx = 2.96460137564761618e-03 * DAYS_PER_YEAR,
|
||||
mutable vy = 2.37847173959480950e-03 * DAYS_PER_YEAR,
|
||||
mutable vz = -2.96589568540237556e-05 * DAYS_PER_YEAR,
|
||||
mass = 4.36624404335156298e-05 * SOLAR_MASS
|
||||
);
|
||||
}
|
||||
|
||||
fn neptune() -> Body.props {
|
||||
let Body.props p;
|
||||
p.x = 1.53e+01;
|
||||
p.y = -2.59e+01;
|
||||
p.z = 1.79e-01;
|
||||
p.vx = 2.68e-03 * DAYS_PER_YEAR;
|
||||
p.vy = 1.62e-03 * DAYS_PER_YEAR;
|
||||
p.vz = -9.51e-05 * DAYS_PER_YEAR;
|
||||
p.mass = 5.15e-05 * SOLAR_MASS;
|
||||
ret p;
|
||||
ret rec(
|
||||
mutable x = 1.53796971148509165e+01,
|
||||
mutable y = -2.59193146099879641e+01,
|
||||
mutable z = 1.79258772950371181e-01,
|
||||
mutable vx = 2.68067772490389322e-03 * DAYS_PER_YEAR,
|
||||
mutable vy = 1.62824170038242295e-03 * DAYS_PER_YEAR,
|
||||
mutable vz = -9.51592254519715870e-05 * DAYS_PER_YEAR,
|
||||
mass = 5.15138902046611451e-05 * SOLAR_MASS
|
||||
);
|
||||
}
|
||||
|
||||
fn sun() -> Body.props {
|
||||
let Body.props p;
|
||||
p.mass = SOLAR_MASS;
|
||||
ret p;
|
||||
ret rec(
|
||||
mutable x = 0.0,
|
||||
mutable y = 0.0,
|
||||
mutable z = 0.0,
|
||||
mutable vx = 0.0,
|
||||
mutable vy = 0.0,
|
||||
mutable vz = 0.0,
|
||||
mass = SOLAR_MASS
|
||||
);
|
||||
}
|
||||
|
||||
impure fn offsetMomentum(mutable Body.props props,
|
||||
float px,
|
||||
float py,
|
||||
float pz) -> Body.props {
|
||||
|
||||
impure fn offsetMomentum(&Body.props props,
|
||||
float px,
|
||||
float py,
|
||||
float pz) -> () {
|
||||
props.vx = -px / SOLAR_MASS;
|
||||
props.vy = -py / SOLAR_MASS;
|
||||
props.vz = -pz / SOLAR_MASS;
|
||||
|
|
|
|||
Loading…
Reference in a new issue