/*
    vanderpol_using_gsl.c, by Timothy Jones, models the Van der Pol equations via gsl 

    Copyright (C) 2010 Timothy Jones

    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/


#include <stdio.h>
#include <gsl/gsl_errno.h>  // GSL_SUCCESS/ERROR 
#include <gsl/gsl_matrix.h> // GSL Matrix
#include <gsl/gsl_math.h>   // basic math functions
#include <gsl/gsl_odeiv.h>  // ODE solvers

 /*    printf("\n");
       printf("vanderpol_using_gsl.c written on Sep 1 2010 by Tim Jones @ Drexel University.");
       printf("\n"); printf("\n");\
      
      

       printf("You are running a.out compiled via gcc vanderpol_using_gsl.c -lgsl -lgslcblas.\n");
       printf("\n"); printf("\n");
       printf("This program integrates the driven van der Pol oscillator using rk8pd from GSL.\n");
       printf("\n"); printf("\n");
       */

/* Dimension of ODEs */
#define DIM 3

static  double A = 0.25;
static  double b=0.7;
static  double c=1.0;
static  double d=10.0;

 int
     func (double t, const double y[], double f[],
           void *params)
     {
       
       double mu = *(double *)params;
       double w=mu;
       f[0] = b*y[1] + (c-d*y[1]*y[1])*y[0];

       f[1] = -y[0] +  A*sin(w*y[2]);
       
       f[2] = 1;

       return GSL_SUCCESS;
     }


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

       double mu = *(double *)params;
       double w=mu;
       gsl_matrix_view dfdy_mat 
         = gsl_matrix_view_array (dfdy, 3, 3);
       gsl_matrix * m = &dfdy_mat.matrix; 
       gsl_matrix_set (m, 0, 0,     (c-d*y[1]*y[1])    );
       gsl_matrix_set (m, 0, 1,     b + 2*d*y[1]*y[0]  );
       gsl_matrix_set (m, 0, 2,       0                );
       gsl_matrix_set (m, 1, 0,      -1                );
       gsl_matrix_set (m, 1, 1,       0                );
       gsl_matrix_set (m, 1, 2,     A*cos(w*y[2])      );
       gsl_matrix_set (m, 2, 0,       0                );
       gsl_matrix_set (m, 2, 1,       0                );
       gsl_matrix_set (m, 2, 2,       0                );
       dfdt[0] = 0.0;
       dfdt[1] = 0.0;
       dfdt[2] = 0.0;
       return GSL_SUCCESS;
     }


 int
     main (void)
     {
       
       double mu=M_PI/2;  // 2*M_PI; 
       system("mv slices2 slices_old; mkdir slices2");

       char temp[512];
       char temp2[1512];
       //system((char *)temp);
       FILE * fd;
       sprintf(temp,"output.dat");
       fd = fopen(temp, "w");
      

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


       gsl_odeiv_system sys = {func, jac, 3, &mu};
     
       double t = 0.0, t1 = 3000.0;
       double h = 1e-6;
       double y[3] = { 0.1, 0.1, 0.0 };
     
       int i=1;
       int ii;

       double tt;
       double ttold;
       
       int LOOPNUM; 
       int DATAPOINTS=200;  // Data points per cycle
       LOOPNUM=t1*DATAPOINTS + 4*DATAPOINTS;//HOW FAR LOOP NEEDS TO GO + compensation for transients.

       FILE * fa;
       int filecount=0;
       double ti=0;

       for (i = 1; i <= (LOOPNUM); i++)
	{
	  //	  double ti = i * t1 / (1.0*LOOPNUM);
          ti = ti + 0.02;  // dt = 4/200, this should be put into the code.
	  //printf("%f\n",ti);
	  while (t < ti)
	    {
	      gsl_odeiv_evolve_apply (e, c, s, 
                                  &sys, 
                                  &t, ti, &h,
                                  y);
	    }
	  if(ti>=4.0){ti=0;t=0;}
          
	  
	  //	  ttold=tt;
	  tt = t;// - floor(t);

	  //if((ttold-tt)>0){filecount=0; t=0.0;}
	  if(t==0){filecount=0; t=0.0;}
         
	  //skip the first set to remove transients. 
	  ii=filecount+100000;
	  if(i>4*DATAPOINTS){  sprintf(temp2, "./slices2/%i.dat",ii);
	                       fa=fopen(temp2, "a");
	                       fprintf(fa, "%.5e %.5e\n", y[0], y[1]); 
			       fclose(fa);
	                       fprintf(fd, "%.5e %.5e %.5e %.5e\n", 2*M_PI*tt, y[0], y[1], y[2]);}
	  filecount++;
	 
    }
 
      printf("\n");
      printf("The program cycled through %i iterations of %f.\n",(int)(floor(t)), 2*M_PI);
      printf("\n"); printf("\n");
      printf("Each of the %i iterations contained %i data points.\n",(int)(floor(t)), DATAPOINTS);
      printf("\n"); printf("\n");
      printf("Output file: %s with four columns containing: theta, y[0] y[1] and y[2]=t.\n",temp);
      printf("\n"); printf("\n");
 
       gsl_odeiv_evolve_free (e);
       gsl_odeiv_control_free (c);
       gsl_odeiv_step_free (s);
       return 0;
     }