program main

!*****************************************************************************80
!
!! P01_OPT01 generates points for TEST_CON problem 1, option 1.
!
!  Discussion:
!
!    We are lucky for this problem because we know an exact solution
!    in which X(2) is the independent variable.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    06 October 2008
!
!  Author:
!
!    John Burkardt
!
  implicit none

  character ( len = 80 ) file_t_name
  character ( len = 80 ) file_x_name
  integer   ( kind = 4 ) file_t_unit
  integer   ( kind = 4 ) file_x_unit
  real      ( kind = 8 ) fx(2)
  integer   ( kind = 4 ) i
  integer   ( kind = 4 ) j
  integer   ( kind = 4 ), parameter :: n = 48
  integer   ( kind = 4 ), parameter :: nvar = 3
  integer   ( kind = 4 ), parameter :: option = 1
  real      ( kind = 8 ) tan(3)
  real      ( kind = 8 ) tan_norm
  real      ( kind = 8 ) x(3)
  real      ( kind = 8 ) x1
  real      ( kind = 8 ) x2
  real      ( kind = 8 ) x2_start
  real      ( kind = 8 ) x2_stop
  real      ( kind = 8 ) x3

  call timestamp ( )
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) 'P01_OPT01'
  write ( *, '(a)' ) '  Generate points for TEST_CON problem 1, option 1'
  write ( *, '(a)' ) ' '
  write ( *, '(a)' ) '     #       X1              X2             X3' // &
    '              TAN1            TAN2            TAN3' // &
    '             FX1             FX2'
  write ( *, '(a)' ) ' '

  x2_start = -2.0D+00
  x2_stop  =  4.0D+00

  file_x_name = 'p01_opt01_x000.txt'
  file_t_name = 'p01_opt01_t000.txt'

  do j = 0, n

    x2 = ( real ( n - j, kind = 8 ) * x2_start   &
         + real (     j, kind = 8 ) * x2_stop  ) &
         / real ( n,     kind = 8 )

    x1 = ( ( ( -11.0D+00 * x2 + 4.0D+00 ) * x2 + 114.0D+00 ) * x2 + 214.0D+00 ) / 6.0D+00

    x3 = ( ( (             x2 - 2.0D+00 ) * x2 -   6.0D+00 ) * x2 +   4.0D+00 ) / 12.0D+00

    x = (/ x1, x2, x3 /)

    call p01_fun ( option, nvar, x, fx )

    tan(1) = ( ( -33.0D+00 * x2 + 8.0D+00 ) * x2 + 114.0D+00 ) / 6.0D+00
    tan(2) = 1.0D+00
    tan(3) = ( (   3.0D+00 * x2 - 4.0D+00 ) * x2 -   6.0D+00 ) / 12.0D+00

    tan_norm = sqrt ( sum ( tan(1:3)**2 ) )
    tan(1:3) = tan(1:3) / tan_norm

    write ( *, '(2x,i4,2x,g14.6,2x,g14.6,2x,g14.6,2x,g14.6,2x,g14.6,2x,g14.6,2x,g14.6,2x,g14.6)' ) &
      j, x1, x2, x3, tan(1), tan(2), tan(3), fx(1), fx(2)

    call get_unit ( file_x_unit )
    open ( unit = file_x_unit, file = file_x_name, status = 'replace' )
    do i = 1, nvar
      write ( file_x_unit, '(g16.8)' ) x(i)
    end do
    close ( unit = file_x_unit )
    call file_name_inc ( file_x_name )

    call get_unit ( file_t_unit )
    open ( unit = file_t_unit, file = file_t_name, status = 'replace' )
    do i = 1, nvar
      write ( file_t_unit, '(g16.8)' ) tan(i)
    end do
    close ( unit = file_t_unit )
    call file_name_inc ( file_t_name )

  end do

  write ( *,'(a)') ' '
  write ( *, '(a)' ) 'P01_OPT01'
  write ( *, '(a)' ) '  Normal end of execution.'

  write ( *, '(a)' ) ' '
  call timestamp ( )

  stop
end
subroutine file_name_inc ( file_name )

!*****************************************************************************80
!
!! FILE_NAME_INC increments a partially numeric filename.
!
!  Discussion:
!
!    It is assumed that the digits in the name, whether scattered or
!    connected, represent a number that is to be increased by 1 on
!    each call.  If this number is all 9's on input, the output number
!    is all 0's.  Non-numeric letters of the name are unaffected.
!
!    If the name is empty, then the routine stops.
!
!    If the name contains no digits, the empty string is returned.
!
!  Example:
!
!      Input            Output
!      -----            ------
!      'a7to11.txt'     'a7to12.txt'
!      'a7to99.txt'     'a8to00.txt'
!      'a9to99.txt'     'a0to00.txt'
!      'cat.txt'        ' '
!      ' '              STOP!
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    14 September 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input/output, character ( len = * ) FILE_NAME.
!    On input, a character string to be incremented.
!    On output, the incremented string.
!
  implicit none

  character c
  integer change
  integer digit
  character ( len = * ) file_name
  integer i
  integer lens

  lens = len_trim ( file_name )

  if ( lens <= 0 ) then
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'FILE_NAME_INC - Fatal error!'
    write ( *, '(a)' ) '  The input string is empty.'
    stop
  end if

  change = 0

  do i = lens, 1, -1

    c = file_name(i:i)

    if ( lge ( c, '0' ) .and. lle ( c, '9' ) ) then

      change = change + 1

      digit = ichar ( c ) - 48
      digit = digit + 1

      if ( digit == 10 ) then
        digit = 0
      end if

      c = char ( digit + 48 )

      file_name(i:i) = c

      if ( c /= '0' ) then
        return
      end if

    end if

  end do

  if ( change == 0 ) then
    file_name = ' '
    return
  end if

  return
