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

/*
 *  symp.cpp
 *      Symplectic integration routines to 4th, 6th, and 8th order.
 *  Copyright (c) 1998  David Whysong,  with alterations by Bill Gray.
 *  (Further revised 2013 April 26 by Bill Gray:  comments added,
 *  'const' declarations added,  more digits shown in results,  revised
 *  to compile cleanly in Linux.)
 *
 *   This program is free software; you can redistribute it and/or modify
 *   it under the terms of the GNU General Public License as published by
 *   the Free Software Foundation, version 2.
 *
 * Reference: "Construction of higher order symplectic integrators" Haruo
 *             Yoshida, Phys. Lett. A 150, pp. 262, 12 Nov. 1990.
 *    http://cacs.usc.edu/education/phys516/Yoshida-symplectic-PLA00.pdf
 *
 *  For the test program,  compile as:
 *
 *  cl -W4 -Ox -DTEST_CODE symp.cpp              (Microsoft C/C++)
 *  wcl386 -W4 -Ox -DTEST_CODE symp.cpp          (Watcom C/C++)
 *  g++ -Wall -o symp -O3 -DTEST_CODE symp.cpp
 *
 *    The test case comes from p. 168,  J M A Danby's _Fundamentals of
 * Celestial Mechanics_;  it corresponds to an object making about
 * one-sixth of a roughly circular orbit.  The test code shows the
 * ending state vector,  plus the error (numerical solution minus
 * analytical solution).  You run the test program with a number of
 * steps and a '4',  '6',  or '8' to indicate the method to be used;
 * for example,
 *
 * symp 19 6
 *
 * to use the symplectic_6 routine,  with 19 integration steps.  You'll
 * notice that doubling the number of steps cuts the error for
 * symplectic_4 about 16-fold;  for symplectic_6,  64-fold;  and for
 * symplectic_8,  about 256-fold (which matches the fact that these
 * routines should be fourth,  sixth,  and eighth-order,  respectively.)
 *
 *    NOTE that at first blush,  the idea of using a sixth or eighth-order
 * integrator is going to look wonderful.  However,  the fourth-order method
 * uses three force-function evaluations;  the sixth-order uses seven;  and
 * the eighth-order uses 15.  Thus,  for example,  all of these :
 *
 * ./symp 28 8
 * ./symp 140 4
 * ./symp 60 6
 *
 *   should require 28*15 = 140*3 = 60*7 = 420 force-function evaluations,
 * and will therefore take about the same amount of computational time.
 * At least in this case,  the sixth-order method provides smaller errors
 * than the other two.  But if you needed higher accuracy (and therefore
 * increased the number of steps),  eventually,  the eighth-order method
 * would be the winner;  and for a sufficiently rough evaluation,  you
 * would find that the fourth-order method worked best.
 *
 *    The 'acceleration' code takes the time and velocity as parameters,
 * but in this particular case,  the acceleration is a function of
 * position only.  I'm pretty sure I've got the time and velocity
 * dependence right,  but can't say for an absolute certainty.  If
 * one wanted to include planetary perturbations (which vary with
 * time) or some sort of drag term (a function of velocity),  this
 * would suddenly matter.
 *
 *   I'll throw in some comments on why this is useful stuff,  derived
 * mostly from e-mail talks with David Whysong.  The big plus with a
 * symplectic method is that quantities such as angular momentum and
 * total energy are pretty well-conserved (more so than with other
 * methods such as Runge-Kutta,  Adams-Bashforth,  Adams-Bashforth-
 * Moulton,  and their kin.)  David says they have only really emerged
 * in the last decade or so,  which explains why I never heard of them
 * in physics courses in 1983-87.  He's using them in research of
 * globular cluster simulations,  where it's necessary to run _very_
 * long integrations.  If,  say,  energy was not conserved over so long
 * an integration,  accumulated errors might cause your simulated
 * globular cluster to undergo a simulated explosion,  or a simulated
 * gravitational collapse.  I _don't_ know if a price is paid in
 * precision for this benefit.  I also have found no references to
 * symplectic integration.
 *
 */

void compute_accel( double *accel, const double *x, const double *v,
                                 const double t)
{
   int i;
   const double r2 = x[0] * x[0] + x[1] * x[1] + x[2] * x[2];
   const double a0 = 1. / (r2 * sqrt( r2));

   for( i = 0; i < 3; i++)
      accel[i] = -x[i] * a0;
}

