#include <stdio.h> //standard libraries to be included
#include <math.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>

#define MAX_N_PLANETS 100 //maximum limits placed on a single simulation
#define MAX_NUMBER_EQUATIONS 6*MAX_N_PLANETS


#define GUESS 0.001 //guess for deltat
#define INTEGRATION_EPSILON 1.e-12
#define EQUATIONS_OF_MOTION(state, deriv) \
  { double fij, vdotr, factor1, factor2; \
    double dx, dy, dz, r2, r; \
    int c1, c2; \
    int i, j; \
    for(i=0; i<n_planets; i++){ \
      c1 = 6*i;	\
      (deriv).x[c1]   = (state).x[c1+3]; \
      (deriv).x[c1+1] = (state).x[c1+4]; \
      (deriv).x[c1+2] = (state).x[c1+5]; \
      (deriv).x[c1+3] = 0.0; \
      (deriv).x[c1+4] = 0.0; \
      (deriv).x[c1+5] = 0.0; \
    } \
    for(i=0; i<n_planets; i++){ c1 = 6*i; \
      for(j=0; j<i; j++){ \
        c2 = 6*j; \
	dx = (state).x[c1]   - (state).x[c2];   \
	dy = (state).x[c1+1] - (state).x[c2+1]; \
	dz = (state).x[c1+2] - (state).x[c2+2]; \
	r2 = dx*dx + dy*dy + dz*dz; \
	fij = 1.0/(r2*sqrt(r2)); \
	(deriv).x[c1+3] += - GM[j]*fij*dx; \
	(deriv).x[c1+4] += - GM[j]*fij*dy; \
	(deriv).x[c1+5] += - GM[j]*fij*dz; \
	(deriv).x[c2+3] +=   GM[i]*fij*dx; \
	(deriv).x[c2+4] +=   GM[i]*fij*dy; \
	(deriv).x[c2+5] +=   GM[i]*fij*dz; \
      } \
    } \
  }

#include "localmath.h" //local files that are used
#include "dynamics.h"

//Global variables
double G; //gravitational constant in AU^3*SM^-1*yr^-2
double GM[MAX_N_PLANETS]; //G * Mplanet for all objects including star
State state; //state.t contains time, state.x[] array contains 3D positions and velocities for everything
int n_planets; //number of planets and particles used in simulation--NOTE: n_planets includes the star
int number_equations; //total number of equations used in simulation for integrator.c
double machine_epsilon;

#include "integrator.c" //local files that are used and require the above to be declared prior to their inclusion
#include "machine-epsilon.c"
#include "kepcart.c"

double determine_machine_epsilon();

int nbody(int argc, void *argv[])
{
  //Inputs
  int j = 0;
  double curr_time = *(double *) argv[j]; j++; //the current time
  double desired_time = *(double *) argv[j]; j++; //the time to integrate forward to
  n_planets = *(int *) argv[j]; j++; //integer # of planets (incl. star) (global variable)
  double *Mplanet = (double *) argv[j]; j++; //array of planet masses  
  double *x = (double *) argv[j]; j++; //array of x coordinates
  double *y = (double *) argv[j]; j++; //array of y coordinates
  double *z = (double *) argv[j]; j++; //array of z coordinates
  double *vx = (double *) argv[j]; j++; //array of vx coordinates
  double *vy = (double *) argv[j]; j++; //array of vy coordinates
  double *vz = (double *) argv[j]; j++; //array of vz coordinates

  
  //Internal variables
  int i, c1;
  double deltat; //a guess at what the time sub-resolution should be
  double zero_deltat;


  //Initialize some global variables
  number_equations=6*n_planets;
  G = 4.0 * PI * PI; //gravitational constant in AU^3*SM^-1*yr^-2
  for(i=0;i<n_planets;i++) GM[i] = G*Mplanet[i];
  

  //Determine machine epsilon
  machine_epsilon = determine_machine_epsilon();

  //The Bulirsch-Stoer integrator expects a "state" structure
  //Copy the input coordinates into the global state structure
  state.t = curr_time;
  for(i=0;i<n_planets;i++) {
    c1 = 6*i;
    state.x[c1] = x[i];
    state.x[c1+1] = y[i];
    state.x[c1+2] = z[i];
    state.x[c1+3] = vx[i];
    state.x[c1+4] = vy[i];
    state.x[c1+5] = vz[i];
  }

  //Integrate the state of the system to the desired time
  deltat = GUESS;
  while(state.t < desired_time) {
    integrate(&state, &deltat);
  }
  zero_deltat = desired_time - state.t;
  while(ABS(zero_deltat)>INTEGRATION_EPSILON){
    integrate(&state, &zero_deltat);
    zero_deltat = desired_time - state.t;
  }

  //Now update the input arrays with the new state
  for(i=0;i<n_planets;i++) {
    c1 = 6*i;
    x[i] = state.x[c1];
    y[i] = state.x[c1+1];
    z[i] = state.x[c1+2];
    vx[i] = state.x[c1+3];
    vy[i] = state.x[c1+4];
    vz[i] = state.x[c1+5];
  }

  
  //Finished
  return 0;
}



