C $Id: ic_gsm_to_geo.f,v 1.2 1998/07/24 21:54:35 asc Exp $
CCCC
C
C  IC_GSM_TO_GEO - return a transformation matrix
C 
C  PURPOSE:  Calculate the transformation matrix from GSM
C            coordinates to GEO coordinates at a given date and time.
C
C  UNIT TYPE:  SUBROUTINE
C
C  INVOCATION METHOD:  CALL IC_GSM_TO_GEO (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 OF TRANSFORMATION 
C                                      MATRIX
C  VELOCITY_REQ           I*4    I    FLAG TO SHOW WHETHER DERIV_MATRIX IS
C                                      COMPUTED.  IF > 0, DERIV_MATRIX IS
C                                      NON-ZERO ON RETURN.
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_GCI_TO_GEO    Returns the GCI to GEO transformation matrix
C    IC_GCI_TO_GSM    Returns the GCI to GSM transformation matrix
C    MXTRP8           Routine that return the transpose of a matrix
C    MULMAT           Routine that multiplies two matrices
C
C  ABNORMAL TERMINATION CONDITIONS, ERROR MESSAGES:
C    If an error occurs, STATUS is returned as one of the following:
C     Returned from IC_GCI_TO_GSM
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/11/90  CODING
C  J. LUBELCZYK                           09/19/91  Added Error handling
C  J. LUBELCZYK ICCR #83, CCR #'S 130, 137 11/91    B3 update
C  C. RAYMOND   CCR 462                   07/17/93  Added code necessary
C                                                    for velocity 
C                                                    transformations
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.  (CCR 462) Any quantity which transforms like velocity in a 
C      non-inertial coordinate frame mixes the input velocity with
C      components from the position vector as well.  In this case,
C      since the derivative matrix is computed by multiplying out
C      matrices from two intermediate calls, the expressions are
C      complicated.  This note gives some details.
C
C      If the matrices returned from a call to IC_GCI_TO_GEO are
C      denoted G1, G1';  and the matrices returned from a call to
C      IC_GCI_TO_GSM are denoted G2, G2';  then the derivative
C      matrix is :
C                            T         T        T
C      DERIV_MATRIX = G1'*(G2) - G1*(G2)*G2'*(G2)
C               T 
C      where (G2) means the transpose of matrix G2.
C
C  2.  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  3.  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.
C
CCCC
C
C  PDL:
C
C  call IC_GCI_TO_GEO to get the GCI to GEO transformation matrix
C
C  call IC_GCI_TO_GSM to get the GCI to GSM transformation matrix
C  if (error getting gci to gsm transfromation matrix) abort to 990
C
C  call MXTRP8 routine to invert the GCI_GSM_MATRIX to GSM_GCI_MATRIX
C
C  call MULMAT to calculate the gsm_to_geo TRANSFORM_MATRIX by multiplying
C       GCI_GEO_MATRIX and GSM_GCI_MATRIX.
C
C  IF ( VELOCITY_REQ .GT. 0 ) 
C
C  THEN
C
C     COMPUTE DERIV_MATRIX
C
C  ELSE
C
C     ZERO DERIV_MATRIX
C
C  ENDIF
C
C990 Continue
C  RETURN
C
CCCC
	subroutine IC_GSM_TO_GEO (ORB_POS_TIME,
     1                            TRANSFORM_MATRIX,
     2                            DERIV_MATRIX,
     3                            VELOCITY_REQ,
     4                            STATUS)
C
	implicit none
C
C*  Calling parameters
C
	integer*4   ORB_POS_TIME(2)	!ISTP TIME FORMAT (YYYYDDD,milli of day)
	real*8	    TRANSFORM_MATRIX(3,3) !transformation matrix
	real*8      DERIV_MATRIX(3,3)   !time derivative of TRANSFORM_MATRIX
	integer*4   VELOCITY_REQ        !flag to compute DERIV_MATRIX
	integer*4   STATUS  !Status of getting transformation matrix
C
C*  Other variables
C
	real*8	    GCI_TO_GEO_MATRIX(3,3)  !gci to geo transformation matrix
	real*8	    GCI_GSM_MATRIX(3,3) !GCI to GSM, GSM to GCI transf matrix
        real*8      TEMP_MATRIX(3,3)  !temporary result matrix
	real*8      I1(3,3)             !Intermediate array
	real*8      I2(3,3)             !Intermediate array
	real*8      I3(3,3)             !Intermediate array
	real*8      G1_prime(3,3)       !returned time derivative from 
c	                                ! IC_GCI_TO_GEO
	real*8      G2_prime(3,3)       !returned time derivative from 
c	                                ! IC_GCI_TO_GSM
	integer*4   I,J                 !Do-loop indices
C
C
C*  Start executable code
C
C*  call IC_GCI_TO_GEO to get the GCI to GEO transformation matrix
C
	call IC_GCI_TO_GEO (ORB_POS_TIME, GCI_TO_GEO_MATRIX,G1_PRIME,
	1                   VELOCITY_REQ)
C
C*  call IC_GCI_TO_GSM to get the GCI to GSM transformation matrix
C
	call IC_GCI_TO_GSM (ORB_POS_TIME, GCI_GSM_MATRIX, G2_PRIME,
	1                   VELOCITY_REQ, STATUS)
	if (STATUS .ne. 0) goto 990
C
C*  call MXTRP8 routine to invert the GCI_GSM_MATRIX to GSM_GCI_MATRIX
C
	call MXTRP8 (GCI_GSM_MATRIX, TEMP_MATRIX, 3, 3)
        GCI_GSM_MATRIX(1,1) = TEMP_MATRIX(1,1)
        GCI_GSM_MATRIX(2,1) = TEMP_MATRIX(2,1)
        GCI_GSM_MATRIX(3,1) = TEMP_MATRIX(3,1)
        GCI_GSM_MATRIX(1,2) = TEMP_MATRIX(1,2)
        GCI_GSM_MATRIX(2,2) = TEMP_MATRIX(2,2)
        GCI_GSM_MATRIX(3,2) = TEMP_MATRIX(3,2)
        GCI_GSM_MATRIX(1,3) = TEMP_MATRIX(1,3)
        GCI_GSM_MATRIX(2,3) = TEMP_MATRIX(2,3)
        GCI_GSM_MATRIX(3,3) = TEMP_MATRIX(3,3)
C
C*  call MULMAT to calculate the gsm_to_geo TRANSFORM_MATRIX by multiplying
C*   GCI_GEO_MATRIX and GSM_GCI_MATRIX.
C
	call MULMAT (GCI_TO_GEO_MATRIX, GCI_GSM_MATRIX, 
     1               3, 3, 3, TRANSFORM_MATRIX)
C
C*  if velocity_req gt 0 --
C*  then
C
	if ( VELOCITY_REQ .GT. 0 ) THEN
	   CALL MULMAT(TRANSFORM_MATRIX, G2_PRIME, 3, 3, 3, I1)
	   CALL MULMAT(I1, GCI_GSM_MATRIX, 3, 3, 3, I2) 
	   CALL MULMAT(G1_PRIME, GCI_GSM_MATRIX, 3, 3, 3, I3)
	   do 100 i=1,3
	   do 100 j=1,3
	      DERIV_MATRIX(I,J) = I3(i,j) - I2(i,j)
100	   continue
C
C*   else
C*      zero-fill DERIV_MATRIX
C
	else
	   do 200 i=1,3
	   do 200 j=1,3
	      DERIV_MATRIX(i,j) = 0.0d0
200	   continue
	endif
990	continue
	return
C
	end
