C $Id: ic_gci_to_gse.f,v 1.2 1998/07/24 21:54:29 asc Exp $
CCCC
C
C  IC_GCI_TO_GSE - return a transformation matrix
C 
C  PURPOSE:  Calculate the transformation matrix between GCI and
C            GSE coordinate systems at a given date and time.
C
C  UNIT TYPE:  SUBROUTINE
C
C  INVOCATION METHOD:  CALL IC_GCI_TO_GSE (orb_pos_time, 
C                                          transform_matrix,
C	                                   deriv_matrix,
C                                          velocity_req,
C				           status)
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  TRANSFORM_MATRIX(3,3)  R*8    O    TRANSFORMATION MATRIX
C  DERIV_MATRIX(3,3)      R*8    O    TIME DERIVATIVE MATRIX
C  VELOCITY_REQ           I*4    I    FLAG TO SHOW WHETHER DERIV_MATRIX NEEDS
C                                      TO BE CALCULATED. IF = 0, THEN DERIV_
C                                      MATRIX IS NULL-FILLED ON RETURN.  IF
C                                      > 0, THEN THE DERIV_MATRIX WILL HAVE
C                                      NON-ZERO VALUES IN IT.
C  STATUS		  I*4	 O    STATUS OF GETTING TRANSFORMATION MATRIX
C
C  FILE/RECORD REFERENCES:  NONE
C
C  EXTERNAL VARIABLES:  NONE
C
C  EXTERNAL REFERENCES:
C    IC_POS_VEL_OF_SUN	 Calculates the position and velocity 
C			 of the sun in GCI coordinates
C    IC_CROSS_PRODUCT    Calculates the cross product of two vectors
C    MULMAT              Multiplies two matrices
C    IC_CONV_MATRIX	 Return the conversion matrix to go from mean of date
C                        to true of date
C
C  ABNORMAL TERMINATION CONDITIONS, ERROR MESSAGES:
C    If an error occurs, STATUS is returned as one of the following:
C       Error status returned by IC_POS_VEL_OF_SUN
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/21/90  Initial PDL
C  J. Lubelczyk                 B1R1      12/10/90  Coding
C  J. Lubelczyk                           09/11/91  Convert to true of date
C  J. LUBELCZYK                           09/19/91  Add error handling
C  J  LUBELCZYK                           11/4/91   Calc. obliquity using y,z
C                                                   comp. of sun pos.; length
C                                                   of all rows and columns of
C                                                   transform matrix = 1
C  J. LUBELCZYK ICCR #83, CCR #'S 130, 137 11/91    B3 update
C  C. RAYMOND   CCR 462         R 3.0     1/7/93    Added velocity 
C                                                    transformation capacity
C  B. SAMUELSON SPOF PORT    NONE         04/11/94  Change: required to port
C     (CSC)     (see notes)                          icss routines to SPOF
C
C  NOTES:
C
C     1.  The changes recorded under ID SPOF-PORT are required to make the
C         ICSS coordinate conversion routines, originally developed under
C         VAX-VMS 5.4 run on the UNIX-based workstations of the SPOF.  (Sun
C         SPARCstations and DEC DECstations).  The changes are as follows:
C            a.  Delete references to ICSS_INC
C            b.  Define Message texts or files to correspond to the messages
C                ICSS_SUCCESSFUL, etc. which are embedded in the error handling.
C            c.  Remove references to the NAG routines F01CRF and F01CKF (matrix
C                transposition and matrix multiplication routines).
C     2.  In addition, to successfully run the software packages, copies of the
C         Solar/Lunar/Planetary (SLP) file and timing coefficients file (TCC)
C         must be ported onto the SPOF.
CCCC
C
C  PDL:
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
C  if error getting the mean position of the sun, abort to 990
C
C  Convert the mean position of the sun to a unit vector
C
C  Calculate magnitude of Sun velocity vector
C
C  call IC_CONV_MATRIX to get the conversion matrix, CMATRIX, 
C    to rotate to true of date
C
C  call MULMAT to get the XAXIS by multiplying 
C    the CMATRIX and the mean of date XAXIS
C
C  Calculate the obliquity
C
C  Derive the Z axis from the obliquity
C                                                        _     _        _
C  call IC_CROSS_PRODUCT to execute the cross product of Z and X to get Y
C    ZxX = (X_ZAXIS, Y_ZAXIS, Z_ZAXIS)x(X_XAXIS, Y_XAXIS, Z_XAXIS) = Y
C
C  Fill the mean of date transformation matrix, MEAN_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
C  call MULMAT to multiply MEAN_MATRIX by CMATRIX to get the transform matrix
C
C  IF ( derivative matrix needed ) 
C  
C  THEN
C   
C     COMPUTE quantities needed for derivative matrix
C
C     BUILD derivative matrix
C
C  ENDIF
C
C990 CONTINUE
C
C  RETURN
C
CCCC
	subroutine IC_GCI_TO_GSE (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
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
        IF ( VELOCITY_REQ .GT. 0 ) THEN
	   MEAN_VEL(1)=PVSUN(4)
	   MEAN_VEL(2)=PVSUN(5)
	   MEAN_VEL(3)=PVSUN(6)
	   CALL MULMAT (CMATRIX, MEAN_VEL, 3, 3, 1, VXAXIS)
           PSLEN= VXAXIS(1)*XAXIS(1) + VXAXIS(2)*XAXIS(2)
     1             + VXAXIS(3)*XAXIS(3) 
	   S(1) = ( VXAXIS(1)/LEN) - XAXIS(1)*(PSLEN/LEN)
	   S(2) = ( VXAXIS(2)/LEN) - XAXIS(2)*(PSLEN/LEN)
	   S(3) = ( VXAXIS(3)/LEN) - XAXIS(3)*(PSLEN/LEN)
	   CALL IC_CROSS_PRODUCT(ZAXIS,S,YAXIS)
	   D_MATRIX(1,1)=S(1)
	   D_MATRIX(1,2)=S(2)
	   D_MATRIX(1,3)=S(3)
	   D_MATRIX(2,1)=YAXIS(1)
	   D_MATRIX(2,2)=YAXIS(2)
	   D_MATRIX(2,3)=YAXIS(3)
	   D_MATRIX(3,1)=0.0D0
	   D_MATRIX(3,2)=0.0D0
	   D_MATRIX(3,3)=0.0D0
	   CALL MULMAT (D_MATRIX, CMATRIX, 3, 3, 3, DERIV_MATRIX)
	ELSE
	   DO 100 K=1,3
	   DO 100 J=1,3
	      DERIV_MATRIX(J,K)=0.0d0

100	   CONTINUE
	ENDIF
990	CONTINUE
	RETURN
C
 	end
