/*
 * The program solves two first-order nonlinear 
 * diffrential equations aka Lotka-Volterra predator-prey model 
 *
 *  r' = 2r - \alpha r f
 *  f' = -f + \alpha r f
 *
 * where \alpha is a constant
 *
 * Invariant of the model: I = log(f^2*r) - \alpha(r + f)
 *
 * The program begins by defining the functions for the derivatives.
 * The program evolves the solution from (r_0, f_0) at t = 0.
 * The step-size h is automatically adjusted by the controller 
 * to maintain an absolute accuracy of 10^{-6} in the function 
 * values (r, f). 
 */

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

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

int func (double t, const double y[], double f[], void *params);
double invariant(const double y[], const double alpha);
 
int func (double t, const double y[], double f[], void *params) 
{
  (void) t;
  double alpha = *(double *) params;
  double arf = alpha * y[0] * y[1];      /* \alpha*r*f */
  f[0] = 2*y[0] - arf;
  f[1] = -y[1] + arf;

  return GSL_SUCCESS;
}

double invariant(const double y[], const double alpha) 
{
  return (log(y[1]*y[1]*y[0]) - alpha*(y[0] + y[1]));
}

int main (void) 
{
  size_t neqs = 2;                       /* number of equations */
  double eps_abs = 1.e-9, eps_rel = 0.;  /* desired precision */
  double stepsize = 1e-6;                /* initial integration step size */
  double alpha = .01;                    /* the parameter of the model */
  double t = 0., t1 = 50.;               /* time interval of the evolution */
  int status;

  /*
   * Explicit embedded Runge-Kutta-Fehlberg (4,5) method. 
   * This method is a good general-purpose integrator.
   */
  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);
    
  gsl_odeiv2_system sys = {func, NULL, neqs, &alpha};
    
  /* 
   * Initial conditions 
   */
  double y[2] = { 300., 150. };

  printf ("% .5e % .5e % .5e % .10e\n", t, y[0], y[1], invariant(y, alpha));
  /*
   * Evolution loop 
   */
  while (t < t1) {
      status = gsl_odeiv2_evolve_apply (e, c, s, &sys, &t, t1, &stepsize, y);

      if (status != GSL_SUCCESS) {
          printf ("Troubles at % .5e % .5e % .5e\n", t, y[0], y[1]);
          break;
      }

      printf ("% .5e % .5e % .5e % .10e\n", t, y[0], y[1], invariant(y, alpha));
  }

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

  return 0;
}