#include <stdio.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_odeiv.h>
     

typedef struct bz_data { 
  int count;
  int jac_count;
  double alpha; 
  double beta; 
  double gamma;
} bz_data; 

int bz (double t, const double y[], double f[], void *params){

  bz_data * b = (bz_data *) params; 

  f[0] = b->alpha*(y[1] + y[0]*(1-b->beta*y[0] - y[1]));
  f[1] = (1/b->alpha)*(y[2] - (1+y[0])*y[1]); 
  f[2] = b->gamma*(y[0] - y[2]); 
  ++b->count;

  return GSL_SUCCESS;
}

int jac_bz (double t, const double y[], double *dfdy, double dfdt[], void *params) {

  bz_data *p = (bz_data *)params;
  double alpha = p->alpha;
  double beta = p->beta;
  double gamma = p->gamma;
  gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 3, 3);
  gsl_matrix * m = &dfdy_mat.matrix; 

  p->jac_count ++ ; 

  gsl_matrix_set (m, 0, 0, alpha*(1-beta*y[0]-y[1] - beta*y[0]));
  gsl_matrix_set (m, 0, 1, alpha*(1 - y[0]));
  gsl_matrix_set (m, 0, 2, 0);
  
  gsl_matrix_set (m, 1, 0, (1/alpha)*(-y[1]));
  gsl_matrix_set (m, 1, 1, -(1/alpha)*(1+y[0]));
  gsl_matrix_set (m, 1, 2, (1/alpha));
  
  gsl_matrix_set (m, 2, 0, gamma*y[0]);
  gsl_matrix_set (m, 2, 1, 0);
  gsl_matrix_set (m, 2, 2, -gamma*y[2]);
  
  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  dfdt[2] = 0.0;

  return GSL_SUCCESS;
}

int main () {
  //const gsl_odeiv_step_type * T = gsl_odeiv_step_rk2imp;
  const gsl_odeiv_step_type * T = gsl_odeiv_step_rk4imp;
  //const gsl_odeiv_step_type * T = gsl_odeiv_step_rkf45;
  //const gsl_odeiv_step_type * T =  gsl_odeiv_step_bsimp;
  //const gsl_odeiv_step_type * T =  gsl_odeiv_step_gear2;
  gsl_odeiv_step * s = gsl_odeiv_step_alloc (T, 3);
  gsl_odeiv_control * c = gsl_odeiv_control_y_new (1e-6, 0.0);
  gsl_odeiv_evolve * e = gsl_odeiv_evolve_alloc (3);
  bz_data dat = {0, 0, 77.27, 8.375E-6,  1.161}; 
  gsl_odeiv_system sys = {bz, jac_bz, 3, &dat};

     
  double t = 0.0, t1 = 360;
  double h = 1e-6;


  double y[] = { 1, 2, 3};

     
  double t2 = t; 
  double interval = 0.05; // suppress excessive output
  while (t < t1)
  {
    int status = gsl_odeiv_evolve_apply (e, c, s, &sys, &t, t1, &h, y);
    if (status != GSL_SUCCESS)
      break;
    if (t > t2 + interval ) {
      printf ("%.5e %.5e %.5e %5e  %d\n", t, y[0], y[1], y[2], dat.count);
      t2 = t; 
    }
  }
     
  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (s);

  printf("Number of Jacobian evaluations = %d\n"
         "Number of Function evaluations = %d\n", dat.jac_count,
  dat.count);

  return 0;
}