C $Id: ic_gci_to_gse_xx.f,v 1.2 1998/07/24 21:54:30 asc Exp $
	subroutine IC_GCI_TO_GSE_xx (ORB_POS_TIME, 
     1                            TRANSFORM_MATRIX,DERIV_MATRIX,
     2                            VELOCITY_REQ,
     3                            STATUS)
C
	implicit none
C
C*  Calling parameters
C
	integer*4   ORB_POS_TIME(2) !time or orb vector, YYYYDDD, milli of day
	real*8	    TRANSFORM_MATRIX(3,3) !transformation matrix
	real*8      DERIV_MATRIX(3,3) !time derivative of transformation matrix
	integer*4   VELOCITY_REQ !flag to compute derivative matrix
	integer*4   STATUS  !Status of getting the transformation matrix
C
C*  Other variables
C
	real*8	    MEAN_MATRIX(3,3)	!mean of date transform matrix
	real*8	    CMATRIX(3,3)	!The conversion matrix to rotate from 
C                                !mean of date to true of date
C*  Other variables
C
 	integer*4   YEAR        !Year - YYYY
 	integer*4   DAY	        !day of year
	real*8      SECS        !Seconds
	real*8	    DXJUL	!Julian ephemeris date at beginning of year
	real*8	    FDAY	!Fraction of a day
	real*8      DJ0         !Julian Day of 2000
        real*8      T0          !TIME IN centuries from 2000
        real*8      UT          !Universal time in hours
        real*8      DTR         ! DEGREES TO RADIANS
C
	REAL*8  M, L, LDOT, DJX, EE, T00, CIRCULO
C
C*  START EXECUTABLE
C
C*  Convert the given millisecond of day [orb_pos_time(2)] to second of day.
C*  Convert the packed form into year and day-of-year
C
	SECS = (real(ORB_POS_TIME(2)))/1000.0D0
        YEAR = ORB_POS_TIME(1)/1000
	DAY  = MOD(ORB_POS_TIME(1),1000)
	FDAY = SECS/86400.0D0
	DJX  = DXJUL(YEAR) + real(DAY)
C
	DJ0 = DXJUL(2000) + real(1) + 0.5
	T0 = (DJX - DJ0)/36525.D0
        T00= (DJX + FDAY - DJ0) /36525.D0
	UT = FDAY*24.D0
C
	DTR = (4.D0*ATAN(1.D0)) / 180.D0
C
        EE = 23.4392911D0 - 0.0130042D0 * T00
C
	M  = 357.5254333D0 + 35999.04975D0*T0 
     1       - 0.000150D0*T0*T0 - 0.0000033D0*T0*T0*T0
     2       + 0.04107D0*UT
	M = CIRCULO (M)
C
	L  = 280.4659025D0 + 36000.76892D0*T0 
     1                             + 0.0003025D0*T0*T0
     2       + 0.04107D0*UT
	L = CIRCULO (L)
C
	L = L 
     1    + 0.00134D0*DCOS(DTR * (351.9841D0+ 22518.7541D0*T00) )         
     2    + 0.00154D0*DCOS(DTR * (254.0782D0+ 45037.5082D0*T00) )        
     3    + 0.00200D0*DCOS(DTR * (157.0477D0+ 32964.3577D0*T00) )        
     4    + 0.00179D0*DSIN(DTR * (297.8528D0+445267.1142D0*T00
     +                                -0.00144D0*T00*T00))
     5    + 0.00178D0*DSIN(DTR * (251.39D0  +    20.20D0*T00) )          
C
	L = CIRCULO (L)
C
	LDOT = L + (1.914657D0-0.004789D0*T00
     1      -0.000014D0*T00*T00)*DSIN(DTR*M)+ 
     2       (0.019994D0-0.000100*T00)*DSIN(DTR*2.*M)
     3       +0.000293D0*DSIN(DTR*3.*M)
	LDOT = CIRCULO (LDOT)
C
	MEAN_MATRIX(1,1) = DCOS(DTR*LDOT)
	MEAN_MATRIX(1,2) = DSIN(DTR*LDOT)
	MEAN_MATRIX(1,3) = 0.D0
	MEAN_MATRIX(2,1) = -MEAN_MATRIX(1,2)
	MEAN_MATRIX(2,2) = MEAN_MATRIX(1,1)
	MEAN_MATRIX(2,3) = 0.D0
	MEAN_MATRIX(3,1) = 0.D0
	MEAN_MATRIX(3,2) = 0.D0
	MEAN_MATRIX(3,3) = 1.D0
C
	CMATRIX(1,1) = 1.D0
	CMATRIX(1,2) = 0.D0
	CMATRIX(1,3) = 0.D0
	CMATRIX(2,1) = 0.D0
	CMATRIX(2,2) = DCOS(DTR*EE)
	CMATRIX(2,3) = DSIN(DTR*EE)
	CMATRIX(3,1) = 0.D0
	CMATRIX(3,2) = -CMATRIX(2,3)
	CMATRIX(3,3) = CMATRIX(2,2)
C
C*  multiply the mean_matrix and cmatrix to get the true of date transf_matrix
C
	CALL MULMAT (MEAN_MATRIX, CMATRIX, 3, 3, 3, TRANSFORM_MATRIX)
C
	RETURN
C
 	end










