C $Id: ic_grnwch_sidereal.f,v 1.2 1998/07/24 21:54:32 asc Exp $
CCCC
C
C  IC_GRNWCH_SIDEREAL - return the greenwich true sidereal time in radians
C 
C  PURPOSE:  Calculate the true of date greenwich sidereal time in radians.
C
C  UNIT TYPE:  SUBROUTINE
C
C  INVOCATION METHOD:  CALL IC_GRNWCH_MEAN_SIDEREAL (orb_pos_time, 
C                                                    gst)
C
C  ARGUMENT LIST:
C
C  NAME	                  TYPE   USE  DESCRIPTION
C  ----                   ----   ---  -----------
C  ORB_POS_TIME(2)        I*4    I    TIME OF ORB. VECTOR, YEAR-DAY-MILLI OF DAY
C  GST		          R*8    O    GREENWICH MEAN SIDEREAL TIME
C
C  FILE/RECORD REFERENCES:  NONE
C
C  EXTERNAL VARIABLES:  NONE
C
C  EXTERNAL REFERENCES:
C	IC_GET_NUT_ANGLES - Returns angles necessary to adjust the Greenwich
C                           hour angle to true of date
C
C  ABNORMAL TERMINATION CONDITIONS, ERROR MESSAGES:  NONE
C
C  ASSUMPTIONS, CONSTRAINTS, RESTRICTIONS:  NONE
C
C  DEVELOPMENT HISTORY
C
C  AUTHOR	CHANGE ID	RELEASE	  DATE	    DESCRIPTION OF CHANGE
C  ------	---------	-------   ----	    ---------------------
C  J. LUBELCZYK                 B1R1      11/28/90  INITIAL PDL
C  J. LUBELCZYK                 B1R1      12/10/90  CODING
C  C. RAYMOND (CSC)                       07/23/91  CHANGED TO J2000 SYSTEM
C                                                    AND SIMPLIFIED THE 
C                                                    COMPUTATION OF TIME
C
C  J. Lubelczyk				  09/12/91  Added code to convert GST
C						    to true of date
C  J. LUBELCZYK ICCR #83, CCR #'S 130, 137 11/91    B3 update
C
C  NOTES:
C  1)  THE ORIGINAL ALGORITHM USED WAS COPIED FROM A SHORT PROGRAM BY 
C      G. D. MEAD, INCLUDED IN 'GEOPHYSICAL COORDINATE 
C                       TRANSFORMATIONS' BY CHRISTOPHER T. RUSSELL
C  2)  THIS VERSION INCORPORATES SEVERAL CHANGES TO CALCULATE THE GREENWICH
C      MEAN SIDEREAL TIME CORRECTLY ON THE J2000 COORDINATE SYSTEM.  THE
C      PREVIOUS VERSION WAS ONLY CORRECT IN THE B1950 COORDINATE SYSTEM.
C  3)  NOW RETURNS THE TRUE OF DATE GREENWICH SIDEREAL TIME ON THE J2000 SYS
C
CCCC
C
C  PDL:
C
C  Convert ORB_POS_TIME to year, day, and seconds of day
C
C  Calculate GST [greenwich sidereal time] in radians
C
C  Call IC_GET_NUT_ANGLES to get angles to adjust the GMST to true of date
C
C  Convert GMST to GST true of date
C
C  RETURN
C
CCCC
	subroutine IC_GRNWCH_SIDEREAL (ORB_POS_TIME, GST)
C
	implicit none
C
C*  Calling parameters
C
	INTEGER*4   orb_pos_time(2) !ISTP time format (YYYYDDD, Milli of Day)
	REAL*8      gst             !Greenwich Sidreal Time
C
C*  Other Variables
C
 	INTEGER*4   year      !4 digit year
	INTEGER*4   day	      !day of year
	REAL*8      secs      !Seconds of day
	REAL*8	    dj        !Julian Day
	REAL*8	    fday      !fraction of day
	REAL*8	    t         !time in centuries
	REAL*8	    half  /0.50D0/
        real*8      C0    /1.7533685592332653D0/     !Polynomial Coef.
        real*8      C1    /628.33197068884084D0/     !Polynomial Coef.
        real*8      C2    /0.67707139449033354D-05/  !Polynomial Coef.
        real*8      C3    /6.3003880989848915D0/     !Polynomial Coef.
        real*8      C4    /-0.45087672343186841D-09/ !Polynomial Coef.
        real*8      TWOPI /6.283185307179586D0/      !Two PI
	real*8	    DELEPS	!Delta Epsilon
	real*8	    DELPSI	!Delta Psi
	real*8	    EPS		!Epsilon
C
C*  Start Executable
C
        YEAR = ORB_POS_TIME(1)/1000
        DAY = MOD(ORB_POS_TIME(1),1000)
C
C*  Convert the given millisecond of day [ORB_POS_TIME(2)] to second of day.
C
	SECS = (real(ORB_POS_TIME(2)))/1000.0d0
C
C**  Begin calculating the greenwich mean sidereal time **
C
	FDAY = SECS/86400.00
	DJ = 365*(YEAR-1900)+(YEAR-1901)/4+DAY-HALF
C
C	      THE NEXT STATEMENT CAUSES THE REFERENCE EPOCH TO BE SHIFTED	
C	   TO THE J2000 REFERENCE EPOCH.
C
	T = (DJ-36525.0D0)/36525.0D0 
        GST   = C0 + T*(C1 + T*(C2 + C4*T)) + C3*FDAY           
        GST  = DMOD( GST,TWOPI)                                          
        IF (GST .LT. 0.0D0) GST = GST + TWOPI
C
C*  Convert GST to true of date by adjusting for nutation
C
	call IC_GET_NUT_ANGLES (T, DELEPS, DELPSI, EPS)
	GST = GST + DELPSI*DCOS(EPS+DELEPS)
        RETURN                                                            
        END                                                               
