!     Last change:  A    26 Jul 2005    1:24 am
PROGRAM superman1
!-----------------------------------------------------------------------------
!                                                           Date: 11-16-04
! Names: Jordon Herzog, Kelly Harris, John Hall, and Alex Beck
!
! Program to find the roots of Lex's function and write them to a file
!
! Variables use are: 
! lowBound     : lower bound of the range containing the root used in 
!                Bisection Method and Secant Method
! highBound    : upper bound of the range containing the root used in 
!                Bisection Method and Secant Method
! initialGuess : initial root guess used in Newton's Method, Householder's
!                Method, Halley's Method, Schroder's Method, and the 
!                polynomial root approximation methods
! bResult      : result of Bisection Method
! nResult      : result of Newton's Method
! hResult      : result of Halley's method
! sResult      : result of Schroder's method
! hhResult     : result of Householder's method
! scResult     : result of Secant method
! poly1Result  : result from the First Polynomial Root Approximation Method 
! poly4Result  : result from the Fourth Polynomial Root Approximation Method
! poly2Result  : result from the Two Polynomial Root Approximation Method
! roots        : an array storing the roots of the Lex's function
! i            : counter variable
! runAgain     : logical value that determines whether to run through 
!                loops again
!
! Input  : initial bounds and guess for a root
! Output : values of the roots in a file
!-----------------------------------------------------------------------------
IMPLICIT NONE

! Declare variables
REAL (KIND=8) :: lowBound, highBound, initialGuess, bResult, nResult
REAL (KIND=8) :: hResult, sResult, hhResult, scResult
REAL (KIND=8) :: poly1Result, poly4Result, poly2Result
REAL (KIND=8), DIMENSION(1:7) :: roots
INTEGER :: i
LOGICAL :: runAgain

runAgain = .TRUE.

WRITE(*,*) "This program is Superman's LAST HOPE! All the kryptonite MUST be found "
WRITE(*,*) "or we will never see Superman again (really he only exists in the movies anyway but...)"
WRITE(*,*)

! Find root -2.99999998395320
WRITE(*,*) "Root 1"
DO
  WRITE(*,*) "TIP: Try looking between -3.5 and -2.5."
  WRITE(*,*)

  ! Get interval for Bisection Method from user
  WRITE(*,*) "What range do you think the root is in? "
  READ(*,*) lowBound, highBound

  ! Get initial guess for Newton's Method from user
  WRITE(*,*) "Within that range, what do you think the value of the root is? "
  READ(*,*) initialGuess

  ! Calculate root approximations
  nResult = newtonsMethod(initialGuess)
  bResult = bisectionMethod(lowBound, highBound)
  hhResult = Householders(initialGuess)

  WRITE(*,*) "Newton's Method result:      ", nResult
  WRITE(*,*) "Bisection Method result:     ", bResult
  WRITE(*,*) "Householder's Method result: ", hhResult
  WRITE(*,*)

  IF (ABS(nResult - hhResult) < 1.0D-6 .AND. ABS(bResult - hhResult) < 1.0D-6) THEN
    WRITE(*,*) "Root 1 was found."
    roots(1) = bResult
    runAgain = .FALSE.
  END IF

  IF (.NOT. runAgain) EXIT

END DO

WRITE(*,*)
WRITE(*,*) "--------------------------------------------------------------------"
WRITE(*,*)

! Find root -1.61364659473
WRITE(*,*) "Root 2"
runAgain = .TRUE.
DO
  WRITE(*,*) "TIP: Try looking between  -2.0 and -1.6135."
  WRITE(*,*)

  ! Get interval for Bisection Method from user
  WRITE(*,*) "What range do you think the root is in? "
  READ(*,*) lowBound, highBound

  ! Get initial guess for Newton's Method from user
  WRITE(*,*) "Within that range, what do you think the value of the root is? "
  READ(*,*) initialGuess

  ! Calculate root approximations
  nResult = newtonsMethod(initialGuess)
  bResult = bisectionMethod(lowBound, highBound)
  poly1Result = poly1(initialGuess)

  WRITE(*,*) "Newton's Method result:      ", nResult
  WRITE(*,*) "Bisection Method result:     ", bResult
  WRITE(*,*) "First Polynomial result:     ", poly1Result
  WRITE(*,*)

 IF (ABS(nResult - poly1Result) < 1.0D-2 .AND. ABS(bResult - poly1Result) < 1.0D-2) THEN
    WRITE(*,*) "Root 2 was found."
    roots(2) = bResult
    runAgain = .FALSE.
  END IF

  IF (.NOT. runAgain) EXIT

END DO

WRITE(*,*)
WRITE(*,*) "--------------------------------------------------------------------"
WRITE(*,*)

