C $Id: ic_conv_matrix.f,v 1.2 1998/07/24 21:54:27 asc Exp $
CCCC
C
C  IC_CONV_MATRIX - Returns the conversion matrix that is necessary to rotate
C		    from mean of date to true of date
C 
C  PURPOSE:  THIS SUBROUTINE CALCULATES, THROUGH APPROPRIATE ANALYTIC
C            EXPRESSIONS, VALUES FOR THE PRECESSION AND NUTATION
C            ANGLES AND THE MATRIX REQUIRED TO ROTATE FROM MEAN OF
C            JULIAN 2000 TO TRUE OF DATE.
C
C  UNIT TYPE:  SUBROUTINE
C
C  INVOCATION METHOD:  CALL IC_CONV_MATRIX (orb_pos_time,
C					    cmatrix)
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  CMATRIX(3,3)        R*8    O    Matrix to rotate from J2000 to true of date
C
C  FILE/RECORD REFERENCES:  NONE
C  
C  EXTERNAL VARIABLES:  NONE
C
C  EXTERNAL REFERENCES:
C	MULMAT - routine that multiplies two matrices
C       IC_GET_NUT_ANGLES          Routine to compute the nutation angles 
C                             
C
C  ABNORMAL TERMINATION CONDITIONS, ERROR MESSAGES:  NONE
C
C  ASSUMPTIONS, CONSTRAINTS, RESTRICTIONS:  NONE
C
C  DEVELOPMENT HISTORY
C
C  AUTHOR	CHANGE ID	RELEASE	  DATE	  DESCRIPTION OF CHANGE
C  ------	---------	-------   ----	  ---------------------
C  C. Raymond                                     Original Algorithm, Coding
C  J. Lubelczyk					  Prolog, Coding
C  J. LUBELCZYK ICCR #83, CCR #'S 130,137 11/91   B3 Update
C  C. Raymond   CCR 568         B4/R2.1   02/92   Code efficiency update; 
C                                                  see note 1.
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 568) This routine had the computationally expensive nutation
C     series embedded in it.  This CCR replaced that computation with a
C     call to the routine IC_GET_NUT_ANGLES which has the identical 
C     nutation series.  However, by changing IC_GET_NUT_ANGLES to only
C     recompute every 1/2 day, substantial CPU time was saved with a
C     minimal cost in accuracy.
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  Calculate the julian date, time in julian centuries from J2000
C
C  Calculate the precession angles
C
C  Compute the transformation matrix between the mean equator and
C   the equinox of 1950 and the mean equator and equinox of date.
C
C  Calculate the true obliquity of date
C
C  Call IC_GET_NUT_ANGLES to calculate the nutation in obliquity
C   and longitude
C
C  Calculate the elements of the nutation matrix
C
C  Call MULMAT to multiply the nutation and precession matrices, the result
C    is the conversion matrix, CMATRIX
C
C  return
C
CCCCC
      SUBROUTINE IC_CONV_MATRIX(ORB_POS_TIME, CMATRIX)
C
        IMPLICIT    NONE
C
C*  Calling Parameters
C
	integer*4   ORB_POS_TIME(2) !ISTP time format (YYYYDDD, Milli of Day)
	REAL*8	    CMATRIX(3,3)    !Conversion matrix
C
C*  Other Variables
C
	REAL*8	    TIME            !Time in Julian Centuries of 36525.0 
C			            !mean solar days from J2000.
	real*8      SECS	    !Seconds
 	integer*4   YEAR	    !Year - YYYY
 	integer*4   DAY		    !day of year
	REAL*8	    NUTMAT(3,3)	    !Nutation matrix
	REAL*8	    PREMAT(3,3)	    !Precession matrix
	REAL*8	    T2		    !TIME squared
	REAL*8	    T3		    !TIME cubed
	REAL*8	    PI		    !3.14159...
	REAL*8	    DTR		    !Conversion factor, Degree to Radians
	REAL*8	    STR		    !Conversion factor, Seconds to Radians
	REAL*8	    R		    !
	REAL*8	    EPSO	    !Mean obliquity of date
	REAL*8	    ZETA	    !Precession angle
	REAL*8	    THETA	    !Precession angle
	REAL*8	    ZEE		    !Precession angle
	REAL*8	    SINZET	    !Sine of ZETA
	REAL*8	    COSZET	    !Cosine of ZETA
	REAL*8	    SINZEE	    !Sine of ZEE
	REAL*8	    COSZEE	    !Cosine of ZEE
	REAL*8	    SINTHE	    !Sine of THETA
	REAL*8	    COSTHE	    !Cosine of THETA
	REAL*8	    DELEPS	    !DELTA EPSILON, Nutation in obliquity
	REAL*8	    DELPSI	    !DELTA PSI, Nutation in longitude
	REAL*8	    EPS		    !EPSILON
	REAL*8	    COSEP	    !Cosine of EPS
	REAL*8	    COSEPO	    !Cosine of OBLIQUITY
	REAL*8	    COSPSI	    !Cosine of PSI
	REAL*8	    SINEP	    !Sine of EPS
	REAL*8	    SINEPO	    !Sine of OBLIQUITY
	REAL*8	    SINPSI	    !Sine of PSI
	real*8	    DXJUL	    !Julian ephemeris date at beginning of year
	real*8	    JDJ2000	    !Julian day of J2000
	real*8	    JUL_DAY	    !Julian date in astronomer's terms
	integer*4   I		    !Year variable for DXJUL 
	real*8	    FDAY	    !Fraction of day
