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 = (dflotj(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) + DFLOTJ(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