! Find root -1.613153416584
WRITE(*,*) "Root 3"
runAgain = .TRUE.
DO
  WRITE(*,*) "TIP: Try looking between -1.6135 and -1.3."
  WRITE(*,*)

  ! Get interval for Bisection Method from user
  WRITE(*,*) "What range do you think the root is in? "
  READ(*,*) lowBound, highBound

  ! Get initial guess for Newton's Method from user
  WRITE(*,*) "Within that range, what do you think the value of the root is? "
  READ(*,*) initialGuess

  ! Calculate root approximations
  nResult = newtonsMethod(initialGuess)
  bResult = bisectionMethod(lowBound, highBound)
  sResult = Schroders(initialGuess)

  WRITE(*,*) "Newton's Method result:      ", nResult
  WRITE(*,*) "Bisection Method result:     ", bResult
  WRITE(*,*) "Schroder's Method result:    ", sResult
  WRITE(*,*)

  IF (ABS(nResult - sResult) < 1.0D-6 .AND. ABS(bResult - sResult) < 1.0D-6) THEN
    WRITE(*,*) "Root 3 was found."
    roots(3) = bResult
    runAgain = .FALSE.
  END IF

  IF (.NOT. runAgain) EXIT

END DO

WRITE(*,*)
WRITE(*,*) "--------------------------------------------------------------------"
WRITE(*,*)

! Find root -1.26863624117952
WRITE(*,*) "Root 4"
runAgain = .TRUE.
DO
  WRITE(*,*) "TIP: Try looking between -1.5 and -1.0."
  WRITE(*,*)

  ! Get interval for Bisection Method from user
  WRITE(*,*) "What range do you think the root is in? "
  READ(*,*) lowBound, highBound

  ! Get initial guess for Newton's Method from user
  WRITE(*,*) "Within that range, what do you think the value of the root is? "
  READ(*,*) initialGuess

  ! Calculate root approximations
  nResult = newtonsMethod(initialGuess)
  bResult = bisectionMethod(lowBound, highBound)
  hResult = halleysMethod(initialGuess)

  WRITE(*,*) "Newton's Method result:      ", nResult
  WRITE(*,*) "Bisection Method result:     ", bResult
  WRITE(*,*) "Halley's Method result:      ", hResult
  WRITE(*,*)

  IF (ABS(nResult - hResult) < 1.0D-4 .AND. ABS(bResult - hResult) < 1.0D-4) THEN
    WRITE(*,*) "Root 4 was found."
    roots(4) = bResult
    runAgain = .FALSE.
  END IF

  IF (.NOT. runAgain) EXIT

END DO

WRITE(*,*)
WRITE(*,*) "--------------------------------------------------------------------"
WRITE(*,*)

! Find root 1.26863624117952
WRITE(*,*) "Root 5"
runAgain = .TRUE.
DO
  WRITE(*,*) "TIP: Try looking between 0.0 and 2.0."
  WRITE(*,*)

  ! Get interval for Bisection Method from user
  WRITE(*,*) "What range do you think the root is in? "
  READ(*,*) lowBound, highBound

  ! Get initial guess for Newton's Method from user
  WRITE(*,*) "Within that range, what do you think the value of the root is? "
  READ(*,*) initialGuess

  ! Calculate root approximations
  nResult = newtonsMethod(initialGuess)
  bResult = bisectionMethod(lowBound, highBound)
  scResult = secantMethod(lowBound, highBound)

  WRITE(*,*) "Newton's Method result:      ", nResult
  WRITE(*,*) "Bisection Method result:     ", bResult
  WRITE(*,*) "Secant Method result:        ", scResult
  WRITE(*,*)

  IF (ABS(nResult - scResult) < 1.0D-6 .AND. ABS(bResult - scResult) < 1.0D-6) THEN
    WRITE(*,*) "Root 5 was found."
    roots(5) = bResult
    runAgain = .FALSE.
  END IF

  IF (.NOT. runAgain) EXIT

END DO

WRITE(*,*)
WRITE(*,*) "--------------------------------------------------------------------"
WRITE(*,*)

! Find root 3.141300595407
WRITE(*,*) "Root 6"
runAgain = .TRUE.
DO
  WRITE(*,*) "TIP: Try looking between 2.5 and 3.1415."
  WRITE(*,*)

  ! Get interval for Bisection Method from user
  WRITE(*,*) "What range do you think the root is in? "
  READ(*,*) lowBound, highBound

  ! Get initial guess for Newton's Method from user
  WRITE(*,*) "Within that range, what do you think the value of the root is? "
  READ(*,*) initialGuess

  ! Calculate root approximations
  nResult = newtonsMethod(initialGuess)
  bResult = bisectionMethod(lowBound, highBound)
  poly2Result = poly2(initialGuess)

  WRITE(*,*) "Newton's Method result:      ", nResult
  WRITE(*,*) "Bisection Method result:     ", bResult
  WRITE(*,*) "Second Polynomial result:    ", poly2Result
  WRITE(*,*)

  IF (ABS(nResult - poly2Result) < 1.0D-2 .AND. ABS(bResult - poly2Result) < 1.0D-2) THEN
    WRITE(*,*) "Root 6 was found."
    roots(6) = bResult
    runAgain = .FALSE.
  END IF

  IF (.NOT. runAgain) EXIT

