#ifndef NBGENERATOR_HH
#define NBGENERATOR_HH

#include <cstdio>
#include <cstdlib>
#include <cmath>

#include "nbio.hh"
#include "nbtypes.hh"

/*
  generate initial condition for the nbody problem
*/
void plummer (int n, vector3 x[], vector3 v[], double m[],
    long int seed = 0)
{
    /* This is a copy from 
       Pit Hut, Jun Makino: The Art of Computational Science, 
       The Kali Code, vol. 5. Initial Conditions: Plummer's Model.
    */ 
    int i;
    double radius,theta,phi,X,Y,velocity,maxr=-1.0;
    const double Pi = 3.141592653589793238462643383279;
    double s[3] = {0.0,0.0,0.0};
    double t[3] = {0.0,0.0,0.0};
    if (seed!=0) srand48(seed);
    for (i=0; i<n; i++)
    {
        m[i] = 1.0/n;
        radius = 1.0/sqrt(pow(drand48(),-2.0/3.0)-1.0);
        if (radius>maxr) maxr=radius;
        theta = acos(-1.0+drand48()*2.0);
        phi = drand48()*2.0*Pi;
        x[i][0] = radius*sin(theta)*cos(phi);
        x[i][1] = radius*sin(theta)*sin(phi);
        x[i][2] = radius*cos(theta);
        s[0] += m[i]*x[i][0];
        s[1] += m[i]*x[i][1];
        s[2] += m[i]*x[i][2];
        X = 0.0;
        Y = 0.1;
        while (Y>X*X*pow(1.0-X*X,3.5))
        {
            X = drand48();
            Y = drand48()*0.1;
        }
        velocity = X*sqrt(2.0)*pow(1.0+radius*radius,-0.25);
        theta = acos(-1.0+drand48()*2.0);
        phi = drand48()*2.0*Pi;
        v[i][0] = velocity*sin(theta)*cos(phi);
        v[i][1] = velocity*sin(theta)*sin(phi);
        v[i][2] = velocity*cos(theta);
        t[0] += m[i]*v[i][0];
        t[1] += m[i]*v[i][1];
        t[2] += m[i]*v[i][2];
    }
    printf("center of mass: %g %g %g\n",s[0],s[1],s[2]);
    for (i=0; i<n; i++)
    {
        x[i][0] -= s[0]; 
        x[i][1] -= s[1]; 
        x[i][2] -= s[2]; 

        v[i][0] -= t[0]; 
        v[i][1] -= t[1]; 
        v[i][2] -= t[2]; 
    }
    s[0]=s[1]=s[2]=0.0;
    for (i=0; i<n; i++)
    {
        s[0] += m[i]*x[i][0];
        s[1] += m[i]*x[i][1];
        s[2] += m[i]*x[i][2];
    }
    printf("new center of mass: %g %g %g\n",s[0],s[1],s[2]);
    printf("maximum radius: %g\n",maxr);
}

void orbitalVelocity(int p1, int p2,
	      double x[][3], double v[][3], double m[])
{
  static const double mass_sun = 1.988435e30;
  static const double pc_in_m = 3.08567758129e16;
  static const double gamma_si = 6.67428e-11;
  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;
}


void collision(int m_num, long int seed,
	      double x[][3], double v[][3], double m[])
{
  if (seed!=0) srand(seed);

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

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

  // initialize 1st galaxy
  for (i=1; i<b2; ++i)
  {
    const double rad = 10;
    double r = 0.1 + .8 * (rad * ((double)rand() / RAND_MAX));
    double a = 2.0*M_PI*((double)rand() / RAND_MAX);
    m[i] = 0.03 + 20 * ((double)rand() / RAND_MAX);
    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] = 100000; // 100 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_num; ++i)
  {
    const double rad = 3;
    double r = 0.1 + .8 * (rad * ((double)rand() / RAND_MAX));
    double a = 2.0*M_PI*((double)rand() / RAND_MAX);
    m[i] = 0.03 + 20 * ((double)rand() / RAND_MAX);
    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];
  }
}

#endif // NBGENERATOR_HH