#include <cstdlib>
#include <cmath>

#include <vector>
#include <tr1/array>
#include <iostream>
#include <map>

#include "physics.hh"

typedef std::tr1::array<double,3> vector3d;

static const int seed = 42;

/**
 * Generate initial condition for the three body problem
 */
void threebody (std::vector<vector3d>& r,
    std::vector<vector3d>& v,
    std::vector<double>& m)
{
    r.resize(3);
    v.resize(3);
    m.resize(3);

    m[0] = 3e12;
    m[1] = 1.0;
    m[2] = 1.0;

    r[0][0] = 0.0; r[0][1] = 0.0; r[0][2] = 0.0;
    r[1][0] = 1.0; r[1][1] = 0.0; r[1][2] = 0.0;
    r[2][0] =-1.0; r[2][1] = 0.0; r[2][2] = 0.0;

    v[0][0] = 0.0; v[0][1] = 0.0; v[0][2] = 0.0;
    v[1][0] = 0.0; v[1][1] = 0.1; v[1][2] = 0.0;
    v[2][0] = 0.0; v[2][1] =-0.1; v[2][2] = 0.0;
}

// hide helper function
namespace {
    /**
     * Compute orbital velocity for given body (helper function)
     */
    void orbitalVelocity (int p1, int p2,
        std::vector<vector3d>& x,
        std::vector<vector3d>& v,
        std::vector<double>& m)
    {
        static const double gamma_1 = gamma_si/(pc_in_m*pc_in_m*pc_in_m)*mass_sun*(365.25*86400)*(365.25*86400);
        double x1 = x[p1][0],
            y1 = x[p1][1],
            m1 = m[p1];
        double x2 = x[p2][0],
            y2 = x[p2][1];

        // calculate distance from the planet p1
        double r[2], dist;
        r[0] = x1 - x2;
        r[1] = y1 - y2;

        // distance in parsec
        dist = sqrt(r[0] * r[0] + r[1] * r[1]);

        // based on the distance from the sun, calculate the velocity needed to maintain a circular orbit
        double abs_v = sqrt(gamma_1 * m1 / dist);

        // calculate a suitable vector perpendicular to r for the velocity of the tracer
        v[p2][0] = ( r[1] / dist) * abs_v;
        v[p2][1] = (-r[0] / dist) * abs_v;
        v[p2][2] = 0.0;
    }
} // end namespace

/**
 * Generate initial condition for colliding galaxies
 */
void collision (std::vector<vector3d>& x,
    std::vector<vector3d>& v,
    std::vector<double>& m)
{
    srand48(seed);

    // initialize particles
    double ratio = 0.8;
    int b1 = 0;
    int b2 = (int) (ratio*m.size());
    int i;

    // initialize 1st black hole
    x[b1][0] = x[b1][1] = x[b1][2] = 0.0;
    v[b1][0] = v[b1][1] = v[b1][2] = 0.0;
    m[b1] = 1e6; // 1 million Sonnenmassen

    // initialize 1st galaxy
    for (i=1; i<b2; ++i)
    {
        const double rad = 10;
        double r = 0.1 + 0.8 * rad * drand48();
        double a = 2.0 * M_PI * drand48();
        m[i] = 0.03 + 20 * drand48();
        x[i][0] = r*sin(a);
        x[i][1] = r*cos(a);
        x[i][2] = 0.0;

        orbitalVelocity(b1,i,x,v,m);
    }

    // initialize 2nd black hole
    x[b2][0] = x[b2][1] = 10.0; x[b2][2] = 0.0;
    m[b2] = 1e5; // 10000 tausend Sonnenmassen

    orbitalVelocity(b1,b2,x,v,m);
    v[b2][0] *= 0.9; v[b2][1] *= 0.9;

    // initialize 2st galaxy
    for (i=b2+1; i<m.size(); ++i)
    {
        const double rad = 3;
        double r = 0.1 + 0.8 * rad * drand48();
        double a = 2.0 * M_PI * drand48();
        m[i] = 0.03 + 20 * drand48();
        x[i][0] = x[b2][0] + r*sin(a);
        x[i][1] = x[b2][1] + r*cos(a);
        x[i][2] = 0.0;

        orbitalVelocity(b2,i,x,v,m);
        v[i][0] += v[b2][0];
        v[i][1] += v[b2][1];
    }
}

/**
 * Choose initial condition
 */
void initialize (std::vector<vector3d>& r,
    std::vector<vector3d>& v,
    std::vector<double>& m,
    std::map<std::string, std::string> options)
{
    std::string init = options["init"];
    if (init == "threebody")
        threebody(r,v,m);
    // if (init == "plummer")
    //     plummer(r,v,m);
    if (init == "collision")
        collision(r,v,m);
}