#include "gravity_support.h"
#include "GraphicsForGravity.h"
#include <stdio.h>
#include <math.h>
#include <string.h>


#define bool	int
#define TRUE	1
#define FALSE	0

#define G		6.673E-11   /* gravitational constant in SI units */

#define InitialTimeStep		1.0
/*
note: this is meant to be "very small" for the problem in question.
(Just in case by bad luck we're at the *start* (far point!) of a "comet-like"
 orbit...in that case all that cleverness in TIME_STEP_UP_BY comes to naught!)
*/

#define TIME_STEP_DOWN_BY	1.2
/*
for turning down the time step if there was an unsatisfactory update.
There's no point making it any smaller as it's MAX_REL_DISP that ultimately
controls accuracy. Any bigger and we waste time being "over-accurate" at
greater computational cost than just choosing a smaller MAX_REL_DISP.
*/

#define TIME_STEP_UP_BY		(1.0 + MAX_REL_DISP * (TIME_STEP_DOWN_BY - 1.0))
/*
for turning up the time step if there was a satisfactory update.
This choice ensures the time step never gets the chance to go spuriously high
at or around the far point of "comet-like" orbits. With a naive choice of
TIME_STEP_UP_BY, say just fixed at 1.01, the judging of safe time steps
by (temporarily very slow) changes in *position* is lulled into a false
sense of security. This definition avoids that problem. It also means at least
about a radian-worth of movement takes place before the time step
judders down; any finer-scale juddering would just be a waste.
*/


static double MAX_REL_DISP;	/* criterion for satisfactory update */
static int    PLOT_HOW_OFTEN;	/* one in this many states are plotted */
static int    PRINT_HOW_OFTEN;	/* one in this many states are printed */

static int state_count = 0;	/* out here for convenience */


/* private forward references */
static void dot (STATE *s, STATE_DOT *s_dot);
static bool WasSatisfactoryUpdate (STATE *CurrentState,
                                   double TimeStep,
                                   STATE *NextState);
static void ReadState (FILE *filehandle, STATE *AState);
static void PutState (STATE *AState, double ATimeStep);
static void PutEnsembleStatistics (STATE *AState);





void main(int argc, char *argv[])
{
  STATE  states[2];			/* TheState and NewState point here */
  STATE  *TheState, *NewState, *temp;	/* temp is just for swapping */
  double TimeStep;
  char   *filename;
  FILE   *filehandle;

  /*
   *	argument checking
   */
  if (argc != 2) {
    fprintf(stderr, "usage: %s <data-file-name>\n", argv[0]);
    exit(1);
  }/*endif*/

  filename=argv[1];

  /*
   *	do we have a readable file?
   */
  filehandle=fopen(filename, "r");
  if (filehandle == NULL) {
    fprintf(stderr, "%s: cannot open file \"%s\" for reading!\n", argv[0], filename);
    exit(1);
  }/*endif*/

  /*
   *	OK, we have a readable file - let's get the data!
   */
  if (fscanf(filehandle,
             "%lf %d %d%*[^\n]\n",
             &MAX_REL_DISP, &PLOT_HOW_OFTEN, &PRINT_HOW_OFTEN) != 3) {
    fprintf(stderr, "%s: bad input data\n", argv[0]);
    exit(1);
  }/*endif*/

  TheState = &states[0];
  NewState = &states[1];

  ReadState(filehandle, TheState);
  *NewState = *TheState;		/* to save having to update names etc */
  TimeStep = InitialTimeStep;

  InitializeGraphics(TheState);


  while (TRUE) {
    if (state_count % PLOT_HOW_OFTEN == 0) {
      DrawState(TheState);
    }/*endif*/;

    if (state_count % PRINT_HOW_OFTEN == 0) {
      PutState(TheState, TimeStep);
    }/*endif*/;

    while (! WasSatisfactoryUpdate(TheState, TimeStep, NewState)) {
      TimeStep /= TIME_STEP_DOWN_BY;
    }/*endwhile*/;
    temp = TheState; TheState = NewState; NewState = temp;
    TimeStep *= TIME_STEP_UP_BY;
    state_count++;
  }/*endwhile*/
}
             




