/*
  This program calculates stroboscopic map for a driven pendulum.
  Equation of motion is: y'' + y'/Q + sin(y) = r cos(at)

  The output is stored in file datapoints.dat

  To compile this code on a GNU or UNIX platform you should type:
     gcc pendulum.c -lgsl -lgslcblas -lm
  This assumes that you have gsl libraries on your library path and 
  gsl headers on your include path. For other platforms (such as Windows)
  please consult your compiler documentation.
*/

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

#define ODE_DIM 3

struct ODE_params
{
  double Q, drive_freq, r;
};

struct PCR_params
{
  int STEPS_PER_SECTION;
  double STROBE_PERIOD;
  struct ODE_params *odeparams;
};

int
func (double t, const double y[], double f[],
      void *params)
{
  struct ODE_params *p 
    = (struct ODE_params *) params;

  double Q = p->Q;
  double a = p->drive_freq;
  double r = p->r;

  f[0] = y[1];
  f[1] = - y[1]/Q -sin(y[0]) + r*cos(y[2]);
  f[2] = a;
  return GSL_SUCCESS;
}

int
jac (double t, const double y[], double *dfdy,
     double dfdt[], void *params)
{
  struct ODE_params *p 
    = (struct ODE_params *) params;

  double Q = p->Q;
  double a = p->drive_freq;
  double r = p->r;

  gsl_matrix_view dfdy_mat
    = gsl_matrix_view_array (dfdy, ODE_DIM, ODE_DIM); 
  gsl_matrix * m = &dfdy_mat.matrix;
  gsl_matrix_set (m, 0, 0, 0.0);
  gsl_matrix_set (m, 0, 1, 1.0);
  gsl_matrix_set (m, 0, 2, 0.0);
  gsl_matrix_set (m, 1, 0, -cos(y[0]));
  gsl_matrix_set (m, 1, 1, -1.0/Q);
  gsl_matrix_set (m, 1, 2, -r*sin(y[2]));
  gsl_matrix_set (m, 2, 0, 0.0);
  gsl_matrix_set (m, 2, 1, 0.0);
  gsl_matrix_set (m, 2, 2, 0.0);
  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  dfdt[2] = 0.0;
  return GSL_SUCCESS;
}

int
Poincare_map_strob (double *y, void *params)
{
  int i,j;
  const gsl_odeiv_step_type * T
    = gsl_odeiv_step_rk8pd;

  gsl_odeiv_step * s
    = gsl_odeiv_step_alloc (T, ODE_DIM);
  gsl_odeiv_control * c
    = gsl_odeiv_control_y_new (1e-12, 0.0);
  gsl_odeiv_evolve * e
    = gsl_odeiv_evolve_alloc (ODE_DIM);

  struct PCR_params *p 
    = (struct PCR_params *) params;
  int N 
    = p->STEPS_PER_SECTION;
  double Tp 
    = p->STROBE_PERIOD;

  double tin=0.0, deltat=Tp/N;
  double t, t1;
  double h = 1e-6;

  gsl_odeiv_system sys = {func, jac, ODE_DIM, p->odeparams};

  gsl_ieee_env_setup();

  t = tin;


  do
    {
      t1 = t + deltat;
      if ( t1 > Tp ) t1 = Tp;

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

      t = t1;
    } 
  while ( t < Tp );

  gsl_odeiv_evolve_free(e);
  gsl_odeiv_control_free(c);
  gsl_odeiv_step_free(s);
  return 0;
} 



int
main(void)
{
  int i,j,k;
  const int NITER=500; 
  FILE *tocke, *fopen();

  /* Initial conditions */
  double y[ODE_DIM]={ 0.0, 1.0, 0.0 };

  /* Set ODE parameters */
  struct ODE_params rod_param={ 2.0,  /*  Quality factor "Q" */
				2.0/3.0, /*  Drive frequency "a" */
				1.5 }; /*  Drive amplitude "r" */

  /* Set Poincare map parameters */
  struct PCR_params rodpcr_param={200, /* Integration steps per section */
				  2.0*M_PI/rod_param.drive_freq, /* strobe period */
				  &rod_param}; /* ODE parameters */

  tocke=fopen("datapoints.dat","w");

  for ( j=0; j < NITER; j++)
    {
      if (fabs(y[0]) > M_PI) y[0] -= 2.*M_PI*y[0]/fabs(y[0]);
      fprintf(tocke,"%g %g\n", y[0], y[1]); 
      /* printf("%g %g\n", y[0], y[1]);  */
      Poincare_map_strob(y,&rodpcr_param);

    }

  fclose(tocke);
  return 0;
}