ic_gci_transf.pro (Documentation for /home/davin/idl/socware/)

;
; IC_GCI_TRANSF
;
; Routines to transform from GCI coordinates to GEO coordinates.
; The following ISTP ICSS routines are appended in this file.
; 
; IC_CONV_MATRIX.FOR
; IC_GCI_TO_GEO.FOR
; IC_GET_NUT_ANGLES.FOR
; IC_GRNWCH_SIDEREAL.FOR
;
; In addition, the NAG routine F01CKF is emulated by a subroutine
; included in this file.  Consequently, this code should not be
; linked with the NAG libraries if available.  The IC_ routines
; have not been modified in any fashion.
;
; A good reference for the routines below is 'An Explanatory
; Supplement to the Astronomical Almanac,' P. Kenneth Seidelmann, ed.
; University Science Books, 1992.  ISBN 0-935702-68-7
;
; G. GERMANY   8/9/95
;
;------------------------------------------------------------------


;
; IC_GCI_TO_GEO - return a transformation matrix
;
; PURPOSE:  Calculate the transformation matrix from GCI 
;           coordinates to GEO coordinates at a given date and time.
;
; UNIT TYPE:  SUBROUTINE
;
; INVOCATION METHOD:  CALL IC_GCI_TO_GEO (orb_pos_time, 
;                                         transform_matrix)
;
; argUMENT LIST:
;
; NAME	                  TYPE   USE  DESCRIPTION
; ----                   ----   ---  -----------
; orb_pos_time(2)        I*4    I    time OF ORB. VECTOR, year-day-MILLI OF day
; transform_matrix(3,3)  R*8    O    TRANSFORMATION MATRIX
;
; FILE/RECORD REFERENCES:  NONE
;
; EXTERNAL VARIABLES:  NONE
;
; EXTERNAL REFERENCES:
;	IC_GRNWCH_SIDEREAL - Returns the Greenwich sidereal time in radians
;	F01CKF - Multiplies two matrices
;	IC_CONV_MATRIX - Returns the conversion matrix to rotate
;                       from mean of date to true of date
;
; ABNORMAL TERMINATION CONDITIONS, ERROR MESSAGES:  NONE
;
; ASSUMPTIONS, CONstrAINTS, REstrICTIONS:  NONE
;
; DEVELOPMENT HISTORY
;
; AUTHOR	CHANGE ID	RELEASE	  DATE	    DESCRIPTION OF CHANGE
; ------	---------	-------   ----	    ---------------------
; J. LUBELCZYK                 B1R1      11/21/90  INITIAL PDL
; J. Lubelczyk			B1R1      12/10/90  CODING
; J. Lubelczyk                           09/18/91  Updated to return true
;                                                  of date trans matrix
; J. LUBELCZYK ICCR #83, CCR #'S 130, 137 11/91    B3 update
;
; NOTES:
;

PRO ic_gci_to_geo, orb_pos_time, transform_matrix


        mean_matrix = DBLARR(3,3)
	transform_matrix = DBLARR(3,3) ;transformation matrix
	cmatrix = DBLARR(3,3)	      ;the conversion matrix to rotate from 
                                      ;mean of date to true of date

	ic_grnwch_sidereal, orb_pos_time, grnwch_sidereal_time

;
;   calculate the sin and cos of the greenwich mean sidereal time
;
	sin_gst = sin(grnwch_sidereal_time)
	cos_gst = cos(grnwch_sidereal_time)

