/* 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[3] p, v;
            imd mass;

            void offsetMomentum(in double[3] pn) pure {
                v[] = -pn[] / 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 {
            double[3] pn = 0.0;

            foreach (const(Body) b; bodies)
                pn[] += b.v[] * b.mass;

            bodies[0].offsetMomentum(pn);
        }


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

            foreach (i, b1; bodies) {
                imd im = b1.mass;
                // b1.v[] ^^ 2 not supported yet
                e += 0.5 * im * (b1.v[0] ^^ 2 + b1.v[1] ^^ 2 + b1.v[2] ^^ 2);

                foreach (b2; bodies[i + 1 .. $]) {
                    // immutable d = b1.p[] - b2.p[]; // can't be immutable yet
                    const double[3] d = b1.p[] - b2.p[];
                    e -= (im * b2.mass) / sqrt(d[0] ^^ 2 + d[1] ^^ 2 + d[2] ^^ 2);
                }
            }

            return e;
        }


        void advance(double dt)() pure {
            enum int NPAIR = NB * (NB - 1) / 2;

            double[3][NPAIR] r = void;
            int k;
            foreach (i; Iota!(0, NB))
                foreach (j; Iota!(i + 1, NB)) {
                    r[k][] = bodies[i].p[] - bodies[j].p[];
                    k++;
                }

            //immutable double[NPAIR] distance = sqrt(sum(r ^^ 2, dim=1));
            double[NPAIR] distance = void;
            foreach (i; Iota!(0, NPAIR))
                distance[i] = sqrt(r[i][0] ^^ 2 + r[i][1] ^^ 2 + r[i][2] ^^ 2);

            //immutable double[NPAIR] mag = dt / distance[] ^^ 3;
            double[NPAIR] mag = void;
            foreach (i; Iota!(0, NPAIR))
                mag[i] = dt / distance[i] ^^ 3;

            k = 0;
            foreach (i; Iota!(0, NB))
                foreach (j; Iota!(i + 1, NB)) {
                    const double[3] rmag = r[k][] * mag[k];
                    bodies[i].v[] -= bodies[j].mass * rmag[];
                    bodies[j].v[] += bodies[i].mass * rmag[];
                    k++;
                }

            foreach (i; Iota!(0, NB))
                bodies[i].p[] += dt * bodies[i].v[];
        }
}


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