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



struct par{
	double mu;
	double om;
	double de;
};

int
func (double t, const double y[], double f[], struct par *pz)
{ 
  
  double mu=pz->mu;
  double om=pz->om;
  double de=pz->de;
  //printf ("%.5e %.5e %.5e\n", pz->mu, pz->om, pz->de);
  
  f[0] = om*y[1];
  f[1] = -y[0] - mu*y[1]*(y[0]*y[0] - 1);
  f[2] = de*y[1]*y[2];
  
  return GSL_SUCCESS;
}

int
jac (double t, const double y[], double *dfdy, 
double dfdt[], struct par *pz)
{
  double mu=pz->mu;
  double om=pz->om;
  double de=pz->de;
  gsl_matrix_view dfdy_mat 
    = gsl_matrix_view_array (dfdy, 3, 3);
  gsl_matrix * m = &dfdy_mat.matrix; 
  gsl_matrix_set (m, 0, 0, 0.0);
  gsl_matrix_set (m, 0, 1, om);
  gsl_matrix_set (m, 0, 2, 0.0);
  gsl_matrix_set (m, 1, 0, -2.0*mu*y[0]*y[1] - 1.0);
  gsl_matrix_set (m, 1, 1, -mu*(y[0]*y[0] - 1.0));
  gsl_matrix_set (m, 1, 2, 0.0);    
  gsl_matrix_set (m, 2, 0, 0.0);
  gsl_matrix_set (m, 2, 1, de*y[2]);
  gsl_matrix_set (m, 2, 2, de*y[1]);
  
  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  dfdt[2] = 0.0;
  return GSL_SUCCESS;
}
/*
void print_para(struct par *pz)
{
	printf ("%.5e %.5e %.5e\n", pz->mu, pz->om, pz->de);
}
*/

int
main (void)
{
  FILE *fp;
  fp = fopen("odetest.dat","w");
  if(fp == NULL)
   {
      fprintf(stderr,"error: can't open odetest.dat!\n");
      exit(EXIT_FAILURE);
   }
   /*
   time_t start, stop;
   double diff;  
   start = time(NULL);
   */
   
  const gsl_odeiv_step_type * T 
    = gsl_odeiv_step_rk8pd;

  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);

  
  struct par z, *pz;
  pz=&z;
  z.mu = 10;
  pz->om = 1.45;
  pz->de = 0.2;
  
  printf ("%.5e %.5e %.5e\n", pz->mu, pz->om, pz->de);
  
  //print_para(pz);
  
  gsl_odeiv_system sys = {func, jac, 3, pz};

  double t = 0.0, t1 = 100.0;
  double h = 1e-6;
  double y[3] = { 2.0, 2.0, 2.0 }; //initial cond

  while (t < t1)
    {
 	int status = gsl_odeiv_evolve_apply (e, c, s, &sys, &t, t1, &h, y);

	 if (status != GSL_SUCCESS)
		break;

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

  gsl_odeiv_evolve_free (e);
  gsl_odeiv_control_free (c);
  gsl_odeiv_step_free (s);
  /*
   stop = time(NULL);
   diff = difftime(start, stop);

   printf("Difference is %e\n", diff);
  */
  return 0;
}