155 lines
3.9 KiB
Plaintext
155 lines
3.9 KiB
Plaintext
var math = use('math/radians')
|
|
var SOLAR_MASS = 4 * pi * pi;
|
|
var DAYS_PER_YEAR = 365.24;
|
|
|
|
function Body(x, y, z, vx, vy, vz, mass) {
|
|
return {x, y, z, vx, vy, vz, mass};
|
|
}
|
|
|
|
function Jupiter() {
|
|
return Body(
|
|
4.84143144246472090e+00,
|
|
-1.16032004402742839e+00,
|
|
-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
|
|
);
|
|
}
|
|
|
|
function Saturn() {
|
|
return Body(
|
|
8.34336671824457987e+00,
|
|
4.12479856412430479e+00,
|
|
-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
|
|
);
|
|
}
|
|
|
|
function Uranus() {
|
|
return Body(
|
|
1.28943695621391310e+01,
|
|
-1.51111514016986312e+01,
|
|
-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
|
|
);
|
|
}
|
|
|
|
function Neptune() {
|
|
return Body(
|
|
1.53796971148509165e+01,
|
|
-2.59193146099879641e+01,
|
|
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
|
|
);
|
|
}
|
|
|
|
function Sun() {
|
|
return Body(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, SOLAR_MASS);
|
|
}
|
|
|
|
var bodies = Array(Sun(), Jupiter(), Saturn(), Uranus(), Neptune());
|
|
|
|
function offsetMomentum() {
|
|
var px = 0;
|
|
var py = 0;
|
|
var pz = 0;
|
|
var size = length(bodies);
|
|
for (var i = 0; i < size; i++) {
|
|
var body = bodies[i];
|
|
var mass = body.mass;
|
|
px += body.vx * mass;
|
|
py += body.vy * mass;
|
|
pz += body.vz * mass;
|
|
}
|
|
|
|
var body = bodies[0];
|
|
body.vx = -px / SOLAR_MASS;
|
|
body.vy = -py / SOLAR_MASS;
|
|
body.vz = -pz / SOLAR_MASS;
|
|
}
|
|
|
|
function advance(dt) {
|
|
var size = length(bodies);
|
|
|
|
for (var i = 0; i < size; i++) {
|
|
var bodyi = bodies[i];
|
|
var vxi = bodyi.vx;
|
|
var vyi = bodyi.vy;
|
|
var vzi = bodyi.vz;
|
|
for (var j = i + 1; j < size; j++) {
|
|
var bodyj = bodies[j];
|
|
var dx = bodyi.x - bodyj.x;
|
|
var dy = bodyi.y - bodyj.y;
|
|
var dz = bodyi.z - bodyj.z;
|
|
|
|
var d2 = dx * dx + dy * dy + dz * dz;
|
|
var mag = dt / (d2 * math.sqrt(d2));
|
|
|
|
var massj = bodyj.mass;
|
|
vxi -= dx * massj * mag;
|
|
vyi -= dy * massj * mag;
|
|
vzi -= dz * massj * mag;
|
|
|
|
var massi = bodyi.mass;
|
|
bodyj.vx += dx * massi * mag;
|
|
bodyj.vy += dy * massi * mag;
|
|
bodyj.vz += dz * massi * mag;
|
|
}
|
|
bodyi.vx = vxi;
|
|
bodyi.vy = vyi;
|
|
bodyi.vz = vzi;
|
|
}
|
|
|
|
for (var i = 0; i < size; i++) {
|
|
var body = bodies[i];
|
|
body.x += dt * body.vx;
|
|
body.y += dt * body.vy;
|
|
body.z += dt * body.vz;
|
|
}
|
|
}
|
|
|
|
function energy() {
|
|
var e = 0;
|
|
var size = length(bodies);
|
|
|
|
for (var i = 0; i < size; i++) {
|
|
var bodyi = bodies[i];
|
|
|
|
e += 0.5 * bodyi.mass * ( bodyi.vx * bodyi.vx +
|
|
bodyi.vy * bodyi.vy + bodyi.vz * bodyi.vz );
|
|
|
|
for (var j = i + 1; j < size; j++) {
|
|
var bodyj = bodies[j];
|
|
var dx = bodyi.x - bodyj.x;
|
|
var dy = bodyi.y - bodyj.y;
|
|
var dz = bodyi.z - bodyj.z;
|
|
|
|
var distance = math.sqrt(dx * dx + dy * dy + dz * dz);
|
|
e -= (bodyi.mass * bodyj.mass) / distance;
|
|
}
|
|
}
|
|
return e;
|
|
}
|
|
|
|
var n = arg[0] || 100000
|
|
|
|
offsetMomentum();
|
|
|
|
log.console(`n = ${n}`)
|
|
log.console(energy().toFixed(9))
|
|
for (var i = 0; i < n; i++)
|
|
advance(0.01);
|
|
log.console(energy().toFixed(9))
|
|
|
|
$stop() |