# include <ctime>
# include <fstream>
# include <iomanip>
# include <iostream>

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

using namespace std;

int main ( );
int vanderpol_func ( double t, const double y[], double f[], void *params );
int vanderpol_jac ( double t, const double y[], double *dfdy, double dfdt[], 
  void *params );
void timestamp ( );

//****************************************************************************80

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    midpoint_gsl_test() demonstrates the GSL midpoint ODE solver.
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    10 November 2023
//
//  Author:
//
//    Taken from the GSL reference manual.
//    This version by John Burkardt.
//
{
  double dt;
  ofstream data;
  string data_filename = "midpoint_gsl_data.txt";
  double epsabs = 1e-6;
  double epsrel = 0.0;
  double h_start = 1e-6;
  int i;
  double mu = 10.0;
  int n = 100;
  double t = 0.0;
  double t0 = 0.0;
  double tstop = 100.0;
  double y[2] = { 1.0, 0.0 };
//
//  Define "sys", the definition of the ODE data.
//
  gsl_odeiv2_system sys = { vanderpol_func, vanderpol_jac, 2, &mu };
//
//  Define "d", a pointer to the ODE driver to be used.
//  Here, "gsl_odeiv2_step_rk2imp" specifies we want to use the GSL midpoint driver.
//
  gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new ( 
    &sys, gsl_odeiv2_step_rk2imp, h_start, epsabs, epsrel );

  timestamp ( );
  cout << "\n";
  cout << "midpoint_gsl_test():\n";
  cout << "  C++ version\n";
  cout << "  gsl_odeiv2_step_rk2imp() solves an ODE.\n";

  cout << "\n";
  cout << "  parameters:\n";
  cout << "    mu    = " << mu << "\n";
  cout << "    t0    = " << t0 << "\n";
  cout << "    y0    = (" << y[0] << "," << y[1] << "\n";
  cout << "    tstop = " << tstop << "\n";
  cout << "    n     = " << n << "\n";

  data.open ( data_filename.c_str ( ) );

  cout << "  " << t 
       << "  " << y[0]
       << "  " << y[1] << "\n";
  data << "  " << t 
       << "  " << y[0]
       << "  " << y[1] << "\n";

  dt = ( tstop - t0 ) / n;

  for ( i = 1; i <= n; i++ )
  {
    double ti = i * dt;
    int status = gsl_odeiv2_driver_apply ( d, &t, ti, y );
    if ( status != GSL_SUCCESS )
    {
      cout << "error, return value = " << status << "\n";
      break;
    }
    cout << "  " << t 
         << "  " << y[0]
         << "  " << y[1] << "\n";
    data << "  " << t 
         << "  " << y[0]
         << "  " << y[1] << "\n";
  }

  data.close();
  cout << "\n";
  cout << "  Solution data stored in '" << data_filename << "'.\n";
//
//  Free memory.
//
  gsl_odeiv2_driver_free ( d );
//
//  Terminate.
//
  cout << "\n";
  cout << "midpoint_gsl_test():\n";
  cout << "  Normal end of execution.\n";
  timestamp ( );

  return 0;
}
//****************************************************************************80

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

//****************************************************************************80
//
//  Purpose:
//
//    vanderpol_func() evaluates the right hand side of vanderpol_ode().
//
//  Discussion:
//
//    This function defines the right hand side of the van der Pol equation:
//
//      u'' + mu u' (u^2-1) + u = 0
//
//    rewritten in first order form as:
//
//      u' = v
//      v' = - u + mu v ( 1 - u^2)
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    10 November 2023
//
//  Author:
//
//    Taken from the GSL reference manual.
//    This version by John Burkardt.
//
//  Input:
//
//    double T: the current time.
//
//    const double Y[]: the current solution value.
//
//    void *PARAMS: user-specified parameters.
//
//  Output:
//
//    int FUNC: returns GSL_SUCCESS indicating a successful computation.
//
//    double F[]: the right hand side of the ODE at T, Y.
//
{
//
//  Declare t to avoid unused parameter warning.
//
  (void) (t);
  double mu = *(double *) params;

  f[0] = y[1];
  f[1] = - y[0] - mu * y[1] * ( y[0] * y[0] - 1.0 );

  return GSL_SUCCESS;
}
//****************************************************************************80

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

//****************************************************************************80
//
//  Purpose:
//
//    vanderpol_jac() evaluates the jacobian of the vanderpol_ode().
//
//  Discussion:
//
//    This function defines the jacobian of the van der Pol equation:
//
//      f(u,v) = ( v              )
//               ( -u+mu*v*(1-u^2 )
//
//      df     = ( 0 )
//      dt       ( 0 )
//
//      df     = (   0          1          )
//      dy       ( -1-2*mu*v*u  mu*(1-u^2) )
//
//  Licensing:
//
//    This code is distributed under the MIT license. 
//
//  Modified:
//
//    10 November 2023
//
//  Author:
//
//    Taken from the GSL reference manual.
//    This version by John Burkardt.
//
//  Input:
//
//    double T: the current time.
//
//    const double Y[]: the current solution value.
//
//    void *PARAMS: user-specified parameters.
//
//  Output:
//
//    int JAC: returns GSL_SUCCESS indicating a successful computation.
//
//    double *DFDY[]: the Jacobian matrix.
//
//    double DFDT[]: the time derivative of F.
//
{
//
//  Declare t to avoid unused parameter warning.
//
  (void) (t);
  double mu = *(double *) params;
//
//  Allocate a 2x2 matrix for dfdy.
//
  gsl_matrix_view dfdy_mat = gsl_matrix_view_array ( dfdy, 2, 2 );
//
//  Let "m" be a pointer to this matrix.
//
  gsl_matrix * m = &dfdy_mat.matrix;
//
//  Set dFdY
//
  gsl_matrix_set ( m, 0, 0, 0.0 );
  gsl_matrix_set ( m, 0, 1, 1.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) );
//
//  Set dFdT
//
  dfdt[0] = 0.0;
  dfdt[1] = 0.0;

  return GSL_SUCCESS;
}
//****************************************************************************80

void timestamp ( )

//****************************************************************************80
//
//  Purpose:
//
//    timestamp() prints the current YMDHMS date as a time stamp.
//
//  Example:
//
//    31 May 2001 09:45:54 AM
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    08 July 2009
//
//  Author:
//
//    John Burkardt
//
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct std::tm *tm_ptr;
  std::time_t now;

  now = std::time ( NULL );
  tm_ptr = std::localtime ( &now );

  std::strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr );

  std::cout << time_buffer << "\n";

  return;
# undef TIME_SIZE
}