end
subroutine get_unit ( iunit )

!*****************************************************************************80
!
!! GET_UNIT returns a free FORTRAN unit number.
!
!  Discussion:
!
!    A "free" FORTRAN unit number is an integer between 1 and 99 which
!    is not currently associated with an I/O device.  A free FORTRAN unit
!    number is needed in order to open a file with the OPEN command.
!
!    If IUNIT = 0, then no free FORTRAN unit could be found, although
!    all 99 units were checked (except for units 5, 6 and 9, which
!    are commonly reserved for console I/O).
!
!    Otherwise, IUNIT is an integer between 1 and 99, representing a
!    free FORTRAN unit.  Note that GET_UNIT assumes that units 5 and 6
!    are special, and will never return those values.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license. 
!
!  Modified:
!
!    18 September 2005
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Output, integer IUNIT, the free unit number.
!
  implicit none

  integer i
  integer ios
  integer iunit
  logical lopen
 
  iunit = 0
 
  do i = 1, 99
 
    if ( i /= 5 .and. i /= 6 .and. i /= 9 ) then

      inquire ( unit = i, opened = lopen, iostat = ios )
 
      if ( ios == 0 ) then
        if ( .not. lopen ) then
          iunit = i
          return
        end if
      end if

    end if
 
  end do

  return
end
subroutine p01_fun ( option, nvar, x, fx )

!*****************************************************************************80
!
!! P01_FUN evaluates the function for problem 1.
!
!  Title:
!
!    The Freudenstein-Roth function
!
!  Description:
!
!    One way to use a continuation code as a nonlinear root finder
!    is to start with a set of nonlinear equations G(X), and an
!    approximate root A, and create a "homotopy" function F(X,Y)
!    with the properties that F(A,0.0) = 0 and F(X,1.0) = G(X).
!    Thus, the homotopy function F has a known exact solution
!    from which we can start with no difficulty.  If the continuation
!    code can take us from Y = 0 to Y = 1, then we have found
!    an X so that F(X,1.0) = 0, so we have found a solution to G(X)=0.
!
!    The Freudenstein-Roth function F(X) is derived in this way
!    from a homotopy of G(X):
!
!      F ( X(1), X(2), X(3) ) =
!        G ( X(1), X(2) ) - ( 1 - X(3) ) * G ( Y1, Y2 )
!
!    where Y1 and Y2 are some fixed values, and
!
!      G(1) = X(1) - X(2)*X(2)*X(2) + 5*X(2)*X(2) -  2*X(2) - 13
!      G(2) = X(1) + X(2)*X(2)*X(2) +   X(2)*X(2) - 14*X(2) - 29
!
!  Options 1, 2, 3:
!
!    The starting point is X0 = ( 15, -2, 0 ).
!
!    A great deal of information is available about the homotopy curve
!    generated by this starting point:
!
!    The function F(X) has the form
!
!      F(1) = X(1) - X(2)**3 + 5*X(2)**2 -  2*X(2) - 13 + 34*(X(3)-1)
!      F(2) = X(1) + X(2)**3 +   X(2)**2 - 14*X(2) - 29 + 10*(X(3)-1)
!
!    There is a closed form representation of the curve in terms of the
!    second parameter:
!
!      X(1) = (-11*X(2)**3 + 4*X(2)**2 + 114*X(2) + 214) /  6
!      X(2) = X(2)
!      X(3) = (    X(2)**3 - 2*X(2)**2 -   6*X(2) +   4) / 12
!
!    The first option simply requests the production of solution points
!    along the curve until a point is reached whose third component is
!    exactly 1.
!
!    Options 2 and 3 use the same starting point, and also stop when the
!    third component is 1.  However, these options in addition search
!    for limit points in the first and third components of the solution,
!    respectively.
!
!    The target solution has X(3) = 1, and is ( 5, 4, 1 ).
!
!    Limit points for X1:
!
!      ( 14.28309, -1.741377,  0.2585779 )
!      ( 61.66936,  1.983801, -0.6638797 )
!
!    Limit points for X3:
!
!     (20.48586, -0.8968053, 0.5875873)
!     (61.02031,  2.230139, -0.6863528)
!
!    The curve has several dramatic bends.
!
!
!  Options 4, 5, and 6:
!
!    The starting point is (4, 3, 0).
!
!    The function F(X) has the form
!
!      F(1) = X(1) - X(2)**3 + 5*X(2)**2 -  2*X(2) - 13 +  3*(X(3)-1)
!      F(2) = X(1) + X(2)**3 +   X(2)**2 - 14*X(2) - 29 - 31*(X(3)-1)
!
!    There is a closed form representation of the curve in terms of the
!    second parameter:
!
!      X(1) = (14*X(2)**3 -79*X(2)**2 +52*X(2) + 245) / 17
!      X(2) = X(2)
!      X(3) = (   X(2)**3 - 2*X(2)**2 - 6*X(2) +   9) / 17
!
!    The correct value of the solution at X(3)=1 is:
!
!      (5, 4, 1)
!
!    In option 5, limit points in the first component are sought,
!    and in option 6, limit points in the third component are
!    sought.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    12 March 1999
!
!  Author:
!
!    John Burkardt
!
!  Reference:
!
!    Ferdinand Freudenstein, Bernhard Roth,
!    Numerical Solutions of Nonlinear Equations,
!    Journal of the Association for Computing Machinery,
!    Volume 10, 1963, Pages 550-556.
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) OPTION, the option index.
!
!    Input, integer ( kind = 4 ) NVAR, the number of variables.
!
!    Input, real ( kind = 8 ) X(NVAR), the argument of the function.
!
!    Output, real ( kind = 8 ) FX(NVAR-1), the value of the function at X.
!
  implicit none

  integer ( kind = 4 ) nvar

  real    ( kind = 8 ) fx(nvar-1)
  real    ( kind = 8 ) gx(2)
  real    ( kind = 8 ) gy(2)
  integer ( kind = 4 ) i
  integer ( kind = 4 ) option
  real    ( kind = 8 ) x(nvar)
  real    ( kind = 8 ) y(3)