/* the procedures at the core of the simulation... */

static void dot (STATE *s, STATE_DOT *s_dot)
/* finds the time derivative of the state. Treats first arg as read-only. */
{
  double mass_i, mass_j, multiplier;
  int    i,j,d;
  VECTOR i_to_j;

  s_dot->n = s->n;

  for (i=0; i<s->n; i++) {
    for (d=0; d<NDIMS; d++) {
      s_dot->bodydata[i].position_dot[d] = s->bodydata[i].velocity[d];
    }/*endfor*/
  }/*endfor*/

  for (i=0; i<s->n; i++) {
    CreateZeroVector (s_dot->bodydata[i].velocity_dot);
    mass_i = s->bodydata[i].mass;
    for (j=0; j<i; j++) {
      mass_j = s->bodydata[j].mass;

      if ((mass_i != 0.0) || (mass_j != 0.0)) {
        for (d=0; d<NDIMS; d++) {
          i_to_j[d] = s->bodydata[j].position[d] - s->bodydata[i].position[d];
        }/*endfor*/;
        multiplier = G / pow(lengthsquared(i_to_j), NDIMS/2.0);

        /* see if body i pulls on body j */
        if (mass_i != 0.0) {
          for (d=0; d<NDIMS; d++) {
            s_dot->bodydata[j].velocity_dot[d] -=
                    multiplier * mass_i * i_to_j[d];
          }/*endfor*/
        }/*endif*/

        /* see if body i gets pulled by body j */
        if (mass_j != 0.0) {
          for (d=0; d<NDIMS; d++) {
            s_dot->bodydata[i].velocity_dot[d] +=
                    multiplier * mass_j * i_to_j[d];
          }/*endfor*/
        }/*endif*/

      }/*endif*/
    }/*endfor*/
  }/*endfor*/
}





static bool WasSatisfactoryUpdate (STATE *CurrentState,
                                   double TimeStep,
                                   STATE *NextState)
/*
attempts update, and reports whether satisfactory.
treats first arg as read-only; assumes "static" parts of *NextState
are already filled in.


We use the Runge-Kutta method, as follows.

We define these intermediate estimates of the mean "speed" (change against time)
that the state changes at, between now and now + TimeStep:

  speed0 = dot(CurrentState)

  speed1 = dot(CurrentState + FRAC00 * TimeStep * speed0)

  speed2 = dot(CurrentState + FRAC10 * TimeStep * speed0
                            + FRAC11 * TimeStep * speed1)

  speed3 = dot(CurrentState + FRAC20 * TimeStep * speed0
                            + FRAC21 * TimeStep * speed1
                            + FRAC22 * TimeStep * speed2)

And so on, say up to speed_n-1.
We then weight these into a final estimate of the mean speed:

  speed  =   FRAC_n-1_0   * speed0
           + FRAC_n-1_1   * speed1
           + FRAC_n-1_2   * speed2
             .
             .
             .
           + FRAC_n-1_n-1 * speed_n-1

And finally, we say NextState = CurrentState + TimeStep * speed.

This agrees with the Taylor series expansion for the state up to some
particular high order in TimeStep (this order is not necessarily the "n" of
speed_n-1, etc, above, which we call the complexity),
leaving errors of the next order higher still
(or *relative* errors [error/TimeStep] of the same particular order).
However, it takes far less computation than any explicit evaluation
of the Taylor series would.

The magic numbers FRACij are given in the static local data structure below.
*/

