/*
 * A simplified mathematical model for atmospheric convection
 * known as Lorenz model
 *
 * dx/dt = sigma*(y - x)
 * dy/dt = x*(rho - z) - y
 * dz/dt = x*y - beta*z
 *
 * It is notable for having chaotic solutions for certain parameter 
 * values and initial conditions.
 *
 * The step-size h is automatically adjusted by the controller 
 * to maintain an absolute accuracy of 10^{-6} in the function 
 * values.
 */

#include <stdio.h>

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

int func (double t, const double y[], double f[], void *params) {
  double sigma = ((double *) params)[0];
  double rho = ((double *) params)[1];
  double beta = ((double *) params)[2];
  f[0] = sigma*(y[1] - y[0]);
  f[1] = y[0]*(rho - y[2]) - y[1];
  f[2] = y[0]*y[1] - beta*y[2];
  return GSL_SUCCESS;
}

int main (void) {
  size_t neqs = 3;                       /* number of equations */
  double eps_abs = 1.e-6, eps_rel = 0.;  /* desired precision */
  double stepsize = 1e-6;                /* initial integration step size */
  double sigma = 10., rho = 28., beta = 8./3.;   /* the parameters of the system */
  double t = 0., t1 = 50.;              /* time interval of the evolution */
  int status;
  double params[3];

  params[0] = sigma;
  params[1] = rho;
  params[2] = beta;

  /*
   * 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, params};
    
  /* 
   * Initial conditions 
   */
  double y[3] = { 0., 1., 0. };

  /*
   * 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 % .5e\n", t, y[0], y[1], y[2]);
          break;
      }

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

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

  return 0;
}