END DO

WRITE(*,*)
WRITE(*,*) "--------------------------------------------------------------------"
WRITE(*,*)

! Find root 3.141792059897
WRITE(*,*) "Root 7"
runAgain = .TRUE.
DO
  WRITE(*,*) "TIP: Try looking between 3.1415 and 4.0."
  WRITE(*,*)

  ! Get interval for Bisection Method from user
  WRITE(*,*) "What range do you think the root is in? "
  READ(*,*) lowBound, highBound

  ! Get initial guess for Newton's Method from user
  WRITE(*,*) "Within that range, what do you think the value of the root is? "
  READ(*,*) initialGuess

  ! Calculate root approximations
  nResult = newtonsMethod(initialGuess)
  bResult = bisectionMethod(lowBound, highBound)
  poly4Result = poly4(initialGuess)

  WRITE(*,*) "Newton's Method result:      ", nResult
  WRITE(*,*) "Bisection Method result:     ", bResult
  WRITE(*,*) "Fourth Polynomial result:    ", poly4Result
  WRITE(*,*)

  IF (ABS(nResult - poly4Result) < 1.0D-2 .AND. ABS(bResult - poly4Result) < 1.0D-2) THEN
    WRITE(*,*) "Root 7 was found."
    roots(7) = bResult
    runAgain = .FALSE.
  END IF

  IF (.NOT. runAgain) EXIT

END DO

WRITE (*,*) "Now go to the Range Calculation program as the first step towards finding the coordinates."

! Write roots to c:\temp\roots.dat
OPEN(3, FILE="c:\temp\roots.dat", STATUS="REPLACE")

DO i = 1,7
  WRITE(3,'(1X, A, I1, A \ )') "Root ", i, ": "
  WRITE(3,*) roots(i)
END DO

CLOSE(3)

STOP

CONTAINS
!-----------------------------------------------------------------------------
! newtonsMethod(X0)
!
! Function that implements Newton's Method
!
! Variables used are:
! X0         : initial root guess
! Xn         : array with the necessary X values for the iteration
! fX         : value of Lex's function evaluated at X
! fPrimeX    : value of the derivative of Lex's function evaluated at X
! root       : root approximation after iterations are complete
! counter    : loop counter
! maxIts     : constant, maximum number of loop iterations
! tolerance1 : constant, error tolerance value for |X0 - X1|
! tolerance2 : constant, error tolerance value for F(X0)
!
! Paremeter  : initial root guess
! Return     : final root approximation
!-----------------------------------------------------------------------------
FUNCTION newtonsMethod(X0) RESULT(root)

  ! Declare variables
  REAL (KIND=8), INTENT(IN) :: X0
  REAL (KIND=8) :: root, fX, fPrimeX
  REAL (KIND=8), DIMENSION(0:2) :: Xn
  INTEGER :: counter
  INTEGER, PARAMETER :: maxIts = 1000

  ! These values needs changing
  REAL (KIND=8), PARAMETER :: tolerance1 = .1D-15
  REAL (KIND=8), PARAMETER :: tolerance2 = .1D-15

  ! Initialize variables
  counter = 0

  ! Put X0 and X1 into Xn array
  Xn(0) = X0


  fX = lexF(Xn(0))
  fPrimeX = lexFprime(Xn(0))

  ! Begin Newton's Method
  IF ((fX /= 0) .AND. (fPrimeX /= 0)) THEN
    ! Newton'st Method iteration
    DO
      ! Calculate F(x) and F'(X)
      fX = lexF(Xn(0))
      fPrimeX = lexFprime(Xn(0))

      ! Compute next root approximation
      Xn(1) = Xn(0)
      Xn(0) = Xn(0) - fX/fPrimeX

      ! Exit condition
      IF ((counter > maxIts) .OR. (ABS(Xn(0) - Xn(1)) < tolerance1) .OR. (ABS(lexF(Xn(0))) < tolerance2)) EXIT

      counter = counter + 1
    END DO

  END IF
  root = Xn(0)

    RETURN
END FUNCTION newtonsMethod

