!***********************************************************************************************************************************
!
!                                                          J D 2 G M S T
!
!  Program:      JD2GMST
!
!  Programmer:   David G. Simpson
!                NASA Goddard Space Flight Center
!                Greenbelt, Maryland  20771
!
!  Date:         February 18, 2007
!
!  Language:     Fortran-90
!
!  Version:      1.00a
!
!  Description:  Convert Julian Day to Greenwich Mean Sidereal Time (GMST).
!
!  Files:        Source files:
!
!                   jd2gmst.f90                   Main program
!
!  Notes:        Algorithm from "Astronomical Algorithms" by Jean Meeus.
!
!***********************************************************************************************************************************

      PROGRAM JD2GMST

      IMPLICIT NONE

      DOUBLE PRECISION :: JD, T, GMST, GMST_H
      CHARACTER(LEN=15) :: GMST_STR, TEMP_STR

      DOUBLE PRECISION :: FRAC, REDUCE_D

!-----------------------------------------------------------------------------------------------------------------------------------

!
!     Get user input (Julian Day).
!

      WRITE (UNIT=*, FMT='(A)', ADVANCE='NO') ' Enter Julian day:  '                ! get user input
      READ (UNIT=*, FMT=*) JD

!
!     Compute Greenwich Mean Sidereal Time (GMST).
!

      T = (JD-2451545.0D0)/36525.0D0                                                ! convert JD to Julian centuries

      GMST = 280.46061837D0 + 360.98564736629D0*(JD-2451545.0D0) +  &               ! GMST (deg)
         3.87933D-4*T**2 - T**3/3.871D7

      GMST = REDUCE_D(GMST, 0.0D0)                                                  ! reduce GMST to range 0-360 deg

      GMST_H = GMST / 15.0D0                                                        ! convert GMST to hours of arc

!
!     Convert GMST hours to hh:mm:ss and print to string GMST_STR.
!

      WRITE (UNIT=TEMP_STR, FMT='(I2.2)') INT(GMST_H)                               ! hours
      GMST_STR(1:2) = TEMP_STR(1:2)

      GMST_STR(3:3) = ':'

      GMST_H = FRAC(GMST_H) * 60.0D0                                                ! minutes
      WRITE (UNIT=TEMP_STR, FMT='(I2.2)') INT(GMST_H)
      GMST_STR(4:5) = TEMP_STR(1:2)

      GMST_STR(6:6) = ':'

      GMST_H = FRAC(GMST_H) * 60.0D0                                                ! seconds
      WRITE (UNIT=TEMP_STR, FMT='(I2.2)') INT(GMST_H)
      GMST_STR(7:8) = TEMP_STR(1:2)

      GMST_STR(9:9) = '.'

      GMST_H = FRAC(GMST_H) * 1.0D6                                                 ! microseconds
      WRITE (UNIT=TEMP_STR, FMT='(I6.6)') INT(GMST_H)
      GMST_STR(10:15) = TEMP_STR(1:6)

!
!     Print results.
!

      WRITE (UNIT=*, FMT='(/A,F13.8,A)') ' GMST = ', GMST, ' deg'                   ! degrees
      WRITE (UNIT=*, FMT='(A)') '      = '//GMST_STR                                ! hh:mm:ss

      STOP

      END PROGRAM JD2GMST





!***********************************************************************************************************************************
!  FRAC
!***********************************************************************************************************************************

      FUNCTION FRAC (X) RESULT (Y)

      IMPLICIT NONE

      DOUBLE PRECISION, INTENT(IN) :: X
      DOUBLE PRECISION :: Y, Z

      Z = ABS(X)
      Y = Z - INT(Z)
      Y = SIGN(Y,X)

      RETURN

      END FUNCTION FRAC





!***********************************************************************************************************************************
!     REDUCE_D
!
!     Reduce a given angle THETA (in radians) to the range [ANGLE_MIN, ANGLE_MIN+360).
!***********************************************************************************************************************************

      FUNCTION REDUCE_D (THETA, ANGMIN) RESULT (Y)

      IMPLICIT NONE

      DOUBLE PRECISION, INTENT(IN) :: THETA, ANGMIN
      DOUBLE PRECISION :: Y

      DOUBLE PRECISION  ANGMAX, REVS


      ANGMAX = ANGMIN + 360.0D0

      IF (THETA .LT. ANGMIN) THEN
         REVS = AINT((ANGMIN-THETA)/360.0D0) + 1
         Y = THETA + REVS*360.0D0
      ELSE IF (THETA .GE. ANGMAX) THEN
         REVS = AINT((THETA-ANGMIN)/360.0D0)
         Y = THETA - REVS*360.0D0
      ELSE
         Y = THETA
      END IF

      RETURN

      END FUNCTION REDUCE_D