C $Id: ic_gci_to_gsm.f,v 1.2 1998/07/24 21:54:30 asc Exp $
CCCC
C
C  IC_GCI_TO_GSM - return a transformation matrix
C 
C  PURPOSE:  Calculate the transformation matrix between GCI and
C            GSM coordinate systems at a given date and time.
C
C  UNIT TYPE:  SUBROUTINE
C
C  INVOCATION METHOD:  CALL IC_GCI_TO_GSM (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 
C                                      NEEDS TO BE CALCULATED. IF=0, THEN
C                                      DERIV_MATRIX IS ZERO-FILLED. IF > 0, 
C                                      THEN THE DERIV-MATRIX WILL HAVE NON-
C                                      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 of the sun in GCI coordinates
C    IC_CROSS_PRODUCT    Calculates the cross product of two vectors
C    IC_CONV_GEO_TO_GCI  Converts a vector in GEO coordinates to 
C                        a vector in GCI coordinates
C    MULMAT              Multiplies two matrices
C    IC_CONV_MATRIX      Returns the conversion matrix to go from mean of J2000
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/23/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 ICCR #83, CCR #'S 130, 137 11/91    B3 update
C  C. Raymond   CCR 568         B4/R2.1   02/92     Efficiency/accuracy 
C                                                    update; see note 2
C  C. Raymond   CCR 523         B4/R2.1   02/92     Change magnetic dipole
C                                                    to IGRF 1990
C  C. Raymond   CCR 462         R 3.0     1/7/93    Added velocity 
C                                                    transformation capacity.
C  B. SAMUELSON SPOF PORT       NONE      05/02/94  Change: required to port
C     (CSC)     (see notes)                          icss routines to SPOF
C  B. SAMUELSON UNIX PORT                 09/18/96  Updating unix version of 
C                                                    code. (see note 9) 
C
C  NOTES:
C
C  1. The position of the Earth's dipole is known once a given model is
C     adopted.  The current version uses the IGRF 1985.0 model.  
C
C  2. (CCR 568) The Sun vector should have been transformed from J2000
C     to true of date before computations were done.  This CCR corrected
C     that deficiency.
C
C  3. (CCR 523) Updated per October 1991 SWG to use the IGRF 1990.0 
C     magnetic reference model for the earth.
C
C  4. (CCR 462) The Sun unit velocity vector is computed as follows:
C       Compute Sun velocity vector magnitude PSLEN from appropriate
C        components of PVSUN.
C       Multiply Sun velocity by CMATRIX to convert to true of date
C        Sun velocity VXAXIS
C       Then form the vector S (LEN is the Sun position vector magnitude):   
C       
C	   S(1) = ( VXAXIS(1)/LEN) - XAXIS(1)*(PSLEN/LEN)
C	   S(2) = ( VXAXIS(2)/LEN) - XAXIS(2)*(PSLEN/LEN)
C	   S(3) = ( VXAXIS(3)/LEN) - XAXIS(3)*(PSLEN/LEN)
C
C  5. (CCR 462) The derivative of the Y-axis is computed from the following
C     expressions:
C       
C      VD = (YAXIS(1)*V_DOT(1) + YAXIS(2)*V_DOT(2) + YAXIS(3)*V_DOT(3))
C      Y_DOT(1) = (V_DOT(1)/NORMY) - YAXIS(1)*(VD/NORMY)
C      Y_DOT(2) = (V_DOT(2)/NORMY) - YAXIS(2)*(VD/NORMY)
C      Y_DOT(3) = (V_DOT(3)/NORMY) - YAXIS(3)*(VD/NORMY)
C
C  6. (CCR 462) The derivative of the Z-axis is computed as follows:
C      
C      CALL IC_CROSS_PRODUCT(S,YAXIS,CP1)
C      CALL IC_CROSS_PRODUCT(XAXIS,Y_DOT,CP2)
C      Vector add the two vectors CP1 and CP2
C
C  7.  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
C  8.  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
C  9.  Updated magnetic dipole from IGRF 1990 to IGRF 1995 and added a
C       call to icss_tilt_angle 
C
CCCC
C
C  PDL:
C
C  call IC_POS_VEL_OF_SUN to get the position of the sun [X axis in GSM] 
C    in the GCI coordinate system
C
C  if not successful getting the pos of the sun abort to 990
C
C  Call IC_CONV_MATRIX to get the CMATRIX to transform from J2000 (GCI)
C   to true of date
C
C  Convert the position of the sun to a unit vector
C
C  Call MULMAT to convert the unit Sun vector to true of date using the
C   CMATRIX
C
C  call IC_CONV_GEO_TO_GCI to convert the earth's dipole in the GEO
C    coordinate system to GCI
C
C  call IC_CROSS_PRODUCT to execute the cross product of the GCI dipole
C    and the position of the sun to get the Y axis. [ZxX=Y]
C
C  Normalize the YAXIS vector to get the required unit vector
C
C  call IC_CROSS_PRODUCT to execute the cross product of X AND Y TO GET Z
C
C  Fill the mean of date transformation matrix, MEAN_MATRIX.
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 ( VELOCITY_REQ > 0 ) 
C  THEN
C     COMPUTE Sun velocity unit vector S using algorithm in note 4
C     CALL IC_GCI_TO_GEO with velocity_request flag set to get derivative
C                                                         matrix GEOD_DOT
C     CALL MXTRP8 to transpose GEOD_DOT
C     CALL MULMAT to pre-multiply the vector GEO_DIPOLE by transposed
C                                 GEOD_DOT matrix forming D_PRIME_DOT
C     CALL IC_CROSS_PRODUCT to get the cross product of D_PRIME_DOT and
C                                                                 XAXIS
C     CALL IC_CROSS_PRODUCT to get the cross product of GCI_DIPOLE and S
C     CALCULATE V_DOT by adding the 2 cross products just computed
C     CALCULATE Y_DOT from the algorithm given in note 5.
C     CALCULATE Z_DOT from the algorithm given in note 6.
C     LOAD row 1 of D_MATRIX from S
C     LOAD row 2 of D_MATRIX from Y_DOT
C     LOAD row 3 of D_MATRIX from Z_DOT
C     CALL MULMAT to post-multiply D_MATRIX by CMATRIX to form DERIV_MATRIX
C  ELSE
C     LOAD DERIV_MATRIX with all zeroes
C  ENDIF
C
C990 Continue
C  RETURN
C
CCCC
	subroutine IC_GCI_TO_GSM (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) !time of orb vector, YYYYDDD, milli of day
	real*8	    TRANSFORM_MATRIX(3,3) !transformation matrix
	real*8      DERIV_MATRIX(3,3) !time derivative matrix
	integer*4   VELOCITY_REQ     !flag for deriv_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      VXAXIS(3)   !True of date Sun velocity 
        real*8      PSLEN       !Sun velocity vector magnitude
        real*8      S(3)        !Sun velocity unit vector
        real*8      GEOD_DOT(3,3) !GCI->GEO derivative matrix/transpose
        real*8      TEMP_MATRIX(3,3) !temporary result matrix
        real*8      DUMMY_MAT(3,3) !Dummy matrix
	real*8      V_DOT(3)    ! working array
	real*8      D_MATRIX(3,3) ! true of date derivativer matrix
	real*8      MEAN_VEL(3) ! Mean of date Sun velocity
	real*8      VD          ! scalar product of YAXIS and V_DOT
        real*8      CP1(3)      !working array
        real*8      CP2(3)      !working array
        real*8      D_PRIME_DOT(3) !working array
        real*8      Y_DOT(3)    ! time derivative of Y-axis
        real*8      Z_DOT(3)    ! time derivative of Z-axis
        real*8      MEAN_XAXIS(3) !used in computing Sun position
	real*8	    YAXIS(3)    !GSE's Y axis in GCI coordinates
	real*8	    ZAXIS(3)    !GSE's Z axis in GCI coordinates
	real*8	    GEO_DIPOLE(3)  !Position of earth's dipole in GEO coord.
	real*8	    GCI_DIPOLE(3)  !Position of earth's dipole in GCI coord.
	real*8	    LEN		!magnitude of the pos of sun vector
	real*8	    NORMY       !Norm of the Y axis
	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
        real*8      tilt
	integer*4   J,K         ! do-loop indices
C
C
C*  Begin executable code
C
C*  The GSM Z axis is in the plane of the dipole and the direction of the sun.
C*  Thus as the earth spins, the Z axis will track the 'wobble' of the dipole.
C*  The position of the earth's dipole as of 1995 (IGRF) 
C*   in the GEO system is the following:
C
	GEO_DIPOLE(1) =  0.57273114D-1    ! CCR 2142 Was 0.6105036340D-1
	GEO_DIPOLE(2) = -0.17362989D0     !          Was -0.178467593D0
	GEO_DIPOLE(3) = 0.983144166D0     !          Was 0.9820499840D0

        call icss_tilt_angle (orb_pos_time, tilt, geo_dipole, status)
C
C*  call IC_POS_VEL_OF_SUN to get the position of the sun [X axis in GSM]
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 position of the sun (XAXIS) 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
C*  Convert the earth's magnetic dipole in GEO coordinates to GCI coordinates
C
	call IC_CONV_GEO_TO_GCI (ORB_POS_TIME, GEO_DIPOLE, GCI_DIPOLE)
C
C*  call IC_CROSS_PRODUCT to get the direction 
C*   of GSM's Y axis in the GCI coordinate system
C*   The result is not a unit vector due to the "wobble" of the earth's dipole
C*   about its axis of rotation.
C
	call IC_CROSS_PRODUCT (GCI_DIPOLE, XAXIS, YAXIS)
C
C*  Normalize the YAXIS vector to get a unit vector
C
	NORMY = DSQRT((YAXIS(1)*YAXIS(1)) + (YAXIS(2)*YAXIS(2))
     -                + (YAXIS(3)*YAXIS(3)))
	YAXIS(1) = YAXIS(1)/NORMY
	YAXIS(2) = YAXIS(2)/NORMY
	YAXIS(3) = YAXIS(3)/NORMY
C
C *EXECUTE THE CROSS PRODUCT OF X AND Y TO GET GSM's Z axis in GCI coordinates
C
	call IC_CROSS_PRODUCT (XAXIS, YAXIS, ZAXIS)
C	
C*  Fill the mean of date transformation matrix.
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*  Get the conversion matrix to rotate to true of date
C
	CALL IC_CONV_MATRIX (ORB_POS_TIME, CMATRIX)
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)
        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_GCI_TO_GEO(ORB_POS_TIME,DUMMY_MAT,GEOD_DOT,
     1                        velocity_req)
	   CALL MXTRP8 (GEOD_DOT,TEMP_MATRIX,3,3)
           GEOD_DOT(1,1) = TEMP_MATRIX(1,1)
           GEOD_DOT(2,1) = TEMP_MATRIX(2,1)
           GEOD_DOT(3,1) = TEMP_MATRIX(3,1)
           GEOD_DOT(1,2) = TEMP_MATRIX(1,2)
           GEOD_DOT(2,2) = TEMP_MATRIX(2,2)
           GEOD_DOT(3,2) = TEMP_MATRIX(3,2)
           GEOD_DOT(1,3) = TEMP_MATRIX(1,3)
           GEOD_DOT(2,3) = TEMP_MATRIX(2,3)
           GEOD_DOT(3,3) = TEMP_MATRIX(3,3)
	   CALL MULMAT (GEOD_DOT,GEO_DIPOLE,3,3,1,D_PRIME_DOT)
           CALL IC_CROSS_PRODUCT(D_PRIME_DOT,XAXIS,CP1)
	   CALL IC_CROSS_PRODUCT(GCI_DIPOLE,S,CP2)
           V_DOT(1)=CP1(1)+CP2(1)
           V_DOT(2)=CP1(2)+CP2(2)
           V_DOT(3)=CP1(3)+CP2(3)
	   VD = (YAXIS(1)*V_DOT(1) + 
	1        YAXIS(2)*V_DOT(2) + YAXIS(3)*V_DOT(3))
	   Y_DOT(1) = (V_DOT(1)/NORMY) - YAXIS(1)*(VD/NORMY)
	   Y_DOT(2) = (V_DOT(2)/NORMY) - YAXIS(2)*(VD/NORMY)
	   Y_DOT(3) = (V_DOT(3)/NORMY) - YAXIS(3)*(VD/NORMY)
	   CALL IC_CROSS_PRODUCT(S,YAXIS,CP1)
	   CALL IC_CROSS_PRODUCT(XAXIS,Y_DOT,CP2)
	   Z_DOT(1)=CP1(1)+CP2(1)
	   Z_DOT(2)=CP1(2)+CP2(2)
	   Z_DOT(3)=CP1(3)+CP2(3)
	   D_MATRIX(1,1)=S(1)
	   D_MATRIX(1,2)=S(2)
	   D_MATRIX(1,3)=S(3)
	   D_MATRIX(2,1)=Y_DOT(1)
	   D_MATRIX(2,2)=Y_DOT(2)
	   D_MATRIX(2,3)=Y_DOT(3)
	   D_MATRIX(3,1)=Z_DOT(1)
	   D_MATRIX(3,2)=Z_DOT(2)
	   D_MATRIX(3,3)=Z_DOT(3)
	   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
C
990	CONTINUE
	RETURN
C
	END
