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

using namespace std;

# include "trapezoidal.hpp"

int main ( );
void predator_prey_test ( double tspan[2], double y0[2], int n );
void predator_prey_dydt ( double t, double y[], double dydt[] );
void stiff_test ( double tspan[2], double y0[1], int n );
void stiff_dydt ( double t, double y[], double dydt[] );
double *stiff_exact ( int n, double t[] );
void plot2 ( int n1, double t1[], double y1[], int n2, double t2[], 
  double y2[], string header, string title );
double *r8vec_linspace_new ( int n, double a_first, double a_last );
void timestamp ( );

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

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    trapezoidal_test() tests trapezoidal().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    17 November 2023
//
//  Author:
//
//    John Burkardt
//
{
  int n;
  double tspan[2];
  double y0[2];

  timestamp ( );
  cout << "\n";
  cout << "trapezoidal_test():\n";
  cout << "  C++ version\n";
  cout << "  Test trapezoidal() on several ODE's.\n";

  tspan[0] = 0.0;
  tspan[1] = 5.0;
  y0[0] = 5000.0;
  y0[1] = 100.0;
  n = 1000;
  predator_prey_test ( tspan, y0, n );

  tspan[0] = 0.0;
  tspan[1] = 1.0;
  y0[0] = 0.0;
  n = 27;
  stiff_test ( tspan, y0, n );
//
//  Terminate.
//
  cout << "\n";
  cout << "trapezoidal_test():\n";
  cout << "  Normal end of execution.\n";
  cout << "\n";
  timestamp ( );

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

void predator_prey_test ( double tspan[2], double y0[2], int n )

//****************************************************************************80
//
//  Purpose:
//
//    predator_prey_test() solves predator-prey ODE using trapezoidal().
//
//  Discussion:
//
//    The physical system under consideration is a pair of animal populations.
//
//    The PREY reproduce rapidly for each animal alive at the beginning of the
//    year, two more will be born by the end of the year.  The prey do not have
//    a natural death rate instead, they only die by being eaten by the predator.
//    Every prey animal has 1 chance in 1000 of being eaten in a given year by
//    a given predator.
//
//    The PREDATORS only die of starvation, but this happens very quickly.
//    If unfed, a predator will tend to starve in about 1/10 of a year.
//    On the other hand, the predator reproduction rate is dependent on
//    eating prey, and the chances of this depend on the number of available prey.
//
//    The resulting differential equations can be written:
//
//      PREY(0) = 5000         
//      PRED(0) =  100
//
//      d PREY / dT =    2 * PREY(T) - 0.001 * PREY(T) * PRED(T)
//      d PRED / dT = - 10 * PRED(T) + 0.002 * PREY(T) * PRED(T)
//
//    Here, the initial values (5000,100) are a somewhat arbitrary starting point.
//
//    The pair of ordinary differential equations that result have an interesting
//    behavior.  For certain choices of the interaction coefficients (such as
//    those given here), the populations of predator and prey will tend to 
//    a periodic oscillation.  The two populations will be out of phase the number
//    of prey will rise, then after a delay, the predators will rise as the prey
//    begins to fall, causing the predator population to crash again.
//
//    There is a conserved quantity, which here would be:
//      E(r,f) = 0.002 r + 0.001 f - 10 ln(r) - 2 ln(f)
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    10 November 2023
//
//  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 tspan[2]: the initial and final times.
//    A reasonable value might be [ 0, 5 ].
//
//    double y0[2]: the initial number of prey and predators.
//    A reasonable value might be [ 5000, 100 ].
//
//    int n: the number of time steps.
//
{
  string command_filename;
  ofstream command;
  string data_filename;
  ofstream data;
  string header = "predator_prey";
  int i;
  const int m = 2;
  double *t;
  double *y;

  cout << "\n";
  cout << "predator_prey_test():\n";
  cout << "  A pair of ordinary differential equations for a population\n";
  cout << "  of predators and prey are solved using trapezoidal().\n";

  t = new double[n+1];
  y = new double[(n+1)*m];

  trapezoidal ( predator_prey_dydt, tspan, y0, n, m, t, y );
//
//  Create the data file.
//
  data_filename = header + "_data.txt";
  data.open ( data_filename.c_str ( ) );

  for ( i = 0; i < n; i++ )
  {
    data << "  " << t[i]
         << "  " << y[0+i*m]
         << "  " << y[1+i*m] << "\n";
  }

  data.close ( );

  cout << "\n";
  cout << "  predator_prey_test: data stored in '" << data_filename << "'.\n";
//
//  Create the command file.
//
  command_filename = header + "_commands.txt";

  command.open ( command_filename.c_str ( ) );

  command << "# " << command_filename << "\n";
  command << "#\n";
  command << "# Usage:\n";
  command << "#  gnuplot < " << command_filename << "\n";
  command << "#\n";
  command << "set term png\n";
  command << "set output '" << header << ".png'\n";
  command << "set xlabel '<-- PREDATOR -->'\n";
  command << "set ylabel '<-- PREY -->'\n";
  command << "set title 'Predator prey: trapezoidal'\n";
  command << "set grid\n";
  command << "set style data lines\n";
  command << "plot '" << data_filename << "' using 2:3 with lines lw 3\n";
  command << "quit\n";

  command.close ( );

  cout << "  predator_prey_test: plot commands stored in '" << command_filename << "'.\n";

  delete [] t;
  delete [] y;

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

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

//****************************************************************************80
//
//  Purpose:
//
//    predator_prey_dydt() evaluates the right hand side of predator_prey_ode().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    10 November 2023
//
//  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 y: 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 stiff_test ( double tspan[2], double y0[1], int n )

//****************************************************************************80
//
//  Purpose:
//
//    stiff_test() uses trapezoidal() on stiff_ode().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    10 November 2023
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double TSPAN(2): the first and last times.
//
//    double Y0: the initial condition.
//
//    int N: the number of steps to take.
//
{
  const int m = 1;
  const int n2 = 101;
  double *t1;
  double *t2;
  double *y1;
  double *y2;

  cout << "\n";
  cout << "stiff_test():\n";
  cout << "  Solve stiff_ode() using trapezoidal().\n";

  t1 = new double[n+1];
  y1 = new double[n+1];
  trapezoidal ( stiff_dydt, tspan, y0, n, m, t1, y1 );

  t2 = r8vec_linspace_new ( n2, tspan[0], tspan[1] );
  y2 = stiff_exact ( n2, t2 );

  plot2 ( n+1, t1, y1, n2, t2, y2, "stiff", "Stiff ODE: trapezoidal method" );

  delete [] t1;
  delete [] t2;
  delete [] y1;
  delete [] y2;

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

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

//****************************************************************************80
//
//  stiff_dydt() evaluates the right hand side of stiff_ode().
//
//  Discussion:
//
//    y' = 50 * ( cos(t) - y )
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    10 November 2023
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    double T, Y[1]: the time and solution value.
//
//  Output:
//
//    double DYDT[1]: the derivative value.
//
{
  dydt[0] = 50.0 * ( cos ( t ) - y[0] );

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

double *stiff_exact ( int n, double t[] )

//****************************************************************************80
//
//  Purpose:
//
//    stiff_exact() evaluates the exact solution of stiff_ode().
//
//  Discussion:
//
//    y' = 50 * ( cos(t) - y )
//    y(0) = 0
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 April 2020
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N: the number of values.
//
//    double T[N]: the evaluation times.
//
//  Output:
//
//    double STIFF_EXACT[N]: the exact solution values.
//
{
  int i;
  double *y;

  y = new double[n];

  for ( i = 0; i < n; i++ )
  {
    y[i] = 50.0 * ( sin ( t[i] ) + 50.0 * cos ( t[i] )
      - 50.0 * exp ( - 50.0 * t[i] ) ) / 2501.0;
  }

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

void plot2 ( int n1, double t1[], double y1[], int n2, double t2[], 
  double y2[], string header, string title )

//****************************************************************************80
//
//  Purpose:
//
//    plot2() plots two curves together.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 April 2020
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    int N1: the size of the first data set.
//
//    double T1(N1), Y1(N1), the first dataset.
//
//    int N2: the size of the second data set.
//
//    double T2(N2), Y2(N2), the secod dataset.
//
//    string HEADER: an identifier for the data.
//
//    string TITLE: a title to appear in the plot.
//
{
  string command_filename;
  ofstream command;
  string data1_filename;
  string data2_filename;
  ofstream data;
  int i;
//
//  Create the data files.
//
  data1_filename = header + "_data1.txt";
  data.open ( data1_filename.c_str ( ) );
  for ( i = 0; i < n1; i++ )
  {
    data << "  " << t1[i] << "  " << y1[i] << "\n";
  }
  data.close ( );

  data2_filename = header + "_data2.txt";
  data.open ( data2_filename.c_str ( ) );
  for ( i = 0; i < n2; i++ )
  {
    data << "  " << t2[i] << "  " << y2[i] << "\n";
  }
  data.close ( );

  cout << "\n";
  cout << "  plot2: data stored in '" << data1_filename
       << "' and '" << data2_filename << "'\n";
//
//  Create the command file.
//
  command_filename = header + "_commands.txt";
  command.open ( command_filename.c_str ( ) );

  command << "# " << command_filename << "\n";
  command << "#\n";
  command << "# Usage:\n";
  command << "#  gnuplot < " << command_filename << "\n";
  command << "#\n";
  command << "set term png\n";
  command << "set output '" << header << ".png'\n";
  command << "set xlabel '<-- T -->'\n";
  command << "set ylabel '<-- Y(T) -->'\n";
  command << "set title '" << title << "'\n";
  command << "set grid\n";
  command << "set style data lines\n";
  command << "plot '" << data1_filename << "' using 1:2 with lines lw 3 lt rgb 'red',\\\n";
  command << "     '" << data2_filename << "' using 1:2 with lines lw 3 lt rgb 'blue'\n";
  command << "quit\n";

  command.close ( );

  cout << "  plot2: plot commands stored in '" << command_filename << "'.\n";

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

double *r8vec_linspace_new ( int n, double a_first, double a_last )

//****************************************************************************80
//
//  Purpose:
//
//    r8vec_linspace_new() creates a vector of linearly spaced values.
//
//  Discussion:
//
//    An R8VEC is a vector of R8's.
//
//    4 points evenly spaced between 0 and 12 will yield 0, 4, 8, 12.
//
//    In other words, the interval is divided into N-1 even subintervals,
//    and the endpoints of intervals are used as the points.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 March 2011
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, int N, the number of entries in the vector.
//
//    Input, double A_FIRST, A_LAST, the first and last entries.
//
//    Output, double R8VEC_LINSPACE_NEW[N], a vector of linearly spaced data.
//
{
  double *a;
  int i;

  a = new double[n];

  if ( n == 1 )
  {
    a[0] = ( a_first + a_last ) / 2.0;
  }
  else
  {
    for ( i = 0; i < n; i++ )
    {
      a[i] = ( ( double ) ( n - 1 - i ) * a_first 
             + ( double ) (         i ) * a_last ) 
             / ( double ) ( n - 1     );
    }
  }
  return a;
}
//****************************************************************************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
}