#include <stdio.h>
#include <stdlib.h>
#include <math.h>

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

/*
 * Global constants
 */

/* set gravity constant to one */
const double G = 1.0;
/* The epsilon^2 for the softening */
const double epsilon2 = 1E-14;

/**
 * Compute acceleration
 *
 * \param[in]  n Number of bodies
 * \param[in]  r Positions of the bodies
 * \param[in]  m Masses of the bodies
 * \param[out] a The computed accelerations of the bodies
 */
void acceleration (int n, vector3 r[], double m[], vector3 a[])
{
}

/**
 * Do one time step
 *
 * \param[in]     n    Number of bodies
 * \param[in]     dt   Time step size
 * \param[in,out] x    IN: old positions, OUT: new positions
 * \param[in,out] v    IN: old velocities, OUT: new velocities
 * \param[in]     m    Masses
 * \param[in,out] a    IN: old accelerations, OUT: new accelerations
 * \param         aold Buffer to hold the old accelerations.  Input and output
 *                     contents is not meaningful.
 */
void euler (const int N, double dt, vector3 r[], vector3 v[], double m[], vector3 a[])
{
}


/**
 * main program
 */
int main ()
{
    /* number of bodies */
    const int N = 50;
    /* number of time steps to do */
    const int timesteps = 1000000;
    /* delta t */
    const double dt = 1E-5;
    /* when to measure MFLOP rate and write output files: every mod time steps */
    const int mod = 10000;

    /* arrays for position, velocity, mass, and acceleration of the bodies */
    vector3 r[N];
    vector3 v[N];
    double  m[N];
    vector3 a[N];
    
    /* loop counter */
    int i;
    /* count the time steps done */
    int k = 0;
    /* accumulated simulation time */
    double t = 0;

    /* get initial values from one of the generator functions */
    plummer(N,r,v,m);
    
    /* initialise the acceleration */
    for (i=0; i<N; i++)
        a[i][0] = a[i][1] = a[i][2] = 0.0;
    acceleration(N,r,m,a);

    /* write out the initial values */
    nbody_write(N,r,v,m,t,dt, k/mod);

    /* do the simulation */
    for (k=1; k<timesteps; k++) {
        /* do one timestep and increment elapsed simulated time */
        euler(N,dt,r,v,m,a);
        t += dt;

        /* write solution if necessary */
        if (k%mod==0) {
            /* write output file */
            nbody_write(N,r,v,m,t,dt, k/mod);
        }

        /* compute new acceleration */
        acceleration(N,r,m,a);
    }

    /* we are done */
    return 0;
}