program main

c*********************************************************************72
c
cc zero_muller_test() tests zero_muller().
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    27 March 2024
c
c  Author:
c
c    John Burkardt
c
      implicit none

      call timestamp ( )
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'zero_muller_test():'
      write ( *, '(a)' ) '  Fortran77 version'
      write ( *, '(a)' ) 
     &  '  Test zero_muller(), which uses Muller''s method,'
      write ( *, '(a)' ) 
     &  '  with complex arithmeic, to solve a nonlinear equation.'

      call test01 ( )
      call test02 ( )
      call test03 ()
c
c  Terminate.
c
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'zero_muller_test():'
      write ( *, '(a)' ) '  Normal end of execution.'
      write ( *, '(a)' ) ''
      call timestamp ( )

      stop 0
      end
      subroutine test01 ( )

c*********************************************************************72
c
cc test01() tests zero_muller() on F(X) = X*X+9.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    19 July 2007
c
c  Author:
c
c    John Burkardt
c
      implicit none

      double precision fatol
      external func01
      double complex fxnew
      integer itmax
      double complex x1
      double complex x2
      double complex x3
      double complex xnew
      double precision xatol
      double precision xrtol

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'test01():'
      write ( *, '(a)' ) '  Demonstrate zero_muller() on F(X) = X*X+9.'

      fatol = 1.0D-05
      itmax = 10
      x1 = dcmplx ( 1.0D+00, 0.0D+00 )
      x2 = dcmplx ( 0.0D+00, 1.0D+00 )
      x3 = dcmplx ( 0.5D+00, 0.5D+00 )
      xatol = 1.0D-05
      xrtol = 1.0D-05

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

      write ( *, '(a)' ) ' '
      write ( *, '(a,f20.10,f20.10)' ) '     X   = ', xnew
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  with function value F(X):'
      write ( *, '(a)' ) ' '
      write ( *, '(a,f20.10,f20.10)' ) '    FX   = ', fxnew
      write ( *, '(a,f20.10)' )        '  ||FX|| = ', abs ( fxnew )

      return
      end
      subroutine test02 ( )

c*********************************************************************72
c
cc test02() tests zero_muller() on F(X) = (X*X+4) * (X-10) * (X+20)
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    19 July 2007
c
c  Author:
c
c    John Burkardt
c
      implicit none

      double precision fatol
      external func02
      double complex fxnew
      integer itmax
      double complex x1
      double complex x2
      double complex x3
      double complex xnew
      double precision xatol
      double precision xrtol

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'test02():'
      write ( *, '(a)' ) '  Demonstrate zero_muller()'
      write ( *, '(a)' ) '  F(X) = (X*X+4)*(X-10)*(X+20).'

      fatol = 1.0D-05
      itmax = 10
      x1 = dcmplx ( 1.0D+00, 0.0D+00 )
      x2 = dcmplx ( 0.0D+00, 1.0D+00 )
      x3 = dcmplx ( 0.5D+00, 0.5D+00 )
      xatol = 1.0D-05
      xrtol = 1.0D-05

      call zero_muller ( func02, fatol, itmax, x1, x2, x3, xatol, 
     &  xrtol, xnew, fxnew )

      write ( *, '(a)' ) ' '
      write ( *, '(a,f20.10,f20.10)' ) '     X   = ', xnew
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) '  with function value F(X):'
      write ( *, '(a)' ) ' '
      write ( *, '(a,f20.10,f20.10)' ) '    FX   = ', fxnew
      write ( *, '(a,f20.10)' )        '  ||FX|| = ', abs ( fxnew )

      return
      end
      subroutine test03 ( )

c*********************************************************************72
c
cc test03() tests zero_muller() on Zhelyazkov's function
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    19 July 2007
c
c  Author:
c
c    John Burkardt
c
      implicit none

      double precision fatol
      external func03
      double complex fxnew
      integer itmax
      integer test
      double complex x1
      double complex x2
      double complex x3
      double complex xnew
      double precision xatol
      double precision xrtol

      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'test03():'
      write ( *, '(a)' ) '  Demonstrate zero_muller()'
      write ( *, '(a)' ) '  f(x) = Zhelyazkov''s function.'

      fatol = 1.0D-07
      itmax = 10

      do test = 1, 2
c
c  First set of starting points.
c  Result is X = ( 1.5705798926, 0.0 )
c
        if ( test == 1 ) then
          x1 = dcmplx ( 1.0D+00, 0.0D+00 )
          x2 = dcmplx ( 0.0D+00, 1.0D+00 )
          x3 = dcmplx ( 0.5D+00, 0.5D+00 )
