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

using namespace std;

int main ( int argc, char *argv[] );
void dtable_close_write ( ofstream &output );
void dtable_data_write ( ofstream &output, int m, int n, double table[] );
void dtable_header_write ( string output_filename, ofstream &output, int m, 
  int n );
void dtable_write ( string output_filename, int m, int n, double table[], 
  bool header );
double r8_epsilon ( );
int s_len_trim ( char *s );
int s_to_i4 ( char *s, int *last, bool *error );
void timestamp ( );
char *timestring ( );

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

int main ( int argc, char *argv[] )

//****************************************************************************80
//
//  Purpose:
//
//    FD_PREDATOR_PREY solves a pair of predator-prey ODE's.
//
//  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:
//
//      PRED(0) = 5000
//      PREY(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.
//
//    In this program, the pair of ODE's is solved with a simple finite difference
//    approximation using a fixed step size.  In general, this is NOT an efficient
//    or reliable way of solving differential equations.  However, this program is
//    intended to illustrate the ideas of finite difference approximation.
//
//    In particular, if we choose a fixed time step size DT, then a derivative
//    such as dPREY/dT is approximated by:
//
//      d PREY / dT = approximately ( PREY(T+DT) - PREY(T) ) / DT
//
//    which means that the first differential equation can be written as
//
//      PREY(T+DT) = PREY(T) + DT * ( 2 * PREY(T) - 0.001 * PREY(T) * PRED(T) ).
//
//    We can rewrite the equation for PRED as well.  Then, since we know the
//    values of PRED and PREY at time 0, we can use these finite difference
//    equations to estimate the values of PRED and PREY at time DT.  These values
//    can be used to get estimates at time 2*DT, and so on.  To get from time
//    T_START = 0 to time T_STOP = 5, we simply take STEP_NUM steps each of size
//    DT = ( T_STOP - T_START ) / STEP_NUM.
//
//    Because finite differences are only an approximation to derivatives, this
//    process only produces estimates of the solution.  And these estimates tend
//    to become more inaccurate for large values of time.  Usually, we can reduce
//    this error by decreasing DT and taking more, smaller time steps.
//
//    In this example, for instance, taking just 100 steps gives nonsensical
//    answers.  Using STEP_NUM = 1000 gives an approximate solution that seems
//    to have the right kind of oscillatory behavior, except that the amplitude
//    of the waves increases with each repetition.  Using 10000 steps, the
//    approximation begins to become accurate enough that we can see that the
//    waves seem to have a fixed period and amplitude.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license.
//
//  Modified:
//
//    21 March 2009
//
//  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.
//
//  Parameters:
//
//    Input, int STEP_NUM, the number of steps.
//
{
  double dt;
  bool header;
  int i;
  bool ierror;
  int length;
  string output_filename;
  int step_num;
  double *tpp;
  double t_start;
  double t_stop;

  timestamp ( );

  cout << "\n";
  cout << "FD_PREDATOR_PREY\n";
  cout << "  C++ version\n";
  cout << "\n";
  cout << "  A finite difference approximate solution of a pair\n";
  cout << "  of ordinary differential equations for a population\n";
  cout << "  of predators and prey.\n";
  cout << "\n";
  cout << "  The exact solution shows wave behavior, with a fixed\n";
  cout << "  period and amplitude.  The finite difference approximation\n";
  cout << "  can provide a good estimate for this behavior if the stepsize\n";
  cout << "  DT is small enough.\n";
//
//  STEP_NUM is an input argument or else read from the user interactively.
//
  if ( 1 < argc )
  {
    step_num = s_to_i4 ( argv[1], &length, &ierror );
  }
  else
  {
    cout << "\n";
    cout << "FD_PREDATOR_PREY:\n";
    cout << "  Please enter the number of time steps:\n";

    cin >> step_num;
  }

  t_start = 0.0;
  t_stop =  5.0;
  dt = ( t_stop - t_start ) / ( double ) ( step_num );
//
//  TPP(1:3,1:STEP_NUM+1) contains TIME, PREY, and PREDATOR values for each step.
//
  tpp = new double[3*(step_num+1)];

  tpp[0+0*3] = t_start;
  tpp[1+0*3] = 5000.0;
  tpp[2+0*3] =  100.0;

  for ( i = 0; i < step_num; i++ )
  {
    tpp[0+(i+1)*3] = tpp[0+i*3] + dt;

    tpp[1+(i+1)*3] = tpp[1+i*3] + dt 
      * (    2.0 * tpp[1+i*3] - 0.001 * tpp[1+i*3] * tpp[2+i*3] );

    tpp[2+(i+1)*3] = tpp[2+i*3] + dt 
      * ( - 10.0 * tpp[2+i*3] + 0.002 * tpp[1+i*3] * tpp[2+i*3] );
  }
  output_filename = "solution_";
  output_filename = output_filename + argv[1];
  output_filename = output_filename + ".txt";

  header = false;
  dtable_write ( output_filename, 3, step_num + 1, tpp, header );

  cout << "\n";
  cout << "  Initial time    = " << t_start << "\n";
  cout << "  Final time      = " << t_stop << "\n";
  cout << "  Initial prey    = " << tpp[1+0*3] << "\n";
  cout << "  Initial pred    = " << tpp[2+0*3] << "\n";
  cout << "  Number of steps = " << step_num << "\n";
  cout << "  Solution data written to \"" << output_filename << "\".\n";

  delete [] tpp;

  cout << "\n";
  cout << "FD_PREDATOR_PREY\n";
  cout << "  Normal end of execution.\n";
  cout << "\n";
  timestamp ( );

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

void dtable_close_write ( ofstream &output )

//****************************************************************************80
//
//  Purpose:
//
//    DTABLE_CLOSE_WRITE closes a file to which a DTABLE was to be written.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    27 September 2006
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, ofstream output, the output file stream.
//
{
  output.close ( );

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

void dtable_data_write ( ofstream &output, int m, int n, double table[] )

//****************************************************************************80
//
//  Purpose:
//
//    DTABLE_DATA_WRITE writes data to a DTABLE file.
//
//  Discussion:
//
//    The file should already be open.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    11 December 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, ofstream &OUTPUT, a pointer to the output stream.
//
//    Input, int M, the spatial dimension.
//
//    Input, int N, the number of points.
//
//    Input, double TABLE[M*N], the table data.
//
{
  int i;
  int j;

  for ( j = 0; j < n; j++ )
  {
    for ( i = 0; i < m; i++ )
    {
      output << setw(10) << table[i+j*m] << "  ";
    }
    output << "\n";
  }

  return;
}
//****************************************************************************80
 
void dtable_header_write ( string output_filename, ofstream &output, int m, 
  int n )
 
//****************************************************************************80
//
//  Purpose:
//
//    DTABLE_HEADER_WRITE writes the header of a DTABLE file.
//
//  Discussion:
//
//    The file should already be open.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 February 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string OUTPUT_FILENAME, the output filename.
//
//    Input, ofstream &OUTPUT, the output stream.
//
//    Input, int M, the spatial dimension.
//
//    Input, int N, the number of points.
//
{
  char *s;

  s = timestring ( );

  output << "#  " << output_filename << "\n";
  output << "#  created by dtable_write in table_io.C" << "\n";
  output << "#  at " << s << "\n";
  output << "#\n";
  output << "#  Spatial dimension M = " << m << "\n";
  output << "#  Number of points N =  " << n << "\n";
  output << "#  EPSILON (unit roundoff) = " << r8_epsilon ( ) << "\n";
  output << "#\n";

  delete [] s;

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

void dtable_write ( string output_filename, int m, int n, double table[], 
  bool header )

//****************************************************************************80
//
//  Purpose:
//
//    DTABLE_WRITE writes information to a DTABLE file.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    23 February 2009
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, string OUTPUT_FILENAME, the output filename.
//
//    Input, int M, the spatial dimension.
//
//    Input, int N, the number of points.
//
//    Input, double TABLE[M*N], the table data.
//
//    Input, bool HEADER, is TRUE if the header is to be included.
//
{
  ofstream output;

  output.open ( output_filename.c_str ( ) );

  if ( !output )
  {
    cerr << "\n";
    cerr << "DTABLE_WRITE - Fatal error!\n";
    cerr << "  Could not open the output file.\n";
    return;
  }

  if ( header )
  {
    dtable_header_write ( output_filename, output, m, n );
  }

  dtable_data_write ( output, m, n, table );

  dtable_close_write ( output );

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

double r8_epsilon ( )

//****************************************************************************80
//
//  Purpose:
//
//    R8_EPSILON returns the R8 roundoff unit.
//
//  Discussion:
//
//    The roundoff unit is a number R which is a power of 2 with the 
//    property that, to the precision of the computer's arithmetic,
//      1 < 1 + R
//    but 
//      1 = ( 1 + R / 2 )
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    01 July 2004
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, double R8_EPSILON, the double precision round-off unit.
//
{
  double r;

  r = 1.0;

  while ( 1.0 < ( double ) ( 1.0 + r )  )
  {
    r = r / 2.0;
  }

  return ( 2.0 * r );
}
//****************************************************************************80

int s_len_trim ( char *s )

//****************************************************************************80
//
//  Purpose:
//
//    S_LEN_TRIM returns the length of a string to the last nonblank.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    26 April 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char *S, a pointer to a string.
//
//    Output, int S_LEN_TRIM, the length of the string to the last nonblank.
//    If S_LEN_TRIM is 0, then the string is entirely blank.
//
{
  int n;
  char *t;

  n = strlen ( s );
  t = s + strlen ( s ) - 1;

  while ( 0 < n ) 
  {
    if ( *t != ' ' )
    {
      return n;
    }
    t--;
    n--;
  }

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

int s_to_i4 ( char *s, int *last, bool *error )

//****************************************************************************80
//
//  Purpose:
//
//    S_TO_I4 reads an I4 from a string.
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    13 June 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Input, char *S, a string to be examined.
//
//    Output, int *LAST, the last character of S used to make IVAL.
//
//    Output, bool *ERROR is TRUE if an error occurred.
//
//    Output, int *S_TO_I4, the integer value read from the string.
//    If the string is blank, then IVAL will be returned 0.
//
{
  char c;
  int i;
  int isgn;
  int istate;
  int ival;

  *error = false;
  istate = 0;
  isgn = 1;
  i = 0;
  ival = 0;

  while ( *s ) 
  {
    c = s[i];
    i = i + 1;
//
//  Haven't read anything.
//
    if ( istate == 0 )
    {
      if ( c == ' ' )
      {
      }
      else if ( c == '-' )
      {
        istate = 1;
        isgn = -1;
      }
      else if ( c == '+' )
      {
        istate = 1;
        isgn = + 1;
      }
      else if ( '0' <= c && c <= '9' )
      {
        istate = 2;
        ival = c - '0';
      }
      else
      {
        *error = true;
        return ival;
      }
    }
//
//  Have read the sign, expecting digits.
//
    else if ( istate == 1 )
    {
      if ( c == ' ' )
      {
      }
      else if ( '0' <= c && c <= '9' )
      {
        istate = 2;
        ival = c - '0';
      }
      else
      {
        *error = true;
        return ival;
      }
    }
//
//  Have read at least one digit, expecting more.
//
    else if ( istate == 2 )
    {
      if ( '0' <= c && c <= '9' )
      {
        ival = 10 * (ival) + c - '0';
      }
      else
      {
        ival = isgn * ival;
        *last = i - 1;
        return ival;
      }

    }
  }
//
//  If we read all the characters in the string, see if we're OK.
//
  if ( istate == 2 )
  {
    ival = isgn * ival;
    *last = s_len_trim ( s );
  }
  else
  {
    *error = true;
    *last = 0;
  }

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

void timestamp ( )

//****************************************************************************80
//
//  Purpose:
//
//    TIMESTAMP prints the current YMDHMS date as a time stamp.
//
//  Example:
//
//    May 31 2001 09:45:54 AM
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    02 October 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    None
//
{
# define TIME_SIZE 40

  static char time_buffer[TIME_SIZE];
  const struct tm *tm;
  size_t len;
  time_t now;

  now = time ( NULL );
  tm = localtime ( &now );

  len = strftime ( time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  cout << time_buffer << "\n";

  return;
# undef TIME_SIZE
}
//****************************************************************************80

char *timestring ( )

//****************************************************************************80
//
//  Purpose:
//
//    TIMESTRING returns the current YMDHMS date as a string.
//
//  Example:
//
//    May 31 2001 09:45:54 AM
//
//  Licensing:
//
//    This code is distributed under the GNU LGPL license. 
//
//  Modified:
//
//    13 June 2003
//
//  Author:
//
//    John Burkardt
//
//  Parameters:
//
//    Output, char *TIMESTRING, a string containing the current YMDHMS date.
//
{
# define TIME_SIZE 40

  const struct tm *tm;
  size_t len;
  time_t now;
  char *s;

  now = time ( NULL );
  tm = localtime ( &now );

  s = new char[TIME_SIZE];

  len = strftime ( s, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm );

  return s;
# undef TIME_SIZE
}