!-----------------------------------------------------------------------------
! bisectionMethod(X0, X1)
!
! Function that implements Bisection Method
!
! Variables used are:
! X0, X1    : initial interval
! Xn        : array with the necessary X values for the iteration
! root      : root approximation after iterations are complete
! counter   : loop counter
! maxIts    : constant, maximum number of loop iterations
! tolerance : constant, error tolerance value for root approximation
!
! Paremeter  : initial root guess
! Return     : final root approximation
!-----------------------------------------------------------------------------
FUNCTION bisectionMethod(X0, X1) RESULT(root)

  ! Declare variables
  REAL (KIND=8), INTENT(IN) :: X0, X1
  REAL (KIND=8) :: root
  REAL (KIND=8), DIMENSION(0:2) :: Xn
  INTEGER :: counter
  INTEGER, PARAMETER :: maxIts = 1000

  ! This value needs changing
  REAL, PARAMETER :: tolerance = .1D-15

  ! Initialize variables
  counter = 0

  ! Put X0 into Xn array
  Xn(0) = X0
  Xn(1) = X1

  ! Bisection Method iteration
  DO
    Xn(2) = (Xn(0) + Xn(1)) / 2

    IF ((lexF(Xn(2)) * lexF(Xn(0))) < 0) THEN
      Xn(1) = Xn(2)
    ELSE
      Xn(0) = Xn(2)
    END IF

    ! Exit condition
    IF ((counter > maxIts) .OR. (ABS(Xn(0) - Xn(1)) < 2 * tolerance) .OR. (lexF(Xn(2)) == 0)) EXIT

    counter = counter + 1
  END DO

  root = Xn(2)

    RETURN
END FUNCTION bisectionMethod

!-----------------------------------------------------------------------------
! lexF(X)
!
! Function that implements Lex's function
!
! Variables use are:
! X      : value at which to evaluate the function
! FofX   : the result of the function evaluated at X
! FofX_1 : part of the function
! FofX_2 : part of the function
! FofX_3 : part of the function
! FofX_4 : part of the function
!
! Paremeter  : value at which to evaluate the function
! Return     : the result of the function evaluated at X
!-----------------------------------------------------------------------------
FUNCTION lexF(X) RESULT(FofX)

  ! Declare variables
  REAL (KIND=8), INTENT(IN) :: X
  REAL (KIND=8) :: FofX, FofX_1, FofX_2, FofX_3, FofX_4

  ! Implement function:  FofX = FofX_1 * FofX_2 * (FofX_3 + FofX_4)

  ! .0005e^(-.004x^2)
  FofX_1 = .0005 * (EXP(-.004 * (X**2)))
  ! 5e^(-x^2) - 1
  FofX_2 = (5 * EXP(-X**2)) -1
  ! x^6 - 54.42192921x^5 - 13.9104024x^4 + 914.7130216x^3
  FofX_3 = (X**6) - 54.42192921 * (X**5) -13.9104024 * (X**4) + 914.7130216 * (X**3)
  ! 502.4509845x^2 - 3846.144357x - 4190.026545
  FofX_4 = 502.4509845 * (X**2) - 3846.144357 * X - 4190.026545

  ! Complete function
  FofX = FofX_1 * FofX_2 * (FofX_3 + FofX_4)

  RETURN
END FUNCTION lexF

!-----------------------------------------------------------------------------
! lexFprime(X)
!
! Function that implements the derivative of Lex's function
!
! Variables use are:
! X                : value at which to evaluate the function
! FpOfX            : the result of the function evaluated at X
! fp0, fp1, fp2    : parts of function
! fp3, fp4, fp5    : parts of function
! fp6, fp7         : parts of function
!
! Paremeter  : value at which to evaluate the function
! Return     : the result of the function evaluated at X
!-----------------------------------------------------------------------------
FUNCTION lexFprime(X) RESULT(FpOfX)

  ! Declare variables
  REAL (KIND=8), INTENT(IN) :: X
  REAL (KIND=8) :: fp0, fp1, fp2, fp3, fp4, fp5, fp6, fp7, FpOfX

  ! Implement derivative: FpOfX = fp0 + fp1 + fp2 + fp3 +fp4 +fp5 + fp6 + fp7

  fp0 = -9.615360893D0 * EXP(-1.004D0 * X**2) + 1.923072179D0 * EXP(-.004D0 * X**2)
  fp1 = (23.54618818D0 * EXP(-1.004D0 * X**2) * X) - (0.5192110907D0 * EXP(-.004D0 * X**2) * X)
  fp2 = (26.16799233D0 * EXP(-1.004D0 * X**2) * X**2) - (1.38745411D0 * EXP(-.004D0 * X**2) * X**2)
  fp3 = -(2.661407966D0 * EXP(-1.004D0 * X**2) * X**3) + (0.02983060874D0 * EXP(-.004D0 * X**2) * X**3)
  fp4 = -(5.272133484D0 * EXP(-1.004D0 * X**2) * X**4) + (0.1397136751D0 * EXP(-.004D0 * X**2) * X**4)
  fp5 = (0.08483022005D0 * EXP(-1.004D0 * X**2) * X**5) - (0.003055641610D0 * EXP(-.004D0 * X**2) * X**5)
  fp6 = (0.2731980846D0 * EXP(-1.004D0 * X**2) * X**6) - (0.0002176877168D0 * EXP(-.004D0 * X**2) * X**6)
  fp7 = -(0.00502D0 * EXP(-1.004D0 * X**2) * X**7) + (0.000004D0 * EXP(-.004D0 * X**2) * X**7)

  ! Complete function
  FpOfX = fp0 + fp1 + fp2 + fp3 +fp4 +fp5 + fp6 + fp7

  RETURN



