C $Id: ic_hseb_to_hsea.f,v 1.2 1998/07/24 21:54:37 asc Exp $
	subroutine IC_HSEB_TO_HSEA (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	    PVSUN(6)	!Position and velocity vectors of the sun
	real*8	    XAXIS(3) 	!pos of sun, GSE's X axis in GCI coordinates
        real*8      YAXIS(3)    !GSE's Y axis in GCI coordinates
        real*8      ZAXIS(3)    !GSE's Z axis in GCI coordinates
	real*8      VXAXIS(3)   !working array
	real*8      MEAN_VEL(3) !working array
	real*8      D_MATRIX(3,3) !working matrix
	real*8      PSLEN       !Dot product Sun velocity and position
	real*8      S(3)        !Time derivative of Sun unit vector
	real*8	    LEN		!magnitude of the pos of sun vector
	real*8	    OBLIQUITY	!Diff between the equatorial
C                               !and ecliptic planes
	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
	integer*4   J,K         !loop indices
        real*8      MEAN_XAXIS(3) !Mean of date xaxis [pos of sun,unit vector]
C
C*  START EXECUTABLE
C
C*  call IC_POS_VEL_OF_SUN to get the position of the sun [X axis in GSE] 
C*    in the GCI coordinate system
C
	call IC_POS_VEL_OF_SUN (ORB_POS_TIME, PVSUN, STATUS)
	if (STATUS .ne. 0) goto 990
C
C*  Convert the mean of date position of the sun 
C*  to a unit vector pointing to the sun
C
	len = dsqrt(PVSUN(1)*PVSUN(1) + PVSUN(2)*PVSUN(2) 
     1              + PVSUN(3)*PVSUN(3))
	MEAN_XAXIS(1) = PVSUN(1)/LEN
	MEAN_XAXIS(2) = PVSUN(2)/LEN
	MEAN_XAXIS(3) = PVSUN(3)/LEN
C
C*  Get the conversion matrix to rotate to true of date
C
	CALL IC_CONV_MATRIX (ORB_POS_TIME, CMATRIX)
C
        call MULMAT (CMATRIX, MEAN_XAXIS, 3, 3, 1, XAXIS)
C
	obliquity = datan(XAXIS(3)/XAXIS(2))
C
C* The GSE Z axis (the direction of the ecliptic pole)
C* is the following in the GCI system.
C
	ZAXIS(1) = 0.0d0
	ZAXIS(2) = -(dsin(OBLIQUITY))
	ZAXIS(3) = dcos(OBLIQUITY)
C
C*  EXECUTE THE CROSS PRODUCT OF Z AND X TO GET Y
C
	call IC_CROSS_PRODUCT (ZAXIS, XAXIS, YAXIS)
C
C*  Fill the mean of date transformation matrix with
C*   the previously calculated X, Y, and Z axes
C*   Row 1 = X axis, Row 2 = Y axis, Row 3 = Z axis
C
	MEAN_MATRIX(1,1) = XAXIS(1)
	MEAN_MATRIX(1,2) = XAXIS(2)
	MEAN_MATRIX(1,3) = XAXIS(3)
	MEAN_MATRIX(2,1) = YAXIS(1)
	MEAN_MATRIX(2,2) = YAXIS(2)
	MEAN_MATRIX(2,3) = YAXIS(3)
	MEAN_MATRIX(3,1) = ZAXIS(1)
	MEAN_MATRIX(3,2) = ZAXIS(2)
	MEAN_MATRIX(3,3) = ZAXIS(3)
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
	CMATRIX(1,1) = 1.D0
	CMATRIX(1,2) = 0.D0
	CMATRIX(1,3) = 0.D0
	CMATRIX(2,1) = 0.D0
	CMATRIX(2,2) = DCOS(OBLIQUITY)
	CMATRIX(2,3) = DSIN(OBLIQUITY)
	CMATRIX(3,1) = 0.D0
	CMATRIX(3,2) = -CMATRIX(2,3)
	CMATRIX(3,3) =  CMATRIX(2,2)
C
        CALL GAUSSJ(CMATRIX,3,3,MEAN_MATRIX,3,3)
C
	CALL MULMAT (TRANSFORM_MATRIX,CMATRIX, 3, 3, 3, MEAN_MATRIX)
C
	TRANSFORM_MATRIX(1,1)=-MEAN_MATRIX(1,1)
	TRANSFORM_MATRIX(1,2)=-MEAN_MATRIX(1,2)
	TRANSFORM_MATRIX(1,3)= MEAN_MATRIX(1,3)
	TRANSFORM_MATRIX(2,1)=-MEAN_MATRIX(2,1)
	TRANSFORM_MATRIX(2,2)=-MEAN_MATRIX(2,2)
	TRANSFORM_MATRIX(2,3)= MEAN_MATRIX(2,3)
	TRANSFORM_MATRIX(3,1)= MEAN_MATRIX(3,1)
	TRANSFORM_MATRIX(3,2)= MEAN_MATRIX(3,2)
	TRANSFORM_MATRIX(3,3)= MEAN_MATRIX(3,3)

	CMATRIX(1,1) = 0.D0
	CMATRIX(1,2) = 0.D0
	CMATRIX(1,3) = 0.D0
	CMATRIX(2,1) = 0.D0
	CMATRIX(2,2) = 0.D0
	CMATRIX(2,3) = 0.D0
	CMATRIX(3,1) = 0.D0
	CMATRIX(3,2) = 0.D0
	CMATRIX(3,3) = 0.D0
C
        CALL GAUSSJ(TRANSFORM_MATRIX,3,3,CMATRIX,3,3)
C
c
c derivative matrix not developed yet
c
	   DO 100 K=1,3
	   DO 100 J=1,3
	      DERIV_MATRIX(J,K)=0.0d0
100	   CONTINUE
990	CONTINUE
	RETURN
C
 	end


















