/*
 * The program solves the restricted three body problem
 *
 * x'' = -(1-\alpha)(x+\alpha)/((x+\alpha)^2+y^2)^{3/2}
 *       -\alpha(x+\alpha-1)/((x+\alpha-1)^2+y^2)^{3/2}
 *       + x + 2 y'
 * y'' = -(1-\alpha)y/((x+\alpha)^2+y^2)^{3/2}
 *       -\alpha y/((x+\alpha-1)^2+y^2)^{3/2}
 *       + y - 2 x'
 *
 * This can be converted into a first order system 
 * by introducing a separate variable for the velocities, 
 * vx = x'(t) and vy = y'(t)
 *
 *   x'  = vx
 *   y'  = vy
 *   vx' = ...
 *   vy' = ...
 *
 * The program begins by defining functions for these 
 * derivatives. The program evolves the solution for the case a = .5
 * (equal masses of the heavy bodies) from (u, v) = (1, 0) at t=0 to t=100. 
 * The step-size h is automatically adjusted by the controller 
 * to maintain an absolute accuracy of 10^{-6} in the function 
 * values (u, v). 
 */

#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv2.h>

int
func (double t, const double yy[], double f[], void *params)
{
  double a = *(double *) params;
  double d1, d2;
  double x = yy[0], y = yy[1], vx = yy[2], vy = yy[3];

  d1 = pow((x + a)*(x + a) + y*y, 1.5);
  d2 = pow((x + a - 1.)*(x + a - 1.) + y*y, 1.5);

  f[0] = vx;
  f[1] = vy;
  f[2] = -(1. - a)*(x + a)/d1 - a*(x + a - 1.)/d2 + x + 2*vy,
  f[3] = -(1. - a)*y/d1 - a*y/d2 + y - 2*vx;

  return GSL_SUCCESS;
}

int
main (void)
{
  size_t neqs = 4;
  double eps_abs = 1.e-6, eps_rel = 0.;

  gsl_odeiv2_step    *s = gsl_odeiv2_step_alloc (gsl_odeiv2_step_rkf45, neqs);
  gsl_odeiv2_control *c = gsl_odeiv2_control_y_new (eps_abs, eps_rel);
  gsl_odeiv2_evolve  *e = gsl_odeiv2_evolve_alloc (neqs);
    
  double a = 0.5;
  gsl_odeiv2_system sys = {func, NULL, neqs, &a};
    
  double t = 0., t1 = 100.;
  double stepsize = 1e-6;
  double yy[4] = { 1., 0., 0., 0. };
    
  while (t < t1) {
      int status = gsl_odeiv2_evolve_apply (e, c, s, &sys, &t, t1, &stepsize, yy);

      if (status != GSL_SUCCESS)
          break;

      printf ("% .5e  % .5e  % .5e  % .5e  % .5e\n", t, yy[0], yy[1], yy[2], yy[3]);
  }

  gsl_odeiv2_evolve_free (e);
  gsl_odeiv2_control_free (c);
  gsl_odeiv2_step_free (s);

  return 0;
}