END FUNCTION lexFprime

!-----------------------------------------------------------------------------
! lexFDoublePrime(X)
!
! Function that implements the second derivative of Lex's function
!
! Variables use are:
! X                : value at which to evaluate the function
! FpOfX            : the result of the function evaluated at X
! fdp0, fdp1, fdp2 : parts of function
! fdp3, fdp4, fdp5 : parts of function
! fdp6, fdp7, fdp8 : parts of function
!
! Paremeter  : value at which to evaluate the function
! Return     : the result of the function evaluated at X
!-----------------------------------------------------------------------------
FUNCTION lexFDoublePrime(X) RESULT(FDpOfX)

  ! Declare variables
  REAL (KIND=8), INTENT(IN) :: X
  REAL (KIND=8) :: FDpOfX
  REAL (KIND=8) :: fdp0, fdp1, fdp2, fdp3, fdp4, fdp5, fdp6, fdp7, fdp8

  ! Implement second derivative:
  ! FDpOfX = fdp0 + fdp1 + fdp2 + fdp3 +fdp4 +fdp5 + fdp6 + fdp7 + fdp8

  fdp0 = -(.5192110906800 * EXP(-0.004 * X**2)) + (23.5461881784000 * EXP(-1.004 * X**2))
  fdp1 = -(2.7902927970840 * EXP(-0.004 * X**2) * X) + (71.6436293404200 * EXP(-1.004*X**2) * X)
  fdp2 = (0.0936455149394400 * EXP(-0.004*X**2) * X**2) - (55.2649697607972000 * EXP(-1.004*X**2) )* X**2
  fdp3 = (.5699543333242240 * EXP(-0.004*X**2) * X**3) - (73.6338625411811200 * EXP(-1.004*X**2) * X**3)
  fdp4 = -(0.01551685291790400 * EXP(-0.004*X**2) * X**4) + (5.76825829634952000 * EXP(-1.004*X**2) * X**4)
  fdp5 = -(0.00242383570193120 * EXP(-0.004*X**2) * X**5) + (12.22563254278765600 * EXP(-1.004*X**2) * X**5)
  fdp6 = (0.00005244513287680 * EXP(-0.004*X**2) * X**6) - (.20547908185638400 * EXP(-1.004*X**2) * X**6)
  fdp7 = -(.548581753945473600 * EXP(-1.004*X**2) * X**7) + (0.000001741501734720 * EXP(-0.004*X**2) * X**7)
  fdp8 = -(0.0000000320 * EXP(-0.004*X**2) * X**8) + (0.0100801600 * EXP(-1.004*X**2) * X**8)

  ! Complete function
  FDpOfX = fdp0 + fdp1 + fdp2 + fdp3 +fdp4 +fdp5 + fdp6 + fdp7 + fdp8

  RETURN

END FUNCTION lexFDoublePrime

!-----------------------------------------------------------------------------
! halleysMethod(X)
!
! Function that implements Halley's Method for the function 5e^(-x^2) - 1
!
! Variables use are:
! X                : value at which to evaluate the function
! Fp2OfX            : the result of the function evaluated at X
! fp0, fp1, fp2    : parts of function
! fp3, fp4, fp5    : parts of function
! fp6, fp7         : parts of function
!
! Paremeter  : value at which to evaluate the function
! Return     : the result of the function evaluated at X
!-----------------------------------------------------------------------------
FUNCTION halleysMethod(X) RESULT(root)

  ! Declare variables
  REAL (KIND=8), INTENT(IN) :: X
  REAL (KIND=8) :: root, f, fp1, fp2
  REAL (KIND=8), DIMENSION(0:1) :: Xn
  INTEGER :: counter
  INTEGER, PARAMETER :: maxIts = 100000
  REAL, PARAMETER :: tolerance = 1.D-15

  ! Initialize counter variable
  counter = 0

  ! Put X in Xn array
  Xn(0) = X

  DO

    ! Define function and its derivatives
    f = (5 * EXP(-Xn(0)**2)) - 1
    fp1 = -10 * Xn(0) * EXP(-Xn(0)**2)
    fp2 = (20 * Xn(0)**2 - 10) * EXP(-Xn(0)**2)

    ! Halley's Method iteration
    Xn(1) = Xn(0) - (2 * f * fp1)/(2 * fp1**2 - f * fp2)

    ! Exit condition
    IF ((counter > maxIts) .OR. (ABS(Xn(0) - Xn(1)) < tolerance)) EXIT

    counter = counter + 1

  END DO

  root = Xn(1)

  RETURN
END FUNCTION halleysMethod