C
	PARAMETER   (JDJ2000 = 2451545.0D0)
C
C       STATEMENT FUNCTION DEFINITION FOR DXJUL -- JULIAN EPHEMERIS
C       DATE AT BEGINNING OF YEAR.
C
	DXJUL(I) = (-32075+1461*(I+4800-13/12)/4
     1             +367*(-1+13/12*12)/12-3
     2             *((I+4900-13/12)/100)/4)-0.5d0
C
      DATA R/1296000.0D0/
C
C*  START EXECUTABLE STATEMENTS
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)
C
C*  Calculate the julian date and the time in Julian centuries from J2000
C
      FDAY = SECS/86400.00D0
      JUL_DAY = DXJUL(YEAR) + real(DAY)+FDAY
      TIME = (JUL_DAY - JDJ2000)/36525.0D0
C
      T2 = TIME*TIME
      T3 = TIME*T2
C
C     CALCULATE CONVERSION FACTORS: DEGREES TO RADANS (DTR), SECONDS TO 
C     RADIANS (STR)                                                     
C                                                                       
      PI=4.0D0*DATAN(1.D0)                                              
      DTR=PI/180.D0                                                     
      STR=DTR/3600.D0                                                   
C
C     CALCULATE PRECESSION ANGLES
C
      ZETA   = 0.11180860865024398D-01*TIME
     1       + 0.14635555405334670D-05*T2
     2       + 0.87256766326094274D-07*T3
      THETA  = 0.97171734551696701D-02*TIME
     1       - 0.20684575704538352D-05*T2
     2       - 0.20281210721855218D-06*T3
      ZEE    = 0.11180860865024398D-01*TIME
     1       + 0.53071584043698687D-05*T2
     2       + 0.88250634372368822D-07*T3
      SINZET = DSIN(ZETA)                                               
      COSZET = DCOS(ZETA)                                               
      SINZEE = DSIN(ZEE)   
      COSZEE = DCOS(ZEE)   
      SINTHE = DSIN(THETA) 
      COSTHE = DCOS(THETA) 
C                                                                       
C     COMPUTE THE TRANSFORMATION MATRIX BETWEEN MEAN EQUATOR AND        
C     EQUINOX OF 1950 AND MEAN EQUATOR AND EQUINOX OF DATE. THIS        
C     MATRIX IS CALLED PREMAT.                                          
C                                                                       
      PREMAT(1,1) = -SINZET*SINZEE  + COSZET*COSZEE*COSTHE              
      PREMAT(2,1) =  COSZEE*SINZET  + SINZEE*COSTHE*COSZET              
      PREMAT(3,1) =  SINTHE*COSZET                                      
      PREMAT(1,2) = -SINZEE*COSZET  - COSZEE*COSTHE*SINZET              
      PREMAT(2,2) =  COSZEE*COSZET  - SINZEE*COSTHE*SINZET              
      PREMAT(3,2) = -SINTHE*SINZET                                      
      PREMAT(1,3) = -COSZEE*SINTHE                                      
      PREMAT(2,3) = -SINZEE*SINTHE                                      
      PREMAT(3,3) =  COSTHE                                             
C
C     CALCULATE MEAN OBLIQUITY OF DATE (EPSO). WHERE TIME IS MEASURED IN
C     JULIAN CENTURIES FROM 2000.0.
C
      EPSO=(1.813D-3*T3-5.9D-4*T2
     *     -4.6815D+1*TIME+8.4381448D+4)*STR
C                                                                       
C     CALL IC_GET_NUT_ANGLES TO COMPUTE NUTATION IN OBLIQUITY AND LONGITUDE
C                                                                       
      CALL IC_GET_NUT_ANGLES(TIME,DELEPS,DELPSI,EPS)
C                                                   
      COSEP=DCOS(EPS)                               
      COSEPO=DCOS(EPSO)                             
      COSPSI=DCOS(DELPSI)                           
      SINEP=DSIN(EPS)                               
      SINEPO=DSIN(EPSO)                             
      SINPSI=DSIN(DELPSI)                           
      NUTMAT(1,1)=COSPSI                            
      NUTMAT(1,2)=-SINPSI*COSEPO                    
      NUTMAT(1,3)=-SINPSI*SINEPO                    
      NUTMAT(2,1)=SINPSI*COSEP                      
      NUTMAT(2,2)=COSPSI*COSEP*COSEPO+SINEP*SINEPO  
      NUTMAT(2,3)=COSPSI*COSEP*SINEPO-SINEP*COSEPO  
      NUTMAT(3,1)=SINPSI*SINEP                      
      NUTMAT(3,2)=COSPSI*SINEP*COSEPO-COSEP*SINEPO  
      NUTMAT(3,3)=COSPSI*SINEP*SINEPO+COSEP*COSEPO  
C                                                                
C     CALCULATE ELEMENTS OF NUTMAT * PREMAT.  THIS MATRIX IS THE 
C     ANALYTICALLY CALCULATED TRANSFORMATION MATRIX, WHICH WILL  
C     TRANSFORM THE MEAN EARTH EQUATOR AND EQUINOX OF J2000 INTO
C     THE TRUE EARTH EQUATOR AND EXQUINOX OF DATE.
C                                                                
      CALL  MULMAT (NUTMAT, PREMAT, 3, 3, 3, CMATRIX)
C
      RETURN
C
      END








