# include <cmath>
# include <complex>
# include <ctime>
# include <iomanip>
# include <iostream>

using namespace std;

# include "zero_muller.hpp"

int main ( );
void test01 ( );
void test02 ( );
void test03 ( );
complex <double> func01 ( complex <double> x );
complex <double> func02 ( complex <double> x );
complex <double> func03 ( complex <double> x );
void timestamp ( );

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

int main ( )

//****************************************************************************80
//
//  Purpose:
//
//    zero_muller_test() tests zero_muller().
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 March 2024
//
//  Author:
//
//    John Burkardt
//
{
  timestamp ( );
  cout << "\n";
  cout << "zero_muller_test():\n";
  cout << "  C++ version\n";
  cout << "  Test zero_muller(), which uses Muller's method,\n";
  cout << "  with complex arithmetic, to solve a nonlinear equation.\n";

  test01 ( );
  test02 ( );
  test03 ( );
//
//  Terminate.
//
  cout << "\n";
  cout << "zero_muller_test():\n";
  cout << "  Normal end of execution.\n";
  cout << "\n";
  timestamp ( );

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

void test01 ( )

//****************************************************************************80
//
//  Purpose:
//
//    test01() tests zero_muller() on F(X) = X*X+9.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 March 2024
//
//  Author:
//
//    John Burkardt
//
{
  double fatol;
  complex <double> fxnew;
  int itmax;
  complex <double> x1;
  complex <double> x2;
  complex <double> x3;
  complex <double> xnew;
  double xatol;
  double xrtol;

  cout << "\n";
  cout << "test01():\n";
  cout << "  Demonstrate zero_muller() on F(X) = X*X+9.\n";

  fatol = 1.0E-05;
  itmax = 10;
  x1 = complex <double> ( 1.0, 0.0 );
  x2 = complex <double> ( 0.0, 1.0 );
  x3 = complex <double> ( 0.5, 0.5 );
  xatol = 1.0E-05;
  xrtol = 1.0E-05;

  zero_muller ( func01, fatol, itmax, x1, x2, x3, xatol, xrtol, 
    xnew, fxnew );

  cout << "\n";
  cout << "  X = " << real ( xnew ) << "  " << imag ( xnew ) << "\n";
  cout << "\n";
  cout << "  with function value F(X):\n";
  cout << "\n";
  cout << "    FX   = " << real ( fxnew ) << "  " << imag ( fxnew ) << "\n";
  cout << "  ||FX|| = " << abs ( fxnew ) << "\n";

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

void test02 ( )

//****************************************************************************80
//
//  Purpose:
//
//    test02() tests zero_muller() on F(X) = (X*X+4) * (X-10) * (X+20)
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 March 2024
//
//  Author:
//
//    John Burkardt
//
{
  double fatol;
  complex <double> fxnew;
  int itmax;
  complex <double> x1;
  complex <double> x2;
  complex <double> x3;
  complex <double> xnew;
  double xatol;
  double xrtol;

  cout << "\n";
  cout << "test02():\n";
  cout << "  Demonstrate zero_muller() on F(X) = (X*X+4)*(X-10)*(X+20).\n";

  fatol = 1.0E-05;
  itmax = 10;
  x1 = complex <double> ( 1.0, 0.0 );
  x2 = complex <double> ( 0.0, 1.0 );
  x3 = complex <double> ( 0.5, 0.5 );
  xatol = 1.0E-05;
  xrtol = 1.0E-05;

  zero_muller ( func01, fatol, itmax, x1, x2, x3, xatol, xrtol, 
    xnew, fxnew );

  cout << "\n";
  cout << "  X = " << real ( xnew ) << "  " << imag ( xnew ) << "\n";
  cout << "\n";
  cout << "  with function value F(X):\n";
  cout << "\n";
  cout << "    FX   = " << real ( fxnew ) << "  " << imag ( fxnew ) << "\n";
  cout << "  ||FX|| = " << abs ( fxnew ) << "\n";

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

void test03 ( )

//****************************************************************************80
//
//  Purpose:
//
//    test03() tests zero_muller() on Zhelyazkov's function
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 March 2024
//
//  Author:
//
//    John Burkardt
//
{
  double fatol;
  complex <double> fxnew;
  int itmax;
  int test;
  complex <double> x1;
  complex <double> x2;
  complex <double> x3;
  complex <double> xnew;
  double xatol;
  double xrtol;

  cout << "\n";
  cout << "test03():\n";
  cout << "  Demonstrate zero_muller() on Zhelyazkov''s function.\n";

  fatol = 1.0E-07;
  itmax = 10;

  for ( test = 1; test <= 2; test++ )
  {
//
//  First set of starting points.
//  Result is X = ( 1.5705798926, 0.0 )
//
    if ( test == 1 )
    {
      x1 = complex <double> ( 1.0, 0.0 );
      x2 = complex <double> ( 0.0, 1.0 );
      x3 = complex <double> ( 0.5, 0.5 );
    }
//
//  Second set of starting points.
//  Result is X = ( -0.5802520567, 0.0 ).
//
    else if ( test == 2 )
    {
      x1 = complex <double> (  0.0, 1.0 );
      x2 = complex <double> (  1.0, 2.0 );
      x3 = complex <double> ( -1.0, 2.0 );
    }

    xatol = 1.0E-07;
    xrtol = 1.0E-07;

    zero_muller ( func01, fatol, itmax, x1, x2, x3, xatol, xrtol, 
      xnew, fxnew );

    cout << "\n";
    cout << "  X = " << real ( xnew ) << "  " << imag ( xnew ) << "\n";
    cout << "\n";
    cout << "  with function value F(X):\n";
    cout << "\n";
    cout << "    FX   = " << real ( fxnew ) << "  " << imag ( fxnew ) << "\n";
    cout << "  ||FX|| = " << abs ( fxnew ) << "\n";
  }

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

complex <double> func01 ( complex <double> x )

//****************************************************************************80
//
//  Purpose:
//
//    func01() evaluates F(X) = X*X+9.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 March 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    complex <double> X, the point at which the function is to
//    be evaluated.
//
//  Output:
//
//    complex <double> FX, the function value at X.
//
{
  complex <double> fx;

  fx = x * x + 9.0;
 
  return fx;
}
//****************************************************************************80

complex <double> func02 ( complex <double> x )

//****************************************************************************80
//
//  Purpose:
//
//    func02() evaluates F(X) = (X*X+4)*(X-1)*(X+2).
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    28 March 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    complex <double> X, the point at which the function is to
//    be evaluated.
//
//  Output:
//
//    complex <double> FX, the function value at X.
//
{
  complex <double> fx;

  fx = ( x * x + 4.0 ) * ( x - 10.0 ) * ( x + 20.0 );
 
  return fx;
}
//****************************************************************************80

complex <double> func03 ( complex <double> z )

//****************************************************************************80
//
//  Purpose:
//
//    func03() evaluates Zhelyazkov's function.
//
//  Licensing:
//
//    This code is distributed under the MIT license.
//
//  Modified:
//
//    29 March 2024
//
//  Author:
//
//    John Burkardt
//
//  Input:
//
//    complex <double> Z, the point at which the function is to
//    be evaluated.
//
//  Output:
//
//    complex <double> FZ, the function value at Z.
//
{
  complex <double> a;
  complex <double> b;
  complex <double> eps = complex <double> ( 0.4, 0.0 );
  complex <double> eta = complex <double> ( 0.64, 0.0 );
  complex <double> fz;
  double me = 0.384;
  double mo = 0.5;
  complex <double> of;
  complex <double> ok;
  complex <double> one = complex <double> ( 1.0, 0.0 );
  double x = 0.5;

  ok = z - me / sqrt ( eta );
  of = z - mo;

  a = of * of + ( ok * ok ) * eta * tanh ( x );

  b = ( of - ok * eta ) / ( of - ok * eta * eta );

  fz = of * of - one + ( eta * ok * ok - one ) * 
    tanh ( x ) - x * x * eps * eps * a * b;

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