!     Last change:  KEH  30 Oct 2004   12:31 pm
PROGRAM superman1
!-----------------------------------------------------------------------------
!                                                           Date:
! Names:
!
! Program to....
!
! Variables use are:
!
! Input  :
! Output :
!-----------------------------------------------------------------------------
IMPLICIT NONE

! Declare variables

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, INTENT(IN) :: X0
  REAL :: root, fX, fPrimeX
  REAL, DIMENSION(0:2) :: Xn
  INTEGER :: counter
  INTEGER, PARAMETER :: maxIts = 3

  ! These values needs changing
  REAL, PARAMETER :: tolerance1 = .0001
  REAL, PARAMETER :: tolerance2 = .0001

  ! Initialize variables
  counter = 0

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

  ! 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)
!
! Function that implements Bisection Method
!
! Variables used are:
! X0        : initial root guess
! 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) RESULT(root)

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

  ! This value needs changing
  REAL, PARAMETER :: tolerance = .0001

  ! Initialize variables
  counter = 0

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

  ! 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, INTENT(IN) :: X
  REAL :: 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
! a                : part of function
! aPrime           : part of function
! b1, b2, b3       : parts of function
! bPrime1, bPrime2 : 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, INTENT(IN) :: X
  REAL :: a, b1, b2, b3, aPrime, bPrime1, bPrime2, FpOfX

  ! Implement derivative: FpOfX = a * (bPrime1 + bPrime2) + (aPrime) * (b1 + b2 + b3)

  ! .0025e^(-1.004x^2) - .0005e^(-.004x^2)
  a = (.0025 * EXP(-1.004 * X**2)) - (.0005 * EXP(-.004 * X**2))

  ! .0025e^(-1.004x^2) - .0005e^(-.004x^2)
  aPrime = (.00502 * X * EXP(-1.004 * X**2)) - (.0000004 * X * EXP(-.004 * X**2))

  ! .0005e^(-.004x^2) + 5e^(-x^2) - 1
  b1 = (.0005 * (EXP(-.004 * X**2))) + ((5 * EXP(-X**2)) -1)
  ! x^6 - 54.42192921x^5 - 13.9104024x^4 + 914.7130216x^3
  b2 = X**6 - (54.42192921 * X**5) - (13.9104024 * X**4) + (914.7130216 * X**3)
  ! 502.4509845x^2 - 3846.144357x - 4190.026545
  b3 = (502.4509845 * X**2) - (3846.144357 * X) - 4190.026545

  ! 6x^5 - 272.1096461x^4 - 55.6416096x^3
  bPrime1 = (6 * X**5) - (272.1096461 * X**4) - (55.6416096 * X**3)
  ! 2744.139065x^2 + 1004.901969x - 3846.144357
  bPrime2 = (2744.139065 * X**2) + (1004.901969 * X) - 3846.144357

  ! Complete function
  FpOfX = a * (bPrime1 + bPrime2) + (aPrime) * (b1 + b2 + b3)

  RETURN
END FUNCTION lexFprime

END PROGRAM superman1