{
  #define RK_IS_BUTCHER_ORDER_6


  #if defined(RK_IS_CLASSICAL_ORDER_4)
    #define RungeKuttaComplexity 4
    static double FRAC[RungeKuttaComplexity][RungeKuttaComplexity] = {
      { 0.5 },
      { 0,     0.5 },
      { 0,     0,     1 },
      { 1.0/6, 2.0/6, 2.0/6, 1.0/6 }
    };
  #elif defined(RK_IS_BUTCHER_ORDER_5)
    #define RungeKuttaComplexity 6
    static double FRAC[RungeKuttaComplexity][RungeKuttaComplexity] = {
      { 0.25 },
      { 0.125,   0.125 },
      { 0,      -0.5,    1 },
      { 3.0/16,  0,      0,        9.0/16 },
      {-3.0/7,   2.0/7,  12.0/7,  -12.0/7,  8.0/7 },
      { 7.0/90,  0,      32.0/90,  12.0/90, 32.0/90, 7.0/90 }
    };
  #elif defined(RK_IS_BUTCHER_ORDER_6)
    #define RungeKuttaComplexity 7
    static double FRAC[RungeKuttaComplexity][RungeKuttaComplexity] = {
      { 1.0/3 },
      { 0,        2.0/3 },
      { 1.0/12,   1.0/3,  -1.0/12 },
      {-1.0/16,   1.125,  -3.0/16,  -0.375 },
      { 0,        1.125,  -0.375,   -0.75,     0.5 },
      { 9.0/44,  -9.0/11,  63.0/44,  18.0/11,  0,      -16.0/11 },
      { 11.0/120, 0,       27.0/40,  27.0/40, -4.0/15, -4.0/15, 11.0/120 }
    };
  #endif


  STATE_DOT speed[RungeKuttaComplexity];
  VECTOR    current_i_to_j, new_i_to_j, change_in_i_to_j;
  int       i,j,d;


  /* compute NextState using Runge-Kutta as explained above */

  for (i=0; i<=RungeKuttaComplexity; i++) {
    copy_state(CurrentState, NextState);
    for (j=0; j<i; j++) {
      incr_state(NextState, FRAC[i-1][j] * TimeStep, &speed[j]);
    }/*endfor*/

    if (i<RungeKuttaComplexity) {
      dot(NextState, &speed[i]);
    }/*endif*/
  }/*endfor*/


  /* OK, that's NextState computed...now see if it's satisfactory! */

  for (i=0; i<CurrentState->n; i++) {
    for (j=0; j<i; j++) {
      if (CurrentState->bodydata[i].mass != 0 ||
          CurrentState->bodydata[j].mass != 0) {

        /* check non-trivial pair of bodies aren't swinging by too fast */
        for (d=0; d<NDIMS; d++) {
          current_i_to_j[d] = CurrentState->bodydata[j].position[d] -
                              CurrentState->bodydata[i].position[d];
          new_i_to_j[d]     = NextState->bodydata[j].position[d] -
                              NextState->bodydata[i].position[d];
          change_in_i_to_j[d] = new_i_to_j[d] - current_i_to_j[d];
        }/*endfor*/;
        if (lengthsquared(change_in_i_to_j) > pow(MAX_REL_DISP, 2.0) *
            lengthsquared(current_i_to_j)) {
          return FALSE;
        }/*endif*/

      }/*endif*/
    }/*endfor*/
  }/*endfor*/;

  return TRUE;
}







/* I/O functions. */

static void ReadState (FILE *filehandle, STATE *AState)
/* specialized to 3-D for now */
{
  AState->Time = 0.0;
  AState->n    = 0;

  while (AState->n < NMAX
         &&
         fscanf(filehandle,
                "%20s %lf  %lf %lf %lf  %lf %lf %lf%*[^\n]\n",
                /* NOTE: 20 in %20s is really NAMELENGTH */
                AState->bodydata[AState->n].name,
                &AState->bodydata[AState->n].mass,

                &AState->bodydata[AState->n].position[0],
                &AState->bodydata[AState->n].position[1],
                &AState->bodydata[AState->n].position[2],

                &AState->bodydata[AState->n].velocity[0],
                &AState->bodydata[AState->n].velocity[1],
                &AState->bodydata[AState->n].velocity[2]) == 8) {
    (AState->n)++;
  }/*endwhile*/
}





