!     Last change:  X    16 Nov 2004    1:34 am
PROGRAM superman2
!-----------------------------------------------------------------------------
!                                                           Date: 11-16-04
! Names: Jordon Herzog, Kelly Harris, John Hall, and Alex Beck
!
! Program to calculate the coordinates of the kryptonite using the roots
! of Lex's function
!
! Variables use are:
! roots           : roots of Lex's function
! x0vals, x1vals  : arrays containing the values of x0 and x1
! yValues         : y-coordinates corresponding to the roots (x-coordinates)
! bisectionResult : result of bisection method run on the roots
! secantResult    : result of secant method run on the roots
! i               : loop counter
! numRoots        : constant, number of roots of Lex's function
!
! Input  : values of a and b and the roots of Lex's function
! Output : coordinate locations of the kryptonite
!-----------------------------------------------------------------------------
IMPLICIT NONE

! Declare variables
INTEGER, PARAMETER :: numRoots = 7
REAL (KIND=8), DIMENSION(1:numRoots) :: roots, yValues, x0vals, x1vals
REAL (KIND=8) :: bisectionResult, secantResult
INTEGER :: i

! Initialize variables and arrays
DO i = 1, numRoots
  x0vals = 0.0D0
  x1vals = 0.0D0
  roots(i) = 0.0D0
  yValues(i) = 0.0D0
END DO

WRITE(*,*) "Coordinate Calculation:"
WRITE(*,*)
WRITE(*,*) "Now we can calculate the final coordinates of the locations of the kryptonite!"
WRITE(*,*) "Using the data we created from the extra program..."
WRITE(*,*)

! Open file for writing
OPEN(3, FILE='c:\temp\coords.dat', STATUS='REPLACE')

DO i = 1, numRoots

  ! Get roots from user
  WRITE(*,'(X,A,I1,A)') "Enter root ", i, " (1 being leftmost, 7, rightmost): "
  READ(*,*) roots(i)

  ! Get X0 and X1 from user
  WRITE(*,*) "Enter the associated set of 'x-a' and 'x+b': "
  READ(*,*) x0vals(i), x1vals(i)


  ! Run bisection and secant methods on roots
  bisectionResult = bisectionMethod(x0vals(i), x1vals(i))
  secantResult = secantMethod(x0vals(i), x1vals(i))

  ! Calculate Y coordinates
  yValues(i) = 100.0D0 * ABS(bisectionResult - secantResult)

  ! Write coordinates to file
  WRITE(*,'(1X, A, I2, A \)') "Writing coordinate set ", i, "..."
    WRITE(3,*) "Location ", i, ":"

 ! Convert coordinate to miles with direction
  IF (roots(i) < 0) THEN
    WRITE(3,*) "Go ", ABS(roots(i)) * 30, " miles West"
  ELSE IF (roots(i) >= 0) THEN
    WRITE(3,*) "Go ", roots(i) * 30, " miles East"
  END IF

  IF (yValues(i) < 0) THEN
    WRITE(3,*) "Go ", ABS(yValues(i)) * 30, " miles South"
  ELSE IF (yValues(i) >= 0) THEN
    WRITE(3,*) "Go ", yValues(i) * 30, " miles North"
  END IF

  WRITE(*,*) "Locked in!"
  WRITE(*,*)
  WRITE(*,*)
END DO

WRITE(*,*) "The coordinates of the kryptonite can be found in 'C:\temp\coords.dat'."
WRITE(*,*) "Depechez-vous! Cherchez des kryptonites! Sauvez Superman! Vites!"

! Close file for writing
CLOSE(3)

STOP
CONTAINS
!-----------------------------------------------------------------------------
! 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 = 1

  ! 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

!-----------------------------------------------------------------------------
! bisectionMethod(X0, X1)
!
! Function that implements Bisection Method
!
! Variables used are:
! X0, X1    : bounds of the initial intervals
! 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 = 3

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

  ! Initialize variables
  counter = 1

  ! 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

END PROGRAM superman2