;
;   Fill the mean of date transformation matrix using the sin and cos
;    of the greenwich mean sidereal time
;
	mean_matrix[0,0] = cos_gst
	mean_matrix[0,1] = -sin_gst
	mean_matrix[0,2] = 0
	mean_matrix[1,0] = sin_gst
	mean_matrix[1,1] = cos_gst
	mean_matrix[1,2] = 0
	mean_matrix[2,0] = 0
	mean_matrix[2,1] = 0
	mean_matrix[2,2] = 1

	ic_conv_matrix,orb_pos_time, cmatrix

        transform_matrix = TRANSPOSE( $
                           TRANSPOSE(mean_matrix) # TRANSPOSE(cmatrix) )

END


;      STATEMENT FUNCTION DEFINITION FOR dxjul -- JULIAN EPHEMERIS
;      DATE AT BEGINNING OF year.
;

FUNCTION dxjul,i

	RETURN,DOUBLE((-32075+1461*(i+4800-13/12)/4  $
                     +367*(-1+13/12*12)/12-3         $
                     *((i+4900-13/12)/100)/4)-0.5	)

END


; IC_CONV_MATRIX - Returns the conversion matrix that is necessary to rotate
;		    from mean of date to true of date
;
; PURPOSE:  THIS SUBROUTINE CALCULATES, THROUGH APPROPRIATE ANALYTIC
;           EXPRESSIONS, VALUES FOR THE PRECESSION AND NUTATION
;           ANGLES AND THE MATRIX REQUIRED TO ROTATE FROM MEAN OF
;           JULIAN 2000 TO TRUE OF DATE.
; NAME	               TYPE   USE  DESCRIPTION
; ----                ----   ---  -----------
; orb_pos_time(2)     I*4    I    time OF ORB. VECTOR, year-day-MILLI OF day
; CMATRIX(3,3)        R*8    O    Matrix to rotate from J2000 to true of date
;
; EXTERNAL REFERENCES:
;	F01CKF - NAG routine that multiplies two matrices
;      IC_GET_NUT_ANGLES          Routine to compute the nutation angles 
;                            

PRO ic_conv_matrix,orb_pos_time, cmatrix

	cmatrix = DBLARR(3,3)    ;conversion matrix
	nutmat = DBLARR(3,3)	    ;Nutation matrix
	premat = DBLARR(3,3)	    ;precession matrix

	Jdj2000 = DOUBLE (2451545.0)
        R = DOUBLE (1296000.0)


;   Convert the given millisecond of day [orb_pos_time[1]] to second of day.
;   Convert the packed form into year and day-of-year

      secs = (DOUBLE(orb_pos_time[1]))/DOUBLE(1000.0)
      year = orb_pos_time[0]/1000
      day  = orb_pos_time[0] MOD 1000

;
;   Calculate the julian date and the time in Julian centuries from J2000
;
      fday = secs/DOUBLE(86400.00)
      jul_day = dxjul(year) + DOUBLE(day)+fday
      time = (jul_day - Jdj2000)/DOUBLE(36525.0)

      T2 = time*time
      T3 = time*T2

;    CALCULATE CONVERSION FACTORS: DEGREES TO RADANS (dtr), SECONDS TO 
;    RADIANS (str)                                                     
;                                                                      
      PI= DOUBLE(4.0 * ATAN(1.0))
      dtr=PI/DOUBLE(180.0)
      str=dtr/DOUBLE(3600.0)


;    CALCULATE PRECESSION ANGLES

      zeta   = DOUBLE(                     	$
               0.11180860865024398D-01*time	$
             + 0.14635555405334670D-05*T2	$
             + 0.87256766326094274D-07*T3 )
      theta  = DOUBLE(                     	$
               0.97171734551696701D-02*time	$
             - 0.20684575704538352D-05*T2	$
             - 0.20281210721855218D-06*T3 )
      zee    = DOUBLE(                     	$
               0.11180860865024398D-01*time	$
             + 0.53071584043698687D-05*T2	$
             + 0.88250634372368822D-07*T3 )

      sinzet = sin(zeta)                                               
      coszet = cos(zeta)                                               
      sinzee = sin(zee)   
      coszee = cos(zee)   
      sinthe = sin(theta) 
      costhe = cos(theta) 
;                                                                      
;    COMPUTE THE TRANSFORMATION MATRIX BETWEEN MEAN EQUATOR AND        
;    EQUINOX OF 1950 AND MEAN EQUATOR AND EQUINOX OF DATE. THIS        
;    MATRIX IS CALLED premat.                                          
;                                                                      
      premat[0,0] = -sinzet*sinzee  + coszet*coszee*costhe              
      premat[0,1] =  coszee*sinzet  + sinzee*costhe*coszet              
      premat[0,2] =  sinthe*coszet                                      
      premat[1,0] = -sinzee*coszet  - coszee*costhe*sinzet              
      premat[1,1] =  coszee*coszet  - sinzee*costhe*sinzet              
      premat[1,2] = -sinthe*sinzet                                      
      premat[2,0] = -coszee*sinthe                                      
      premat[2,1] = -sinzee*sinthe                                      
      premat[2,2] =  costhe                                             

;    CALCULATE MEAN OBLIQUITY OF DATE (epso). WHERE TIME IS MEASURED IN
;    JULIAN CENTURIES FROM 2000.0.

      epso=DOUBLE( (1.813E-3*T3-5.9E-4*T2	$
                    -4.6815E+1*time+8.4381448E+4)*str )


;    CALL IC_GET_NUT_ANGLES TO COMPUTE NUTATION IN OBLIQUITY AND LONGITUDE

      ic_get_nut_angles,time,deleps,delpsi,eps

      cosep=cos(eps)                               
      cosepO=cos(epso)                             
      cospsi=cos(delpsi)                           
      sinep=sin(eps)                               
      sinepO=sin(epso)                             
      sinpsi=sin(delpsi)                           

      nutmat[0,0]=cospsi                            
      nutmat[1,0]=-sinpsi*cosepO                    
      nutmat[2,0]=-sinpsi*sinepO                    
      nutmat[0,1]=sinpsi*cosep                      
      nutmat[1,1]=cospsi*cosep*cosepO+sinep*sinepO  
      nutmat[2,1]=cospsi*cosep*sinepO-sinep*cosepO  
      nutmat[0,2]=sinpsi*sinep                      
      nutmat[1,2]=cospsi*sinep*cosepO-cosep*sinepO  
      nutmat[2,2]=cospsi*sinep*sinepO+cosep*cosepO  

;    CALCULATE ELEMENTS OF nutmat * premat.  THIS MATRIX IS THE 
;    ANALYTICALLY CALCULATED TRANSFORMATION MATRIX, WHICH WILL  
;    TRANSFORM THE MEAN EARTH EQUATOR AND EQUINOX OF J2000 INTO
;    THE TRUE EARTH EQUATOR AND EXQUINOX OF DATE.

;     cmatrix = nutmat # premat
      cmatrix = TRANSPOSE(TRANSPOSE(premat) # TRANSPOSE(nutmat))

END




; IC_GET_NUT_ANGLES - Returns angles that are necessary to adjust the
;		       Greenwich Hour angle to true of date
;
; PURPOSE : THIS SUBROUTINE CALCULATES, THROUGH APPROPRIATE ANALYTIC
;           EXPRESSIONS, VALUES FOR THE NUTATION
;           ANGLES TO ROTATE FROM MEAN OF
;           JULIAN 2000 TO TRUE OF DATE.
;
; NAME	      TYPE    USE	DESCRIPTION
; ----       ----    ---       -----------
; time       R*8     I         time IN JULIAN CENTURIES OF 36525.0
;                              MEAN SOLAR dayS FROM J2000. (NOTE: THIS
;                              CAN BE POSITIVE OR NEGATIVE.)
; deleps     R*8     O	        DELTA epsILON, Nutation in obliquity
; delpsi     R*8     O         DELTA PSI, Nutation in longitude
; eps        R*8     O         epsILON



PRO ic_get_nut_angles,time,deleps,delpsi,eps


	TOL = DOUBLE(0.5)/DOUBLE(36525.0)


        isinco = FLTARR(2,106)   ;Array used in nutation calculations
	icosco = FLTARR(2,106)   ;Array used in nutation calculations
	fund = DBLARR(1,5)	    ;The fundamental arguments
	T = DBLARR(1,2)	    ;Time
	ifunar = intarr(5,106)   ;array used in nutation calculations
        oldtim = DOUBLE(0.0)



;                                                                     
;    INITIALIZE VALUES OF IFUNAR, isinco, AND icosco ARRAYS FOR USE IN
;    NUTATION CALCULATIONS.                                           
;                                                                     
;    THE 1980 IAU THEORY OF NUTATION,CONTAINED IN JPL
;    DE200 PLANETARY EPHEMERIS.


      ifunar[0,0:79] = [0,0,2,-2,2, 0,2,-2,2,0, 2,2,2,0,2, 0,0,2,0,0, $
                2,0,2,0,0, -2,-2,0,0,2, 2,0,2,2,0, 2,0,0,0,2, $
                2,2,0,2,2,  2,2,0,0,2,  0,2,2,2,0, 2,0,2,2,0, $
                0,2,0,-2,0, 0,2,2,2,0,  2,2,2,2,0, 0,0,2,0,0]
      ifunar[0,80:105] = [2,2,0,2,2,  2,4,0,2,2,  0,4,2,2,2, 0,-2,2,0,-2, $
                2,0,-2,0,2, 0]

      ifunar[1,0:79] = [1,2,1,0,2, 0,1,1,2,0,  2,2,1,0,0, 0,1,2,1,1, $
                1,1,1,0,0, 1,0,2,1,0,  2,0,1,2,0, 2,0,1,1,2, $
                1,2,0,2,2, 0,1,1,1,1,  0,2,2,2,0, 2,1,1,1,1, $
                0,1,0,0,0, 0,0,2,2,1,  2,2,2,1,1, 2,0,2,2,0]
      ifunar[1,80:105] = [2,2,0,2,1, 2,2,0,1,2,  1,2,2,0,1, 1,1,2,0,0, $
                1,1,0,0,2, 0]

      ifunar[2,0:79] = [0,0,0,0,0, -1,-2,0,0,1,  1,-1,0,0,0, 2,1,2,-1,0, $
               -1,0,1,0,1,  0,1,1,0,1,   0,0,0,0,0,  0,0,0,0,0, $
                0,0,0,0,0,  0,0,0,0,0,	 1,1,-1,0,0, 0,0,0,0,0, $
               -1,0,1,0,0,  1,0,-1,-1,0, 0,-1,1,0,0, 0,0,0,0,0]
      ifunar[2,80:105] = [0,0,0,1,0,  0,0,-1,0,0,  0,0,0,0,1, -1,0,0,1,0,	$
               -1,1,0,0,0,  1]

      ifunar[3,0:79] = [0,0,-2,2,-2, 1,0,2,0,0,    0,0,0,2,0,   0,0,0,0,-2, $
                0,2,0,1,2,   0,0,0,-1,0,   0,1,0,1,1,  -1,0,1,-1,-1, $
                1,0,2,1,2,   0,-1,-1,1,-1, 1,0,0,1,1,   2,0,0,1,0, $
                1,2,0,1,0,   1,1,1,-1,-2,  3,0,1,-1,2,  1,3,0,-1,1]
      ifunar[3,80:105] = [-2,-1,2,1,1, -2,-1,1,2,2,   1,0,3,1,0,  -1,0,0,0,1, $
                0,1,1,2,0,   0]

      ifunar[4,0:79] = [0,0,0,0,0,    -1,-2,0,-2,0, -2,-2,-2,-2,-2, 0,0,-2,0,2, $
               -2,-2,-2,-1,-2, 2,2,0,1,-2,   0,0,0,0,-2,    0,2,0,0,2, $
                0,2,0,-2,0,    0,0,2,-2,2,  -2,0,0,2,2,    -2,2,2,-2,-2, $
                0,0,-2,0,1,    0,0,0,2,0,    0,2,0,-2,0,    0,0,1,0,-4]
      ifunar[4,80:105] = [2,4,-4,-2,2,   4,0,-2,-2,2,  2,-2,-2,-2,0,  2,0,-1,2,-2, $
                0,-2,2,2,4,    1]


      isinco[0,0:79] = [-171996.,2062.,46.,11.,-3.,   -3.,-2.,1.,-13187.,1426.,	$
                -517.,217.,129.,48.,-22.,     17.,-15.,-16.,-12.,-6., $
                -5.,4.,4.,-4.,1.,              1.,-1.,1.,1.,-1., $
                -2274.,712.,-386.,-301.,-158., 123.,63.,63.,-58.,-59., $
                -51.,-38.,29.,29.,-31.,        26.,21.,16.,-13.,-10., $
                -7.,7.,-7.,-8.,6.,             6.,-6.,-7.,6.,-5., $
                 5.,-5.,-4.,4.,-4.,           -3.,3.,-3.,-3.,-2., $
                -3.,-3.,2.,-2.,2.,            -2.,2.,2.,1.,-1.]
       isinco[0,80:105] = [1.,-2.,-1.,1.,-1.,           -1.,1.,1.,1.,-1., $
                -1.,1.,1.,-1.,1.,              1.,-1.,-1.,-1.,-1., $
                -1.,-1.,-1.,1.,-1.,            1.]

       isinco[1,0:38] = [-174.2,.2,0.,0.,0.,            0.,0.,0.,-1.6,-3.4, $
                 1.2,-.5,.1,0.,0.,            -.1,0.,.1,0.,0., $
                 0.,0.,0.,0.,0.,               0.,0.,0.,0.,0., $
               -.2,.1,-.4,0.,0.,               0.,0.,.1,-.1]

      icosco[0,0:79] = [92025.,-895.,-24.,0.,1.,   0.,1.,0.,5736.,54., $
                224.,-95.,-70.,1.,0.,      0.,9.,7.,6.,3., $
                3.,-2.,-2.,0.,0.,          0.,0.,0.,0.,0., $
                977.,-7.,200.,129.,-1.,   -53.,-2.,-33.,32.,26., $
                27.,16.,-1.,-12.,13.,     -1.,-10.,-8.,7.,5., $
                0.,-3.,3.,3.,0.,          -3.,3.,3.,-3.,3., $
                0.,3.,0.,0.,0.,            0.,0.,1.,1.,1., $
                1.,1.,-1.,1.,-1.,          1.,0.,-1.,-1.,0.]
      icosco[0,80:105] = [-1.,1.,0.,-1.,1.,           1.,0.,0.,-1.,0., $
                0.,0.,0.,0.,0.,            0.,0.,0.,0.,0., $
                0.,0.,0.,0.,0.,            0.]

      icosco[1,0:33] = [8.9,.5,0.,0.,0.,           0.,0.,0.,-3.1,-.1, $
                -.6,.3,0.,0.,0.,           0.,0.,0.,0.,0., $
                 0.,0.,0.,0.,0.,           0.,0.,0.,0.,0.,	$
               -.5,0.,0.,-.1]

      R = DOUBLE(1296000.0)

      IF ( ABS ( time - oldtim ) LE TOL ) THEN BEGIN
         deleps = olddep 
         delpsi = olddps 
         eps    = oldeps
         RETURN
      ENDIF

      T2 = time*time
      T3 = time*T2

;    CONVERT IFUNAR, isinco, AND icosco ARRAYS TO REAL*8 ARRAYS FUNarg,
;    sincof, AND coscof, RESPECTIVELY.                                 
;                                                                      
      funarg = DOUBLE(ifunar)
      sincof = DOUBLE(isinco)
      coscof = DOUBLE(icosco)


;    CALCULATE CONVERSION FACTORS: DEGREES TO RADANS (dtr), SECONDS TO 
;    RADIANS (str)                                                     

      PI= DOUBLE(4.0 * ATAN(1.0))
      dtr=PI/DOUBLE(180.0)
      str=dtr/DOUBLE(3600.0)

;    BEGIN COMPUTATION OF NUTATION IN OBLIQUITY AND LONGITUDE          

;    CALCULATE FUNDAMENTAL argUMENTS FOR USE IN NUTATION CALCULATIONS
;    time IS REFERENCED TO J2000.0.
;    fund(1,1)= F
;    fund(2,1)= OMEGA
;    fund(3,1)= L PRIME
;    fund(4,1)= L
;    fund(5,1)= D
	
      fund[0,0]=str*(335778.877E0+(1342.0E0*R+295263.137E0)*time	$
      -13.257E0*T2+1.1E-2*T3)
      fund[0,1]=str*(450160.280E0-(5.E0*R+482890.539E0)*time+	$
      7.455E0*T2+8.0E-3*T3)
      fund[0,2]=str*(1287099.804E0+(99.0E0*R+1292581.224E0)*time-	$
      5.77E-1*T2-1.2E-2*T3)
      fund[0,3]=str*(485866.733E0+(1325.0E0*R+715922.633E0)*time+	$
      31.310E0*T2+6.4E-2*T3)
      fund[0,4]=str*(1072261.307E0+(1236.0E0*R+1105601.328E0)*time-	$
      6.891E0*T2+1.9E-2*T3)


;    CALCULATE MEAN OBLIQUITY OF DATE (epso). WHERE time IS MEASURED IN
;    JULIAN CENTURIES FROM 2000.0.

      epso=(1.813E-3*T3-5.9E-4*T2-4.6815E+1*time+8.4381448E+4)*str


;                                                                      
;    CALCULATE NUTATION IN LONGITUDE (delpsi) AND NUTATION IN OBLIQUITY
;    (deleps).  THIS IS A THREE STEP PROCESS:                          
;    (1) CALCULATE argUMENTS OF sinE (FOR delpsi) AND COsinE (FOR deleps)
;        THESE ARE OF THE FORM                                        
;                                                                     
;        arg = SUMMATION ( A(I) * fund(I,1) ), I = 1,5                 
;                                                                      
;        WHERE THE A(I)'S ARE ELEMENTS OF FUNarg.                      
;                                                          
;      arg = funarg # fund
;      arg = fund # funarg
      arg = TRANSPOSE(TRANSPOSE(funarg) # TRANSPOSE(fund))

;                                                                      
;    (2) CALCULATE COEFFICIENTS OF sinE AND COsinE, WHICH ARE THE PRODUCTS
;        OF sincof * T AND coscof * T.  THESE COEFFICIENTS ARE IN UNITS
;        OF 0.0001 SECONDS OF ARC.                                     
;                                                                      
      T[0,0]=DOUBLE(1.0)
      T[0,1]=time       
                         
;      cofcos = coscof # T
;      cofsin = sincof # T                       
;      cofcos = T # coscof
;      cofsin = T # sincof
      cofcos = TRANSPOSE(TRANSPOSE(coscof) # TRANSPOSE(T))
      cofsin = TRANSPOSE(TRANSPOSE(sincof) # TRANSPOSE(T))

      cofcos=cofcos*DOUBLE(1.E-4)
      cofsin=cofsin*DOUBLE(1.E-4)

;                                                                      
;    (3) CALCULATE THE sinES AND COsinES OF THE argUMENTS AND MULTIPLY 
;        BY THEIR COEFFICIENTS, THEN ADD.  COMPUTE delpsi AND deleps.  
;                                                                      
      sumpsi=DOUBLE(0.0)
      sumeps=DOUBLE(0.0)

      sinp=sin(arg)
      cose=cos(arg)

      FOR E=0,105 DO BEGIN
         prodps=cofsin[0,E]*sinp[0,E]                                    
         prodep=cofcos[0,E]*cose[0,E]                                    
         sumpsi=sumpsi+prodps                                              
         sumeps=sumeps+prodep  
      ENDFOR

      deleps=sumeps*str                                                 
      delpsi=sumpsi*str                  

;                                        
;    CALCULATE TRUE OBLIQUITY OF DATE (eps).       
;                                                  
      eps=epso+deleps
      olddep = deleps
      olddps = delpsi
      oldeps = eps
      oldtim = time

END



; IC_GRNWCH_SIDEREAL - return the greenwich true sidereal time in radians
;
; PURPOSE:  Calculate the true of date greenwich sidereal time in radians.
;
; NAME	                  TYPE   USE  DESCRIPTION
; ----                   ----   ---  -----------
; orb_pos_time(2)        I*4    I    time OF ORB. VECTOR, year-day-MILLI OF day
; gst		          R*8    O    GREENWICH MEAN SIDEREAL time
; EXTERNAL REFERENCES:
;	IC_GET_NUT_ANGLES - Returns angles necessary to adjust the Greenwich
;                          hour angle to true of date
; NOTES:
; 1)  THE ORIGINAL ALGORITHM USED WAS COPIED FROM A SHORT PROGRAM BY 
;     G. D. MEAD, INCLUDED IN 'GEOPHYSICAL COORDINATE 
;                      TRANSFORMATIONS' BY CHRISTOPHER T. RUSSELL
; 2)  THIS VERSION INCORPORATES SEVERAL CHANGES TO CALCULATE THE GREENWICH
;     MEAN SIDEREAL time CORRECTLY ON THE J2000 COORDINATE SYSTEM.  THE
;     PREVIOUS VERSION WAS ONLY CORRECT IN THE B1950 COORDINATE SYSTEM.
; 3)  NOW RETURNS THE TRUE OF DATE GREENWICH SIDEREAL time ON THE J2000 SYS

PRO ic_grnwch_sidereal,orb_pos_time, gst

	half  = DOUBLE(0.50)
        C0    = DOUBLE(1.7533685592332653)     ;Polynomial Coef.
        C1    = DOUBLE(628.33197068884084)     ;Polynomial Coef.
        C2    = DOUBLE(0.67707139449033354E-05)  ;Polynomial Coef.
        C3    = DOUBLE(6.3003880989848915)     ;Polynomial Coef.
        C4    = DOUBLE(-0.45087672343186841E-09) ;Polynomial Coef.
        TWOPI = DOUBLE(6.283185307179586)      ;Two PI


        year = orb_pos_time[0]/1000
        day  = orb_pos_time[0] MOD 1000

;
;   Convert the given millisecond of day [orb_pos_time[1]] to second of day.
;
	secs = (DOUBLE(orb_pos_time[1]))/DOUBLE(1000.0)

;
;    Begin calculating the greenwich mean sidereal time **
;
	fday = secs/86400.00
	dj = DOUBLE(365*(year-1900)+(year-1901)/4+day-half)


;
;	      THE NEXT STATEMENT CAUSES THE REFERENCE EPOCH TO BE SHIFTED	
;	   TO THE J2000 REFERENCE EPOCH.
;

	T = (dj-DOUBLE(36525.0))/DOUBLE(36525.0)
        gst = DOUBLE(C0 + T*(C1 + T*(C2 + C4*T)) + C3*fday)
        gst = DOUBLE(gst MOD TWOPI)
        IF (gst LT DOUBLE(0.0)) THEN gst = gst + TWOPI

;
;   Convert gst to true of date by adjusting for nutation
;
	ic_get_nut_angles, T, deleps, delpsi, eps
	gst = gst + delpsi*cos(eps+deleps)

        END