!
!  Get the starting point, Y.
!
  call p01_start ( option, nvar, y )
!
!  G is the function value at the starting point,
!  F the function value at the current point.
!
  call p01_gx ( y, gy )

  call p01_gx ( x, gx )
!
!  The parameter X3 generates the homotopy curve.
!
  fx(1:nvar-1) = gx(1:nvar-1) + ( x(3) - 1.0D+00 ) * gy(1:nvar-1)

  return
end
subroutine p01_gx ( x, g )

!*****************************************************************************80
!
!! P01_GX evaluates the underlying function for problem 1.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    11 December 1998
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, real ( kind = 8 ) X(2), the point at which the function is to
!    be evaluated.
!
!    Output, real ( kind = 8 ) G(2), the value of the function at X.
!
  implicit none

  real    ( kind = 8 ) g(2)
  real    ( kind = 8 ) x(2)

  g(1) = x(1) - ( ( x(2) - 5.0D+00 ) * x(2) + 2.0D+00 ) * x(2) - 13.0D+00
  g(2) = x(1) + ( ( x(2) + 1.0D+00 ) * x(2) - 14.0D+00 ) * x(2) - 29.0D+00

  return
end
subroutine p01_start ( option, nvar, x )

!*****************************************************************************80
!
!! P01_START returns a starting point for problem 1.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    10 March 1999
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer ( kind = 4 ) OPTION, the option index.
!
!    Input, integer ( kind = 4 ) NVAR, the number of variables.
!
!    Output, real ( kind = 8 ) X(NVAR), the starting point.
!
  implicit none

  integer ( kind = 4 ) nvar

  integer ( kind = 4 ) option
  real    ( kind = 8 ) x(nvar)

  if ( option == 1 .or. option == 2 .or. option == 3 ) then
    x(1:3) = (/ 15.0D+00, -2.0D+00, 0.0D+00 /)
  else if ( option == 4 .or. option == 5 .or. option == 6 ) then
    x(1:3) = (/ 4.0D+00, 3.0D+00, 0.0D+00 /)
  else
    write ( *, '(a)' ) ' '
    write ( *, '(a)' ) 'P01_START - Fatal error!'
    write ( *, '(a,i8)' ) '  Unrecognized option  = ', option
    stop
  end if

  return
end
subroutine timestamp ( )

!*****************************************************************************80
!
!! TIMESTAMP prints the current YMDHMS date as a time stamp.
!
!  Example:
!
!    May 31 2001   9:45:54.872 AM
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    31 May 2001
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    None
!
  implicit none

  character ( len = 8 ) ampm
  integer ( kind = 4 ) d
  character ( len = 8 ) date
  integer ( kind = 4 ) h
  integer ( kind = 4 ) m
  integer ( kind = 4 ) mm
  character ( len = 9 ), parameter, dimension(12) :: month = (/ &
    'January  ', 'February ', 'March    ', 'April    ', &
    'May      ', 'June     ', 'July     ', 'August   ', &
    'September', 'October  ', 'November ', 'December ' /)
  integer ( kind = 4 ) n
  integer ( kind = 4 ) s
  character ( len = 10 )  time
  integer ( kind = 4 ) values(8)
  integer ( kind = 4 ) y
  character ( len = 5 ) zone

  call date_and_time ( date, time, zone, values )

  y = values(1)
  m = values(2)
  d = values(3)
  h = values(5)
  n = values(6)
  s = values(7)
  mm = values(8)

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

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

  return
end