void symplectic_4( double *x, double *v, double t, const double dt)
{
   int i, j;
   const double b = 1.25992104989487319066654436028;      /* 2^1/3 */
   const double a = 2 - b;
   const double x0 = -b / a;
   const double x1 = 1. / a;
   const static double d4[3] = { x1, x0, x1 };
   const static double c4[4] = { x1/2, (x0+x1)/2, (x0+x1)/2, x1/2 };

   for( i = 0; i < 4; i++)
      {
      double accel[3];
      const double step = dt * c4[i];

      for( j = 0; j < 3; j++)
         x[j] += step * v[j];
      t += step;
      if( i != 3)
         {
         compute_accel( accel, x, v, t);
         for( j = 0; j < 3; j++)
            v[j] += dt * d4[i] * accel[j];
         }
      }
}

void symplectic_6( double *x, double *v, double t, const double dt)
{
   int i, j;
   const double w1 = -0.117767998417887E1;
   const double w2 = 0.235573213359357E0;
   const double w3 = 0.784513610477560E0;
   const double w0 = (1-2*(w1+w2+w3));
   const static double d6[7] = { w3, w2, w1, w0, w1, w2, w3 };
   const static double c6[8] = { w3/2, (w3+w2)/2, (w2+w1)/2, (w1+w0)/2,
                         (w1+w0)/2, (w2+w1)/2, (w3+w2)/2, w3/2 };

   for( i = 0; i < 8; i++)
      {
      double accel[3];
      const double step = dt * c6[i];

      for( j = 0; j < 3; j++)
         x[j] += step * v[j];
      t += step;
      if( i != 7)
         {
         compute_accel( accel, x, v, t);
         for( j = 0; j < 3; j++)
            v[j] += dt * d6[i] * accel[j];
         }
      }
}


void symplectic_8( double *x, double *v, double t, const double dt)
{
   int i, j;
   const double W1 = -0.161582374150097E1;
   const double W2 = -0.244699182370524E1;
   const double W3 = -0.716989419708120E-2;
   const double W4 =  0.244002732616735E1;
   const double W5 =  0.157739928123617E0;
   const double W6 =  0.182020630970714E1;
   const double W7 =  0.104242620869991E1;
   const double W0 = (1-2*(W1+W2+W3+W4+W5+W6+W7));
   const static double d8[15] = { W7, W6, W5, W4, W3, W2, W1, W0,
                                  W1, W2, W3, W4, W5, W6, W7 };
   const static double c8[16] = { W7/2, (W7+W6)/2, (W6+W5)/2, (W5+W4)/2,
                         (W4+W3)/2, (W3+W2)/2, (W2+W1)/2, (W1+W0)/2,
                         (W1+W0)/2, (W2+W1)/2, (W3+W2)/2, (W4+W3)/2,
                         (W5+W4)/2, (W6+W5)/2, (W7+W6)/2,  W7/2 };

   for( i = 0; i < 16; i++)
      {
      double accel[3];
      const double step = dt * c8[i];

      for( j = 0; j < 3; j++)
         x[j] += step * v[j];
      t += step;
      if( i != 15)
         {
         compute_accel( accel, x, v, t);
         for( j = 0; j < 3; j++)
            v[j] += dt * d8[i] * accel[j];
         }
      }
}

#ifdef TEST_CODE
int main( const int argc, const char **argv)
{
   double x[6] = {1., .1, -.1, -.1, 1., .2};
   const double analytic[6] = {
          .4660807846539234, .9006112349077236, .1140459806673234,
         -.8765944368525084, .4731566049794786, .1931594924439623 };
   const int n_steps = atoi( argv[1]);
   int i;

   for( i = 0; i < n_steps; i++)
      switch( argv[2][0])
         {
         case '4':
            symplectic_4( x, x + 3, 0., 1. / (double)n_steps);
            break;
         case '6':
            symplectic_6( x, x + 3, 0., 1. / (double)n_steps);
            break;
         case '8':
            symplectic_8( x, x + 3, 0., 1. / (double)n_steps);
            break;
         default:
            printf( "Option '%s' not understood\n", argv[2]);
            return( -1);
            break;
         }
   printf( "position:  %18.15lf %18.15lf %18.15lf \n", x[0], x[1], x[2]);
   printf( "velocity:  %18.15lf %18.15lf %18.15lf \n", x[3], x[4], x[5]);
   printf( "posn err:  %18.15lf %18.15lf %18.15lf \n",
         x[0] - analytic[0], x[1] - analytic[1], x[2] - analytic[2]);
   printf( "vel err:   %18.15lf %18.15lf %18.15lf \n",
         x[3] - analytic[3], x[4] - analytic[4], x[5] - analytic[5]);
   return( 0);
}
#endif