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

using namespace std;

# include "bisection.hpp"

int main ( );
void bisection_example ( double a, double b, double f ( double x ), 
  string f_string );
double fcos ( double x );
double fpoly ( double x );
double kepler ( double x );
void timestamp ( );

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

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    bisection_test() tests bisection().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    05 December 2023
//
//  Author:
//
//    John Burkardt
//
{
  double a;
  double b;

  timestamp ( );
  cout << "\n";
  cout << "bisection_test():\n";
  cout << "  C++ version\n";
  cout << "  Test bisection()\n";

  a = 0.0;
  b = 8.0;
  bisection_example ( a, b, fpoly, "x^2 - 2*x - 15" );

  a = 0.0;
  b = 1.0;
  bisection_example ( a, b, fcos, "cos(x) - x" );

  a = 0.0;
  b = 10.0;
  bisection_example ( a, b, kepler, "Kepler function"  );
//
//  Terminate.
//
  cout << "\n";
  cout << "bisection_test():\n";
  cout << "  Normal end of execution.\n";
  timestamp ( );

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

void bisection_example ( double a, double b, double f ( double x ), 
  string f_string )

//****************************************************************************80
//
//  Purpose:
//
//    bisection_example() applies bisection() to a particular example.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    04 December 2023
//
//  Author:
//
//    John Burkardt
//
{
  double a2;
  double b2;
  double fa;
  double fb;
  double fx;
  int it;
  double tol;
  double x;

  a2 = a;
  b2 = b;

  tol = 10.0 * DBL_EPSILON * ( b2 - a2 );
  bisection ( a2, b2, tol, f, it );
  x = ( a2 + b2 ) / 2.0;
  fa = f(a2);
  fb = f(b2);
  fx = f(x);

  cout << "\n";
  cout << "  Function = " << f_string << "\n";
  cout << "  a = " << a2 << ", f(a) = " << fa << "\n";
  cout << "  b = " << b2 << ", f(b) = " << fb << "\n";
  cout << "  Interval tolerance = " << tol << "\n";
  cout << "  Number of bisections = " << it << "\n";
  cout << "  x = " << x << ", f(x) = " << fx << "\n";

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

double fcos ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    fcos() evaluates the function cos(x)-x.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    05 December 2023
//
//  Input:
//
//    double x, the argument.
//
//  Output:
//
//    double fcos: the function value.
//
{
  double value;

  value = cos ( x ) - x;

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

double fpoly ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    fpoly() evaluates a polynomial function.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    05 December 2023
//
//  Input:
//
//    double x, the argument.
//
//  Output:
//
//    double fpoly: the function value.
//
{
  double value;

  value = x * x - 2.0 * x - 15.0;

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

double kepler ( double x )

//****************************************************************************80
//
//  Purpose:
//
//    kepler() evaluates a version of Kepler's equation.
//
//  Discussion:
//
//    Kepler's equation relates the mean anomaly M, the eccentric anomaly E,
//    and the eccentricity e of a planetary orbit.
//
//    Typically, e is a fixed feature of the orbit, the value of M is determined
//    by observation, and the value of E is desired.
//
//    Kepler's equation states that:
//      M = E - e sin(E)
//
//    Suppose we have an orbit with e = 2, and we have observed M = 5.  What is
//    the value of E?  The equation becomes:
//      5 = E - 2 sin ( E ).
//
//    To solve for E, we need to rewrite this as a function:
//      F(E) = 5 - E + 2 sin ( E )
//    and then use a nonlinear equation solver to solve for the value of E
//    such that F(E)=0.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    05 December 2023
//
//  Input:
//
//    double x, the current estimate for the value of E.
//
//  Output:
//
//    double kepler: the Kepler equation residual F(E).
//
{
  double value;

  value = 5.0 - x + 2.0 * sin ( x );

  return value;
}
//****************************************************************************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
}