Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- import std.stdio; [59/187]
- // import std.math;
- import core.stdc.stdio;
- import core.stdc.math;
- enum pi = 3.141592653589793;
- enum solar_mass = (4 * pi * pi);
- enum days_per_year = 365.24;
- struct Planet {
- double x, y, z;
- double vx, vy, vz;
- double mass;
- }
- void advance(Planet[] bodies, double dt)
- {
- foreach (i; 0 .. bodies.length) {
- auto b = &(bodies[i]);
- foreach (j; i+1 .. bodies.length) {
- auto b2 = &(bodies[j]);
- auto dx = b.x - b2.x;
- auto dy = b.y - b2.y;
- auto dz = b.z - b2.z;
- auto distance = sqrt(dx * dx + dy * dy + dz * dz);
- auto mag = dt / (distance * distance * distance);
- b.vx -= dx * b2.mass * mag;
- b.vy -= dy * b2.mass * mag;
- b.vz -= dz * b2.mass * mag;
- b2.vx += dx * b.mass * mag;
- b2.vy += dy * b.mass * mag;
- b2.vz += dz * b.mass * mag;
- }
- }
- foreach (i; 0 .. bodies.length) {
- auto b = &(bodies[i]);
- b.x += dt * b.vx;
- b.y += dt * b.vy;
- b.z += dt * b.vz;
- }
- }
- double energy(Planet[] bodies)
- {
- double e = 0.0;
- foreach (i; 0 .. bodies.length) {
- auto b = &(bodies[i]);
- e += 0.5 * b.mass * (b.vx * b.vx + b.vy * b.vy + b.vz * b.vz);
- foreach (j; i+1 .. bodies.length) {
- auto b2 = &(bodies[j]);
- auto dx = b.x - b2.x;
- auto dy = b.y - b2.y;
- auto dz = b.z - b2.z;
- auto distance = sqrt(dx * dx + dy * dy + dz * dz);
- e -= (b.mass * b2.mass) / distance;
- }
- }
- return e;
- }
- void offset_momentum(Planet[] bodies)
- {
- double px = 0.0, py = 0.0, pz = 0.0;
- foreach (i; 0 .. bodies.length) {
- px += bodies[i].vx * bodies[i].mass;
- py += bodies[i].vy * bodies[i].mass;
- pz += bodies[i].vz * bodies[i].mass;
- }
- bodies[0].vx = - px / solar_mass;
- bodies[0].vy = - py / solar_mass;
- bodies[0].vz = - pz / solar_mass;
- }
- enum NBODIES = 5;
- Planet[NBODIES] bodies = [
- { /* sun */
- x: 0.0, y: 0.0, z: 0.0, vx: 0.0, vy: 0.0, vz: 0.0, mass:solar_mass
- },
- { /* jupiter */
- 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
- },
- { /* saturn */
- 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
- },
- { /* uranus */
- 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
- },
- { /* neptune */
- 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
- }
- ];
- void main(string[] args)
- {
- import std.conv;
- int n = to!(int)(args[1]);
- offset_momentum(bodies);
- printf("%.9f\n", energy(bodies));
- foreach (i; 0 .. n)
- advance(bodies, 0.01);
- printf("%.9f\n", energy(bodies));
- }
Add Comment
Please, Sign In to add comment