!-----------------------------------------------------------------------------
! poly6(X)
!
! Derelict function that we used to do preliminary testing for convergence
! to roots
!
! Variables use are:
! x       : value at which to evaluate the function
! x0      : the iteration dummy variable
! a, b, c : coefficients of x^5, x^4, and x^3, respectively
! d, e    : coefficients of x^2, x, respectively
! g       : constant term
! f       : result of the function evaluated at x0
!
! Paremeter  : value at which to evaluate the function
! Return     : the result of the function evaluated at X
!-----------------------------------------------------------------------------
FUNCTION poly6(x) RESULT(f)

   IMPLICIT NONE
   REAL(KIND=8), INTENT(IN) :: x
   REAL(KIND=8) :: f, a, b, c, d, e, g, x0

   a = 54.42192921D0
   b = 13.91040240D0
   c = 914.7130216D0
   d = 502.4509845D0
   e = 3846.144357D0
   g = 4190.026545D0

   x0 = x

   DO
      f = ((a * (x0 ** 5)) + (b * (x0 ** 4)) - (c * (x0 ** 3)) - (d * (x0 ** 2)) + (e * x0) + g) ** (1.0D0 / 6.0D0)
      IF (ABS(x0 - f) <= 0.1D-15) EXIT
      x0 = f
   END DO

RETURN
END FUNCTION poly6

!-----------------------------------------------------------------------------
! poly5(X)
!
! Derelict function that we used to do preliminary testing for convergence
! to roots
!
! Variables use are:
! x       : value at which to evaluate the function
! x0      : the iteration dummy variable
! a, b, c : coefficients of x^5, x^4, and x^3, respectively
! d, e    : coefficients of x^2, x, respectively
! g       : constant term
! f       : result of the function evaluated at x0
!
! Paremeter  : value at which to evaluate the function
! Return     : the result of the function evaluated at X
!-----------------------------------------------------------------------------
FUNCTION poly5(x) RESULT(f)

   IMPLICIT NONE
   REAL(KIND=8), INTENT(IN) :: x
   REAL(KIND=8) :: f, a, b, c, d, e, g, x0

   a = 54.42192921D0
   b = 13.91040240D0
   c = 914.7130216D0
   d = 502.4509845D0
   e = 3846.144357D0
   g = 4190.026545D0

   x0 = x

   DO
      f= ((((x0 ** 6)) - (b * (x0 ** 4)) + (c * (x0 ** 3)) + (d * (x0 ** 2)) - (e * x0) - g) / a) ** (1.0D0 / 5.0D0)
      IF (ABS(x0 - f) <= 0.1D-7) EXIT
      x0 = f
   END DO

RETURN
END FUNCTION poly5

!-----------------------------------------------------------------------------
! poly4(X)
!
! Root finding function derived from the fourth degree iteration map
! of Lex's function
!
! Variables use are:
! x       : value at which to evaluate the function
! x0      : the iteration dummy variable
! a, b, c : coefficients of x^5, x^4, and x^3, respectively
! d, e    : coefficients of x^2, x, respectively
! g       : constant term
! f       : result of the function evaluated at x0
!
! Paremeter  : value at which to evaluate the function
! Return     : the result of the function evaluated at X
!-----------------------------------------------------------------------------
FUNCTION poly4(x) RESULT(f)

   IMPLICIT NONE
   REAL(KIND=8), INTENT(IN) :: x
   REAL(KIND=8) :: f, a, b, c, d, e, g, x0

   a = 54.42192921D0
   b = 13.91040240D0
   c = 914.7130216D0
   d = 502.4509845D0
   e = 3846.144357D0
   g = 4190.026545D0

   x0 = x

   DO
      f= ((((x0 ** 6)) - (a * (x0 ** 5)) + (c * (x0 ** 3)) + (d * (x0 ** 2)) - (e * x0) - g) / b) ** (1.0D0 / 4.0D0)
      IF (ABS(x0 - f) <= 0.1D-7) EXIT
      x0 = f
   END DO

RETURN
END FUNCTION poly4

!-----------------------------------------------------------------------------
! poly3(X)
!
! Derelict function that we used to do preliminary testing for convergence
! to roots
!
! Variables use are:
! x       : value at which to evaluate the function
! x0      : the iteration dummy variable
! a, b, c : coefficients of x^5, x^4, and x^3, respectively
! d, e    : coefficients of x^2, x, respectively
! g       : constant term
! f       : result of the function evaluated at x0
!
! Paremeter  : value at which to evaluate the function
! Return     : the result of the function evaluated at X
!-----------------------------------------------------------------------------
FUNCTION poly3(x) RESULT(f)

   IMPLICIT NONE
   REAL(KIND=8), INTENT(IN) :: x
   REAL(KIND=8) :: f, a, b, c, d, e, g, x0

   a = 54.42192921D0
   b = 13.91040240D0
   c = 914.7130216D0
   d = 502.4509845D0
   e = 3846.144357D0
   g = 4190.026545D0

   x0 = x

   DO
      f= ((((a * (x0 ** 5)) - (x0 ** 6)) + (b * (x0 ** 4)) - (d * (x0 ** 2)) + (e * x0) + g) / c) ** (1.0D0 / 3.0D0)
      IF (ABS(x0 - f) <= 0.1D-15) EXIT
      x0 = f
   END DO

