#include <stdio.h>
#include <string.h>

#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_odeiv2.h>

#define DBL_ZERO_MEMSET(dest,n) memset((dest),0,(n)*sizeof(double))

/* Henon-Heiles Problem */

static int
rhs_henon (double t, const double y[], double f[], void *p)
{
  (void) t;
  (void) p;

  f[0] = -y[2] - 2.0 * y[2] * y[3];
  f[1] = -y[3] - y[2] * y[2] + y[3] * y[3];
  f[2] = y[0];
  f[3] = y[1];

  return GSL_SUCCESS;
}

static int
jac_henon (double t, const double y[], double *dfdy, double dfdt[], void *p)
{
  (void) t;
  (void) p;

  DBL_ZERO_MEMSET (dfdy, 16);

  dfdy[2] = -1.0 - 2.0 * y[3];
  dfdy[3] = -2.0 * y[2];
  dfdy[6] = dfdy[3];
  dfdy[7] = -1.0 + 2.0 * y[3];
  dfdy[8] = 1.0;
  dfdy[13] = 1.0;

  DBL_ZERO_MEMSET (dfdt, 4);

  return GSL_SUCCESS;
}

static const gsl_odeiv2_system sys_henon = { rhs_henon, jac_henon, 4, NULL };

static int
int_henon (double t, const double y[], double I[], void *p)
{
  (void) t;
  (void) p;

  I[0] = 0.0;
  for (size_t i = 0; i < 4; i++)
    {
      I[0] += y[i] * y[i];
    }
  I[0] /= 2;

  I[0] += (y[2] * y[2] - y[3] * y[3] / 3.0) * y[3];

  return GSL_SUCCESS;
}

int
main (int argc, char *argv[])
{
  double t = 0.0, h = 0.1, t1 = 1e4;
  double y[4] = {-0.1 -0.07, 0.05, -0.01};

  double y_err[4], dydt_in[4], dydt_out[4];

  int ret = 0;

  const gsl_odeiv2_step_type *T = gsl_odeiv2_step_imprk2;
  gsl_odeiv2_step *s = gsl_odeiv2_step_alloc (T, 4);
  gsl_odeiv2_control *c = gsl_odeiv2_control_y_new (1e-10, 0.0);
  gsl_odeiv2_step_set_control (s, c);

  /* initialise dydt_in from system parameters */
  GSL_ODEIV_FN_EVAL(&sys_henon, t, y, dydt_in);

  do
    {
      {
        double H;

        int_henon (t, y, &H, 0);

        static int first = 1;
        static double H0;

        if (first)
          {
            H0 = H;

            first = !first;
          }

        fprintf (stdout, "%le\t%+le\n", t, H - H0);
      }

      {
        int status = gsl_odeiv2_step_apply (s, t, h,
                                           y, y_err,
                                           dydt_in,
                                           dydt_out,
                                           &sys_henon);

        if (status != GSL_SUCCESS)
          {
            ret = status;

			fprintf (stderr,"%d\n",status);

            break;
          }

        memcpy (dydt_in, dydt_out, 4 * sizeof (double));

        t += h;
      }
    }
  while (t < t1);

  gsl_odeiv2_step_free (s);
  gsl_odeiv2_control_free (c);

  return ret;
}

/*
  Local Variables:
  compile-command: "c99 -lgsl -lgslcblas bug.c -o bug"
  End:
*/