c
c  Second set of starting points.
c  Result is X = ( -0.5802520567, 0.0 ).
c
        else if ( test == 2 ) then
          x1 = dcmplx (  0.0D+00, 1.0D+00 )
          x2 = dcmplx (  1.0D+00, 2.0D+00 )
          x3 = dcmplx ( -1.0D+00, 2.0D+00 )
        end if

        xatol = 1.0D-07
        xrtol = 1.0D-07

        call zero_muller ( func03, fatol, itmax, x1, x2, x3, xatol, 
     &    xrtol, xnew, fxnew )

        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) ' '
        write ( *, '(a,f20.10,f20.10)' ) '     X   = ', xnew
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) '  with function value F(X):'
        write ( *, '(a)' ) ' '
        write ( *, '(a,f20.10,f20.10)' ) '    FX   = ', fxnew
        write ( *, '(a,f20.10)' )        '  ||FX|| = ', abs ( fxnew )

      end do

      return
      end
      subroutine func01 ( x, fx )

c*********************************************************************72
c
cc func01() evaluates F(X) = X*X+9.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    19 July 2007
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    double complex X, the point at which the function is to
c    be evaluated.
c
c  Output:
c
c    double complex FX, the function value at X.
c
      implicit none

      double complex fx
      double complex x

      fx = x * x + 9.0D+00
 
      return
      end
      subroutine func02 ( x, fx )

c*********************************************************************72
c
cc func02() evaluates F(X) = (X*X+4)*(X-1)*(X+2).
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    19 July 2007
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    double complex X, the point at which the function is to
c    be evaluated.
c
c  Output:
c
c    double complex FX, the function value at X.
c
      implicit none

      double complex fx
      double complex x

      fx = ( x * x + 4.0D+00 ) * ( x - 10.0D+00 ) * ( x + 20.0D+00 )
 
      return
      end
      subroutine func03 ( z, fz )

c*********************************************************************72
c
cc func03() evaluates Zhelyazkov's function.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    19 July 2007
c
c  Author:
c
c    John Burkardt
c
c  Input:
c
c    double complex Z, the point at which the function is to
c    be evaluated.
c
c  Output:
c
c    double complex FZ, the function value at Z.
c
      implicit none

      double complex a
      double complex b
      double complex, parameter :: eps = dcmplx ( 0.4D+00, 0.0D+00 )
      double complex, parameter :: eta = dcmplx ( 0.64D+00, 0.0D+00 )
      double complex fz
      double precision, parameter :: me = 0.384D+00
      double precision, parameter :: mo = 0.5D+00
      double complex of
      double complex ok
      double complex, parameter :: one = dcmplx ( 1.0D+00, 0.0D+00 )
      double precision, parameter :: x = 0.5D+00
      double complex z

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

      a = of * of + ( ok * ok ) * eta * dtanh ( x )

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

      fz = of * of - one + ( eta * ok * ok - one ) * 
     &  dtanh ( x ) - x * x * eps * eps * a * b

      return
      end
      subroutine timestamp ( )

c*********************************************************************72
c
cc timestamp() prints the YMDHMS date as a timestamp.
c
c  Licensing:
c
c    This code is distributed under the MIT license.
c
c  Modified:
c
c    12 June 2014
c
c  Author:
c
c    John Burkardt
c
      implicit none

      character * ( 8 ) ampm
      integer d
      character * ( 8 ) date
      integer h
      integer m
      integer mm
      character * ( 9 ) month(12)
      integer n
      integer s
      character * ( 10 ) time
      integer y

      save month

      data month /
     &  'January  ', 'February ', 'March    ', 'April    ', 
     &  'May      ', 'June     ', 'July     ', 'August   ', 
     &  'September', 'October  ', 'November ', 'December ' /

      call date_and_time ( date, time )

      read ( date, '(i4,i2,i2)' ) y, m, d
      read ( time, '(i2,i2,i2,1x,i3)' ) h, n, s, mm

      if ( h .lt. 12 ) then
        ampm = 'AM'
      else if ( h .eq. 12 ) then
        if ( n .eq. 0 .and. s .eq. 0 ) then
          ampm = 'Noon'
        else
          ampm = 'PM'
        end if
      else
        h = h - 12
        if ( h .lt. 12 ) then
          ampm = 'PM'
        else if ( h .eq. 12 ) then
          if ( n .eq. 0 .and. s .eq. 0 ) then
            ampm = 'Midnight'
          else
            ampm = 'AM'
          end if
        end if
      end if

      write ( *, 
     &  '(i2,1x,a,1x,i4,2x,i2,a1,i2.2,a1,i2.2,a1,i3.3,1x,a)' ) 
     &  d, trim ( month(m) ), y, h, ':', n, ':', s, '.', mm, 
     &  trim ( ampm )

      return
      end