static void PutState (STATE *AState, double ATimeStep)
/* treats first arg as read-only */
{
  int i, d;

  printf("STATE NO. %d OF THIS RUN\n", state_count);

  printf("Time=%18g\n", AState->Time);
  printf("Time step to be tried next = %18g\n", ATimeStep);
  printf("Max relative displacement tolerated = %18g\n", MAX_REL_DISP);
  printf("[note: plotting one state in %d, printing one state in %d]\n",
           PLOT_HOW_OFTEN, PRINT_HOW_OFTEN);
  printf("-------------------- ------------------ ------------------ ------------------\n");
  printf("        BODY                MASS             POSITION           VELOCITY     \n");
  printf("-------------------- ------------------ ------------------ ------------------\n");

  for (i=0; i<AState->n; i++) {
    printf("%-20s %18g ", AState->bodydata[i].name, AState->bodydata[i].mass);
    for (d=0; d<NDIMS; d++) {
      printf("%18g %18g\n                                        ",
             AState->bodydata[i].position[d], AState->bodydata[i].velocity[d]);
    }/*endfor*/;
    printf("\n");
  }/*endfor*/;
  PutEnsembleStatistics(AState);
  printf("\n\n\n\n");
}





static void PutEnsembleStatistics (STATE *AState)
/* treats arg as read-only. Potential energy is specialized to 3-D for now. */
{
  double total_mass;
  VECTOR total_moment;
  VECTOR total_momentum;
  double total_potential_energy;
  double total_kinetic_energy;
  VECTOR centre_of_mass_position;
  VECTOR centre_of_mass_velocity;
  double centre_of_mass_virtual_kinetic_energy;
  int    i, j, d;
  VECTOR i_to_j;

  total_mass = 0.0;
  CreateZeroVector(total_moment);
  CreateZeroVector(total_momentum);
  total_potential_energy = 0.0;
  total_kinetic_energy = 0.0;
  for (i=0; i<AState->n; i++) {
    if (AState->bodydata[i].mass != 0.0) {
      total_mass += AState->bodydata[i].mass;
      for (d=0; d<NDIMS; d++) {
        total_moment[d]   += AState->bodydata[i].mass * AState->bodydata[i].position[d];
        total_momentum[d] += AState->bodydata[i].mass * AState->bodydata[i].velocity[d];
      }/*endfor*/;
      for (j=i+1; j<AState->n; j++) {
        for (d=0; d<NDIMS; d++) {
          i_to_j[d] = AState->bodydata[j].position[d] - AState->bodydata[i].position[d];
        }/*endfor*/;
        total_potential_energy +=
         - G * AState->bodydata[i].mass * AState->bodydata[j].mass / sqrt(lengthsquared(i_to_j));
      }/*endfor*/;
      total_kinetic_energy +=
       0.5 * AState->bodydata[i].mass * lengthsquared(AState->bodydata[i].velocity);
    }/*endif*/
  }/*endfor*/

  for (d=0; d<NDIMS; d++) {
    centre_of_mass_position[d] = total_moment[d] / total_mass;
    centre_of_mass_velocity[d] = total_momentum[d] / total_mass;
  }/*endfor*/;
  centre_of_mass_virtual_kinetic_energy =
   lengthsquared(total_momentum) / (2.0 * total_mass);

  printf("ENSEMBLE\n");
  printf(" (AS CENTRE OF MASS) ");
  printf("%18g ", total_mass);
  for (d=0; d<NDIMS; d++) {
    printf("%18g %18g\n                                        ",
            centre_of_mass_position[d], centre_of_mass_velocity[d]);
  }/*endfor*/;
  printf("\n");
  
  printf("FURTHER ENSEMBLE STATISTICS:\n");
  printf(" centre-of-mass virtual kinetic energy = %18g\n\n",
           centre_of_mass_virtual_kinetic_energy);
  printf(" IN STANDARD COORDINATE FRAME:\n");
  printf("  potential energy = %18g\n", total_potential_energy);
  printf("  kinetic energy   = %18g\n", total_kinetic_energy);
  printf("  total energy     = %18g\n\n", total_potential_energy + total_kinetic_energy);
  printf(" IN ENSEMBLE'S REST FRAME:\n");
  printf("  potential energy = %18g\n", total_potential_energy);
  printf("  kinetic energy   = %18g\n", total_kinetic_energy -
                 centre_of_mass_virtual_kinetic_energy);
  printf("  total energy     = %18g\n\n",
            total_potential_energy + total_kinetic_energy -
                 centre_of_mass_virtual_kinetic_energy);
}