# include <cmath>
# include <cstdlib>
# include <cstring>
# include <ctime>
# include <iostream>

using namespace std;

# include "midpoint_adaptive.hpp"

int main ( );
void lotka_dydt ( double t, double y[], double dy[] );
void lotka_test ( );
void timestamp ( );

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

int main ( )

//****************************************************************************80
//
// Purpose:
//
//    midpoint_adaptive_test() tests midpoint_adaptive().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    10 July 2024
//
//  Author:
//
//    John Burkardt
//
{
  timestamp ( );
  cout << "\n";
  cout << "midpoint_adaptive_test():\n";
  cout << "  C++ version\n";
  cout << "  midpoint_adaptive() solves ODE's using an adaptive\n";
  cout << "  version of the midpoint method.\n";

  lotka_test ( );
//
//  Terminate.
//
  cout << "\n";
  cout << "midpoint_adaptive_test():\n";
  cout << "  Normal end of execution.\n";
  cout << "\n";
  timestamp ( );

  return ( 0 );
}
//****************************************************************************80

void lotka_test ( )

//****************************************************************************80
//
//  Purpose:
//
//    lotka_test() solves lotka_ode() using midpoint_adaptive().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    15 July 2024
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    George Lindfield, John Penny,
//    Numerical Methods Using MATLAB,
//    Second Edition,
//    Prentice Hall, 1999,
//    ISBN: 0-13-012641-1,
//    LC: QA297.P45.
//
{
  double abstol;
  string label;
  int m;
  int n_fsolve;
  int n_rejected;
  int nmax;
  int nstep;
  double reltol;
  double *t;
  double t0;
  double tau0;
  double tmax;
  double *y;
  double y0[2];

  cout << "\n";
  cout << "lotka_test():\n";
  cout << "  A pair of ordinary differential equations for a population\n";
  cout << "  of predators and prey are solved using midpoint().\n";
  cout << "\n";
  cout << "  Steady state at r = 5000, f = 2000.\n";
  label = "lotka";

  t0 = 0.0;
  tmax = 20.0;
  m = 2;
  y0[0] = 5000.0;
  y0[1] = 100.0;
  tau0 = 2.5E-2;
  reltol = 1.0E-3;
  abstol = 1.0E-3;
  nmax = 400;
  nstep = 0;
  t = new double[nmax+1];
  y = new double[(nmax+1)*m];

  mad_compute ( lotka_dydt, t0, tmax, m, y0, tau0, reltol, abstol, nmax, 
    nstep, n_rejected, n_fsolve, t, y );

  mad_phase_plot ( nstep, m, t, y, label );

  mad_solution_plot ( nstep, m, t, y, label );

  mad_timestep_plot ( nstep, t, label );

  mad_stats ( nstep, n_rejected, n_fsolve, t );
//
//  Free memory.
//
  delete [] y;
  delete [] t;

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

void lotka_dydt ( double t, double y[], double dydt[] )

//****************************************************************************80
//
//  Purpose:
//
//    lotka_dydt() evaluates the right hand side of lotka_ode().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    10 July 2024
//
//  Author:
//
//    John Burkardt
//
//  Reference:
//
//    George Lindfield, John Penny,
//    Numerical Methods Using MATLAB,
//    Second Edition,
//    Prentice Hall, 1999,
//    ISBN: 0-13-012641-1,
//    LC: QA297.P45.
//
//  Input:
//
//    double t: the current time.
//
//    double y[2]: the current solution variables, rabbits and foxes.
//
//  Output:
//
//    double dydt[2]: the right hand side of the 2 ODE's.
//
{
  dydt[0] =    2.0 * y[0] - 0.001 * y[0] * y[1];
  dydt[1] = - 10.0 * y[1] + 0.002 * y[0] * y[1];

  return;
}
//****************************************************************************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:
//
//    19 March 2018
//
//  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
}