# include <float.h>
# include <math.h>
# include <stdio.h>
# include <stdlib.h>
# include <time.h>

# include "bisection.h"

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

/******************************************************************************/

int main ( )

/******************************************************************************/
/*
  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 ( );
  printf ( "\n" );
  printf ( "bisection_test():\n" );
  printf ( "  C version\n" );
  printf ( "  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.
*/
  printf ( "\n" );
  printf ( "bisection_test():\n" );
  printf ( "  Normal end of execution.\n" );
  timestamp ( );

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

void bisection_example ( double a, double b, double f ( double x ), 
  char *f_string )

/******************************************************************************/
/*
  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);

  printf ( "\n" );
  printf ( "  Function = '%s'\n", f_string );
  printf ( "  a = %g, f(a) = %g\n", a2, fa );
  printf ( "  b = %g, f(b) = %g\n", b2, fb );
  printf ( "  Interval tolerance = %g\n", tol );
  printf ( "  Number of bisections = %d\n", it );
  printf ( "  x = %g, f(x) = %g\n", x, fx );

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

double fcos ( double x )

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

double fpoly ( double x )

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

double kepler ( double x )

/******************************************************************************/
/*
  Purpose:

    kepler() evaluates a version of Kepler's equation.

  Discussion:

    Kepler's equation relates the mean anomaly M, the eccentric anomaly E,
    andthe 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;
}
/******************************************************************************/

void timestamp ( )

/******************************************************************************/
/*
  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:

    24 September 2003

  Author:

    John Burkardt
*/
{
# define TIME_SIZE 40

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

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

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

  fprintf ( stdout, "%s\n", time_buffer );

  return;
# undef TIME_SIZE
}