RETURN
END FUNCTION poly3

!-----------------------------------------------------------------------------
! poly2(X)
!
! Root finding function derived from the second degree iteration map
! of Lex's function
!
! Variables use are:
! x       : value at which to evaluate the function
! x0      : the iteration dummy variable
! a, b, c : coefficients of x^5, x^4, and x^3, respectively
! d, e    : coefficients of x^2, x, respectively
! g       : constant term
! f       : result of the function evaluated at x0
!
! Paremeter  : value at which to evaluate the function
! Return     : the result of the function evaluated at X
!-----------------------------------------------------------------------------
FUNCTION poly2(x) RESULT(f)

   IMPLICIT NONE
   REAL(KIND=8), INTENT(IN) :: x
   REAL(KIND=8) :: f, a, b, c, d, e, g, x0

   a = 54.42192921D0
   b = 13.91040240D0
   c = 914.7130216D0
   d = 502.4509845D0
   e = 3846.144357D0
   g = 4190.026545D0

   x0 = x

   DO
      f= ((((a * (x0 ** 5)) - (x0 ** 6)) + (b * (x0 ** 4)) - (c * (x0 ** 3)) + (e * x0) + g) / d) ** (1.0D0 / 2.0D0)
      IF ((ABS(x0 - f) <= 0.1D-15)) EXIT
      x0 = f
   END DO

RETURN
END FUNCTION poly2

!-----------------------------------------------------------------------------
! poly1(X)
!
! Root finding function derived from the first degree iteration map
! of Lex's function
!
! Variables use are:
! x       : value at which to evaluate the function
! x0      : the iteration dummy variable
! a, b, c : coefficients of x^5, x^4, and x^3, respectively
! d, e    : coefficients of x^2, x, respectively
! g       : constant term
! f       : result of the function evaluated at x0
!
! Paremeter  : value at which to evaluate the function
! Return     : the result of the function evaluated at X
!-----------------------------------------------------------------------------
FUNCTION poly1(x) RESULT(f)

   IMPLICIT NONE
   REAL(KIND=8), INTENT(IN) :: x
   REAL(KIND=8) :: f, a, b, c, d, e, g, x0

   a = 54.42192921D0
   b = 13.91040240D0
   c = 914.7130216D0
   d = 502.4509845D0
   e = 3846.144357D0
   g = 4190.026545D0

   x0 = x

   DO
      f = ((x ** 6) - (a * (x ** 5)) - (b * (x ** 4)) + (c * (x ** 3)) + (d * (x ** 2)) - g) / e
      IF (ABS(x0 - f) <= 0.1D-7) EXIT
      x0 = f
   END DO

RETURN
END FUNCTION poly1

!------------------------------------------------------------------------------
! Householder's Method
!
! Function that implements Householder's Method of root finding.
!
! Variables used are:
! X0         : initial root guess
! Xn         : array with the necessary X values for the iteration
! root       : root approximation after iterations are complete
! fX         : value of the function evaluated at Xn(0)
! fPrimeX    : value of the derivative of the function evaluated at Xn(0)
! fDPrimeX   : value of the second derivative of the function evaluated at Xn(0)
! counter    : loop counter
! maxIts     : constant, maximum number of loop iterations
! tolerance1 : constant, error tolerance value for difference
!              between the current and previous approximations
! tolerance2 : constant, error tolerance value for the function
!              evaluted at the current approximation
!
! Paremeter  : initial root guesses
! Return     : final root approximation
!-----------------------------------------------------------------------------

FUNCTION Householders(X0) RESULT(root)

 ! Declare variables
  REAL(KIND=8), INTENT(IN) :: X0
  REAL(KIND=8) :: root, fX, fPrimeX, fDPrimeX
  REAL(KIND=8), DIMENSION(0:2) :: Xn
  INTEGER :: counter
  INTEGER, PARAMETER :: maxIts = 1000

  ! Tolerance values
  REAL(KIND=8), PARAMETER :: tolerance1 = .1D-15
  REAL(KIND=8), PARAMETER :: tolerance2 = .1D-15

  ! Initialize variables
  counter = 0

  ! Put X0 into Xn array
  Xn(0) = X0

  fX = lexF(Xn(0))
  fPrimeX = lexFprime(Xn(0))
  fDPrimeX = lexFDoublePrime(Xn(0))

  ! Begin Schroder's Method
  IF ((fX /= 0) .AND. (fPrimeX /= 0) .AND. (fDPrimeX /= 0)) THEN
    ! Schroder's Method iteration
    DO
      ! Calculate F(x), F'(X), and F''(X)
      fX = lexF(Xn(0))
      fPrimeX = lexFprime(Xn(0))
      fDPrimeX = lexFDoublePrime(Xn(0))

      ! Compute next root approximation
      Xn(1) = Xn(0)
      Xn(0) = Xn(0) - (fX/fPrimeX) * (1 + (fX *  fDPrimeX) / 2 * (fPrimeX**2))

      ! Exit condition
      IF ((counter > maxIts) .OR. (ABS(Xn(0) - Xn(1)) < tolerance1) .OR. (ABS(lexF(Xn(0))) < tolerance2)) EXIT

      counter = counter + 1
    END DO

  END IF
  root = Xn(0)

