/* The Great Computer Language Shootout
   http://shootout.alioth.debian.org/
   compile: dmd -O -inline -release nbody.d
*/

import core.stdc.stdio, std.conv, std.math, std.typetuple;

template Iota(int start, int stop) {
    static if (stop <= start)
        alias TypeTuple!() Iota;
    else
        alias TypeTuple!(Iota!(start, stop-1), stop-1) Iota;
}


struct BodySystem {
    alias immutable(double) imd;
    private:
        enum double SOLAR_MASS = 4 * PI ^^ 2;
        enum double DAYS_PER_YEAR = 365.24;

        static struct Body {
            double x, y, z, vx, vy, vz;
            imd mass;

            void offsetMomentum(in double px, in double py, in double pz) pure nothrow {
                vx = -px / SOLAR_MASS;
                vy = -py / SOLAR_MASS;
                vz = -pz / SOLAR_MASS;
            }
        }

        auto bodies = TypeTuple!(
            Body(0,0,0, 0,0,0, SOLAR_MASS),

            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),

            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),

            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),

            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));

            enum int NB = bodies.length;

    public:
        void offset() pure nothrow {
            double px = 0.0, py = 0.0, pz = 0.0;

            foreach (const(Body) b; bodies) {
                px += b.vx * b.mass;
                py += b.vy * b.mass;
                pz += b.vz * b.mass;
            }

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


        double energy() const pure nothrow {
            double e = 0.0;

            foreach (i, b1; bodies) {
                imd im = b1.mass;
                e += 0.5 * im * (b1.vx ^^ 2 + b1.vy ^^ 2 + b1.vz ^^ 2);

                foreach (b2; bodies[i + 1 .. $]) {
                    imd dx = b1.x - b2.x;
                    imd dy = b1.y - b2.y;
                    imd dz = b1.z - b2.z;

                    e -= (im * b2.mass) / sqrt(dx ^^ 2 + dy ^^ 2 + dz ^^ 2);
                }
            }

            return e;
        }


        void advance(double dt)() pure nothrow {
            foreach (i; Iota!(0, NB))
                foreach (j; Iota!(i + 1, NB)) {
                    enum double im = bodies[i].mass;
                    enum double jm = bodies[j].mass;

                    imd dx = bodies[i].x - bodies[j].x;
                    imd dy = bodies[i].y - bodies[j].y;
                    imd dz = bodies[i].z - bodies[j].z;

                    imd dSquared = dx ^^ 2 + dy ^^ 2 + dz ^^ 2;
                    imd mag = dt / (dSquared * sqrt(dSquared));

                    bodies[i].vx -= dx * jm * mag;
                    bodies[i].vy -= dy * jm * mag;
                    bodies[i].vz -= dz * jm * mag;

                    bodies[j].vx += dx * im * mag;
                    bodies[j].vy += dy * im * mag;
                    bodies[j].vz += dz * im * mag;
                }

            foreach (i; Iota!(0, NB))
                with (bodies[i]) {
                    x += dt * vx;
                    y += dt * vy;
                    z += dt * vz;
                }
        }
}


void main(string[] args) {
    immutable int n = args.length > 1 ? to!int(args[1]) : 100;

    BodySystem bodies;
    bodies.offset();
    printf("%0.9f\n", bodies.energy());

    foreach (i; 0 .. n)
        bodies.advance!(0.01)();

    printf("%0.9f\n", bodies.energy());
}