RETURN
END FUNCTION Householders

!------------------------------------------------------------------------------
! Schroder's Method
!
! Function that implements Schroder's Method of root finding.
!
! Variables used are:
! X0         : initial root guess
! Xn         : array with the necessary X values for the iteration
! root       : root approximation after iterations are complete
! fX         : value of the function evaluated at Xn(0)
! fPrimeX    : value of the derivative of the function evaluated at Xn(0)
! fDPrimeX   : value of the second derivative of the function evaluated at Xn(0)
! counter    : loop counter
! maxIts     : constant, maximum number of loop iterations
! tolerance1 : constant, error tolerance value for difference
!              between the current and previous approximations
! tolerance2 : constant, error tolerance value for the function
!              evaluted at the current approximation
!
! Paremeter  : initial root guesses
! Return     : final root approximation
!-----------------------------------------------------------------------------
FUNCTION Schroders(X0) RESULT(root)

 ! Declare variables
  REAL(KIND=8), INTENT(IN) :: X0
  REAL(KIND=8) :: root, fX, fPrimeX, fDPrimeX
  REAL(KIND=8), DIMENSION(0:2) :: Xn
  INTEGER :: counter
  INTEGER, PARAMETER :: maxIts = 1000

  ! Tolerance values
  REAL(KIND=8), PARAMETER :: tolerance1 = .1D-15
  REAL(KIND=8), PARAMETER :: tolerance2 = .1D-15

  ! Initialize variables
  counter = 0

  ! Put X0 and X1 into Xn array
  Xn(0) = X0


  fX = lexF(Xn(0))
  fPrimeX = lexFprime(Xn(0))
  fDPrimeX = lexFDoublePrime(Xn(0))

  ! Begin Schroder's Method
  IF ((fX /= 0) .AND. (fPrimeX /= 0) .AND. (fDPrimeX /= 0)) THEN
    ! Schroder's Method iteration
    DO
      ! Calculate F(x), F'(X), and F''(X)
      fX = lexF(Xn(0))
      fPrimeX = lexFprime(Xn(0))
      fDPrimeX = lexFDoublePrime(Xn(0))

      ! Compute next root approximation
      Xn(1) = Xn(0)
      Xn(0) = Xn(0) - ((fX)*(fPrimeX))/(((fPrimeX)**2)-((fX)*(fDPrimeX)))

      ! Exit condition
      IF ((counter > maxIts) .OR. (ABS(Xn(0) - Xn(1)) < tolerance1) .OR. (ABS(lexF(Xn(0))) < tolerance2)) EXIT

      counter = counter + 1
    END DO

  END IF
  root = Xn(0)

RETURN
END FUNCTION Schroders

!-----------------------------------------------------------------------------
! secantMethod(X0, X1)
!
! Function that implements Secant Method
!
! Variables used are:
! X0, X1    : initial root guesses
! Xn        : array with the necessary X values for the iteration
! tempX     : temporary X value
! root      : root approximation after iterations are complete
! counter   : loop counter
! maxIts    : constant, maximum number of loop iterations
! tolerance : constant, error tolerance value for root approximation
!
! Paremeter  : initial root guesses
! Return     : final root approximation
!-----------------------------------------------------------------------------
FUNCTION secantMethod(X0, X1) RESULT(root)

  ! Declare variables
  REAL (KIND=8), INTENT(IN) :: X0, X1
  REAL (KIND=8) :: tempX, root
  REAL (KIND=8), DIMENSION(0:2) :: Xn
  INTEGER :: counter
  INTEGER, PARAMETER :: maxIts = 3

  ! This value needs changing
  REAL (KIND=8), PARAMETER :: tolerance = .0001

  ! Initialize variables
  counter = 0

  ! Put X0 and X1 into Xn array
  Xn(0) = X0
  Xn(1) = X1

  IF (ABS(lexF(Xn(0))) < ABS(lexF(Xn(1)))) THEN
    tempX = Xn(0)
    Xn(0) = Xn(1)
    Xn(1) = tempX
  END IF

  ! Secant Method iteration
  DO
    Xn(2) = Xn(1) - lexF(Xn(1)) * (Xn(0) - Xn(1))/(lexF(Xn(0)) - lexF(Xn(1)))
    Xn(0) = Xn(1)
    Xn(1) = Xn(2)

    ! Exit condition
    IF ((counter > maxIts) .OR. (ABS(lexF(Xn(2))) < tolerance)) EXIT

    counter = counter + 1
  END DO

  root = Xn(2)

    RETURN
END FUNCTION secantMethod

END PROGRAM superman1

