SUBROUTINE DIRNA IMPLICIT NONE C C Add to TOC variables of new equation of equinox procedure C C PROGRAMMER: C 93.09.01 Norbert Zacharias - Equation of the equinoxes GAST update C 94.01.04 David Gordon - Cleaned up TOC entries. C 98.01.20 David Gordon - New TOC entries for New definition of the C equation of the equinoxes. Removed Lcodes C for previous EQE (calc 8.x) contributions. C INCLUDE 'ccon.i' C VARIABLES 'FROM': C 1. KDIUC - THE DIURNAL SPIN UTILITY ROUTINE FLOW CONTROL FLAG. C = 0 (default) Use new definition of the equation of C the equinoxes. C = 1 Use old definition of the equation of the C equinoxes. C C Lcode for difference in GAST (UT1) C [Even though we change the default to use the New definition, we keep this C Lcode the same as before so we don't have to change the name.] CALL ADDR (2,'EQE DIFF','GAST diff. (New-Old) from EQE ', * 1,1,1) C C Add contribution Lcode for the delay and rate contribution due to the C change in the definition of the equation of the equinox. The default is C now to use the New definition. Also remove previous EQE contribution C Lcodes that might be in the database. CALL DELR(2,'EQE CONT') C For the default (New definition), put in Lcode to convert back to the old C definition. If(KDIUC .eq. 0) Then CALL ADDR (2,'OLDEQCON','Old Eqn. of Equinoxes Contrib. ', * 2,1,1) CALL DELR(2,'NEWEQCON') Endif C C For the Old definition, put in Lcode to convert to the new definition. If(KDIUC .eq. 1) Then CALL ADDR (2,'NEWEQCON','New Eqn. of Equinoxes Contrib. ', * 2,1,1) CALL DELR(2,'OLDEQCON') Endif C Return END C C*************************************************************************** SUBROUTINE DIRNL ( DATDCT, DPSI, DUT1AT, EPS, DEPS, FA, UT1, 1 XJD, DIURNV, GAST, GMST, GASTD, RS ) IMPLICIT None C C 1. DIRNL C C 1.1 DIRNL PROGRAM SPECIFICATION C C 1.1.1 DIRNL is the utility routine which computes the diurnal spin portion of C the complete crust fixed to J2000.0 rotation matrix and its first two CT C time derivatives. DIRNL also computes the diurnal angular velocity of C the Earth, the Greenwich Mean Sidereal Time, and the Greenwich Apparent C Siderial Time (GAST) and its CT time derivative. C C 1.1.2 RESTRICTIONS - NONE C C 1.1.3 REFERENCES - 1) 'THE EXPLANATORY SUPPLEMENT TO THE AMERICAN C EPHEMERIS AND NAUTICAL ALMANAC", P.72-76, C 2) AOKI, S. ET AL., "THE NEW DEFINITION OF UNIVERSAL C TIME", ASTRON. ASTROPHYS., ????,1980. C 3) McCarthy, D., IERS Technical Note 13, Paris 1992 C C 1.2 DIRNL PROGRAM INTERFACE C C 1.2.1 CALLING SEQUENCE - C C INPUT VARIABLES: C 1. DATDCT - THE PARTIAL DERIVATIVE OF ATOMIC TIME WITH C RESPECT TO COORDINATE TIME. (SEC/SEC) C 2. DPSI(2) - THE NUTATION IN LONGITUDE AND ITS CT C TIME DERIVATIVE. (RAD, RAD/SEC) C 3. DUT1AT - THE PARTIAL DERIVATIVE OF UT1 TIME WITH C RESPECT TO ATOMIC TIME. (SEC/SEC) C 4. EPS(2) - THE TRUE OBLIQUITY OF THE ECLIPTIC AND ITS C CT TIME DERIVATIVE. (RAD, RAD/SEC) C 5. DEPS(2) - THE NUTATION IN OBLIQUITY AND ITS CT TIME C DERIVATIVE COMPUTED FROM WAHR OR TAKEN FROM C THE DATA BASE. (RAD, RAD/SEC) C 6. FA(5) - The fundamental nutation arguments (see NUTFA) C 6. OMEGA - The mean longitude of the ascending node of the C Moon. (arcsec!) C 7. UT1 - THE UT1 TIME OF THE DAY. (SEC) C 8. XJD - THE JULIAN DATE AT ZERO HOURS UTC OF THE C DATE IN QUESTION. (DAYS) C C OUTPUT VARIABLES: C 1. DIURNV - THE DIURNAL ANGULAR VELOCITY OF THE EARTH. C (RAD/SEC) C 2. GAST(2) - THE GREENWICH APPARENT SIDEREAL TIME AND C ITS CT TIME DERIVATIVE. (RAD, RAD/SEC) C 3. GMST - THE GREENWICH MEAN SIDEREAL TIME. (RAD) C 4. GASTD - The difference in GAST (new-old version) C 5. RS(3,3,3) - THE DIURNAL SPIN PORTION OF THE COMPLETE C CRUST FIXED TO 2000.0 ROTATION MATRIX C AND THE FIRST TWO CT TIME DERIVATIVES OF C THAT MATRIX. (UNITLESS, 1/SEC, 1/SEC**2) C Real*8 DATDCT, DPSI(2), DUT1AT, EPS(2), DEPS(2), FA(5), UT1, . XJD, DIURNV, GAST(2), GMST, GASTD, RS(3,3,3) C C 1.2.2 COMMON BLOCKS USED - C Real*8 PI, TWOPI, HALFPI, CONVD, CONVDS, CONVHS, SECDAY COMMON / CMATH / PI, TWOPI, HALFPI, CONVD, CONVDS, CONVHS, SECDAY EXTERNAL CMATHB C C VARIABLES 'FROM': C 1. CONVHS - THE CONVERSION FACTOR OF RADIANS PER TIME-SECOND. C (RAD/TIME-SEC) C INCLUDE 'ccon.i' C VARIABLES 'FROM': C 1. KDIUC - THE DIURNAL SPIN UTILITY ROUTINE FLOW CONTROL FLAG. C 2. KDIUD - THE DIURNAL SPIN UTILITY ROUTINE DEBUG OUTPUT FLAG. C 3. KNUTC - THE NUTATION MODULE FLOW CONTROL FLAG. USED TO C DETERMINE WHETHER OR NOT DPSI SHOULD BE ZEROED. C Real*8 GMSTC(4), SIDVEL(3), DPSIS(2), epsmd, omega, CENTJ, DJ2000, . GMSTMN, DAYSJ, CENT Integer*4 I C DATA CENTJ / 36525.D0 /, 1 DJ2000 / 2451545.D0 /, 2 GMSTC / + 24110.54841D0, + 8640184.812866D0, 3 + 0.093104D0, - 6.2D-6 /, 4 SIDVEL / + 1.002737909350795D0, + 5.9006D-11, - 5.9D-15 / C C 1.2.4 DATA BASE ACCESS - NONE C C 1.2.5 EXTERNAL INPUT/OUTPUT - C INPUT VARIABLES - NONE C OUTPUT VARIABLES - POSSIBLE DEBUG OUTPUT C C 1.2.6 SUBROUTINE INTERFACE - C CALLER SUBROUTINES: DRIVR C CALLED SUBROUTINES: DDROT, DROTT, ROTAT, EQEG C C 1.2.7 CONSTANTS USED - CENTJ, CONVHS, DJ2000, GMSTC(4), SIDVEL(3) C C 1.2.7.1 CONSTANTS INITIALIZED IN THIS UTILITY ROUTINE - C 1. CENTJ - THE NUMBER OF COORDINATE TIME DAYS PER JULIAN C CENTURY (CENTJ = 36525.D0 DAYS/CENTURY) C 2. DJ2000 - THE JULIAN DATE OF THE EPOCH JANUARY 1.5, 2000. C (DJ2000 = 2451545.D0 DAYS) C 4. GMSTC(4) - THE CONSTANTS APPEARING IN TERMS 1-4 IN THE C CALCULATION OF THE GREENWICH MEAN SIDEREAL TIME AT C ZERO HOURS UT1 OF THE DATE IN QUESTION. (TIME-SEC) C (SEE REFERENCES) C ( GMSTC(1) = + 24110.54841D0, C GMSTC(2) = + 8640184.812866D0, C GMSTC(3) = + 0.093104D0, C GMSTC(4) = - 6.2D-6 ) C 4. SIDVEL(3) - THE CONSTANTS APPEARING IN TERMS 1-3 IN THE C CALCULATION OF THE EARTH'S DIURNAL ANGULAR VELOCITY. C (TIME SEC/SEC) (SEE REFERENCES) C ( SIDVEL(1) = + 1.002737909350795D0, C SIDVEL(2) = + 5.9006D-11, C SIDVEL(3) = - 5.9D-15 ) C C 1.2.8 PROGRAM VARIABLES - C 1. CENT - THE NUMBER OF JULIAN CENTURIES ELAPSED SINCE THE C EPOCH JANUARY 1.5, 2000. (CENTURIES) C 2. DAYSJ - THE NUMBER OF JULIAN DAYS ELAPSED SINCE THE EPOCH C JANUARY 1.5, 2000. (DAYS) C 3. GMSTMN - THE GREENWICH MEAN SIDEREAL TIME AT ZERO HOURS UT1 C OF THE DATE IN QUESTION. (RAD) C 4. DPSIS(2) - THE VALUE OF NUTATION USED IN THE CALCULATION OF C SIDEREAL TIME. C 5. EPSMD - Mean obliquity of the ecliptic at date. (rad) C 6. OMEGA - Mean longitude of the ascending node of the Moon C (arcsec) C C 1.2.9 PROGRAMMER C 770207 DALE MARKHAM C 770718 PETER DENATALE C 771128 BRUCE SCHUPLER C 790201 BRUCE SCHUPLER C 810804 CHOPO MA C 870604 SAVITA GOEL C 890726 Jim Ryan: Documentation simplified. C 930901 Norbert Zacharias: Equation of equinoxes GAST update C C DIRNL Program Structure C C Compute mean obliquity of date, req. for EQE update in DIRNL epsmd = eps (1) - deps (1) C Mean longitude of ascending node of the Moon omega = fa (5) C Compute the number of Julian days elapsed since the epoch C JANUARY 1.5, 2000 to 0 hours UT1 of the date in question. DAYSJ = XJD - DJ2000 C C Compute the number of Julian centuries elapsed since Epoch January 1.5, 2000 CENT = DAYSJ / CENTJ C C Compute the diurnal angular velocity of the Earth. DIURNV = ( SIDVEL(1) + SIDVEL(2) * CENT + SIDVEL(3) * CENT**2 ) 1 * CONVHS C C Compute the Greenwich Mean Sidereal Time at 0 hours UT1 of the date C in question. (radians) GMSTMN = ( GMSTC(1) + GMSTC(2) * CENT + GMSTC(3) * CENT**2 + . GMSTC(4) * CENT**3 ) * CONVHS C C Compute the Greenwich Mean Sidereal Time. GMST = GMSTMN + DIURNV * UT1 GMST = DMOD(GMST,TWOPI) C C Check to see if the nutation in longitude is to be turned off. DO I=1,2 IF (KNUTC .EQ. 1) DPSIS(I) = 0.0D0 IF (KNUTC .NE. 1) DPSIS(I) = DPSI(I) ENDDO C C Compute the Greenwich Apparent Sidereal Time by calling the geometry part of C the equation of equinox procedure. CALL EQEG (gmst, dpsis(1), eps(1), epsmd, omega, . gast, gastd) C C Compute the time derivative of Greenwich Apparent Siderial Time GAST(2) = DIURNV * DUT1AT * DATDCT 1 + DCOS ( EPS(1) ) * DPSIS(2) 2 - DSIN ( EPS(1) ) * DPSIS(1) * EPS(2) C C Construct the diurnal spin matrix and its first two CT time derivatives. C C Construct the diurnal spin matrix. CALL ROTAT ( -GAST(1), 3, RS(1,1,1) ) C C Construct the first CT time derivative of the diurnal spin matrix CALL DROTT ( -GAST(1), -GAST(2), 3, RS(1,1,2) ) C C Construct the second CT time derivative of the diurnal spin matrix. CALL DDROT ( -GAST(1), GAST(2)**2, 3, RS(1,1,3) ) C C**************************************************************************** C Check KDIUD for debug output. IF ( KDIUD .NE. 0 ) THEN WRITE ( 6, 9) 9 FORMAT (1X, "Debug output for utillity DIRNL." ) WRITE(6,8)' CENT ',CENT 8 FORMAT(A,4D25.16/(7X,5D25.16)) WRITE(6,8)' CENTJ ',CENTJ WRITE(6,8)' DJ2000 ',DJ2000 WRITE(6,8)' DPSIS ',DPSIS WRITE(6,8)' CONVHS ',CONVHS WRITE(6,8)' DAYSJ ',DAYSJ WRITE(6,8)' GMSTC ',GMSTC WRITE(6,8)' GMSTMN ',GMSTMN WRITE(6,8)' SIDVEL ',SIDVEL C WRITE ( 6, 9200 ) DATDCT, DPSI, DUT1AT, EPS, UT1, XJD, 1 DIURNV, GAST, GMST, RS 9200 FORMAT (1X, "DATDCT = ", D30.16, /, 1X, 1 "DPSI = ", 2 ( D30.16, 10X ), /, 1X, 2 "DUT1AT = ", D30.16, /, 1X, 3 "EPS = ", 2 ( D30.16, 10X ), /, 1X, 4 "UT1 = ", D30.16, /, 1X, 5 "XJD = ", D30.16, /, 1X, 6 "DIURNV = ", D30.16, /, 1X, 7 "GAST = ", 2 ( D30.16, 10X ), /, 1X, 8 "GMST = ", D30.16, /, 1X, 9 "RS = ", 9 ( 3 ( D30.16, 10X ), /, 1X ) ) C ENDIF ! debug output flag C Return END C C****************************************************************************** SUBROUTINE EQEG (gmst,dpsi,eps,epsmd,omega, gast,gastd) IMPLICIT NONE C C This subroutine computes the geometry section of the equation of equinox C procedure. C C References: D.McCarthy, IERS Technical Note 13, Paris 1992 C C Calling sequence: C input variables: C gmst = Greenwich mean siderial time (radian) C dpsi = nutation in longitude (radian) C eps = true obliquity of the ecliptic (radian) C epsmd = obliquity of the mean ecliptic of date (radian) C omega = mean longitude of the asc. node of the moon (arcsec) C C output variables: C GAST(1)=Greenwich apparent siderial time (New or Old depending on KDIUC) C (radians) C GASTD = Difference in GAST of two versions (New-Old) (radians) C REAL*8 gmst, dpsi, eps, epsmd, omega, gast(2), gastd REAL*8 old_gast, aswas_gast, new_gast, . eqe_const(2), omega_rad REAL*8 pi, twopi, halfpi, convd, convds, convhs, secday COMMON /CMATH/ pi, twopi, halfpi, convd, convds, convhs, secday EXTERNAL CMATHB C C Variables from: convds = convert arcsec to radian C INCLUDE 'ccon.i' C Variables from: C KDIUC - Flow control flag for EQE routines C = 0, GAST defined by new equation of equinoxes definition C = 1, GAST defined by old equation of equinoxes definition C KDIUD - Debug flag for EQE routines C C Data base access: None C C Subroutine interface: C caller subr.: DIRNL C called subr.: None C C Program variables: C eqe_const = Constants in new version of GAST formula C depending on sine of omega (arcsec) C aswas_gast = GAST as calculated incorrectly before (CALC 7.6 and C earlier). (radian) C old_gast = GAST as calculated before (CALC 7.6 and earlier) but C corrected to use mean obliquity of date, to be used until C 1997.02.26. (radians) C omega_rad = Omega (longit. of asc. node of moon) (radians) C new_gast = GAST for use after 1997.02.26 as recommended by the IERS. C C Programmer: C 930901 Norbert Zacharias: first implementation, equat.of eq. GAST C update. C 94.01.04 David Gordon - Clean up documentation, change sign of GAST C correction. C 98.01.20 David Gordon - Updated to use new definition EQE C DATA eqe_const / 0.00264D0, 0.000063D0 / ! in arc sec C C compute old GAST as before, use true obliquity eps as before aswas_gast = gmst + dpsi * DCOS (eps) C C compute old GAST (formula in effect until 19970226) C but with correct epsilon (mean obliquity of date) old_gast = gmst + dpsi * DCOS (epsmd) C C compute new GAST as recomm. in IERS techn.note omega_rad = omega * convds new_gast = gmst + dpsi * DCOS (epsmd) . + eqe_const(1) * convds * DSIN (omega_rad) . + eqe_const(2) * convds * DSIN (2.0d0*omega_rad) C C Calc 9, use new definition If (KDIUC.eq.0) gast(1) = new_gast C For diehards, allow using old definition If (KDIUC.eq.1) gast(1) = old_gast C Define difference, Gastd, as (New - Old) gastd = new_gast - old_gast C C check debug output flag IF (KDIUD.NE.0) THEN WRITE (6,'(1x,a)') 'debug output for EQEG' WRITE (6,'(1x,a,d25.16)') 'gmst = ', gmst WRITE (6,'(1x,a,d25.16)') 'dpsi = ', dpsi WRITE (6,'(1x,a,d25.16)') 'eps = ', eps WRITE (6,'(1x,a,d25.16)') 'epsmd = ', epsmd WRITE (6,'(1x,a,d25.16)') 'omega = ', omega WRITE (6,'(1x,a,d25.16)') 'aswas_gast = ', aswas_gast WRITE (6,'(1x,a,d25.16)') 'old_gast = ', old_gast WRITE (6,'(1x,a,d25.16)') 'new_gast = ', new_gast WRITE (6,'(1x,a,d25.16)') 'gast(1) = ', gast(1) WRITE (6,'(1x,a,d25.16)') 'aswas-old = ', aswas_gast-old_gast WRITE (6,'(1x,a,d25.16)') 'gastd = ', gastd ENDIF C Return END C **************************************************************************** SUBROUTINE DIRNC (gastd, dut1p) IMPLICIT NONE C C This subr. computes the contributions to the delay and delay rate due to the C different versions of the equation of equinox procedure. The New version C (to be used after 26 Feb 1997) becomes the default in Calc 9.0. To get back C to the old definition, simply add the contribution Lcode OLDEQCON during C SOLVE analysis. C C Users can still use the old EQE definition by setting KDIUC=1. In that case, C a diferent contribution goes into the database, 'NEWEQCON'. When NEWEQCON C is added to the delays and rates in this case, they will be converted to C the new definition. Don't count on this being supported in SOLVE though. C C References: D.McCarthy, IERS Technical Note 13, Paris 1992 C E-Mail Memo, 'UT1 Errors', D. Gordon, 10 April, 1997 C C Calling sequence: C Input Variables: C GASTD = Differ. in GAST of the two versions (new-old) (radian) C DUT1P(2,2) = Partial derivatives of delay and rate wrt (A1-UT1) C (First index - delay and rate, second index - first and C second derivatives) (s/s, s/s**2) C [Note: The partial we want here is actually the delay and rate C derivatives wrt GAST. We approximate this by dividing C dut1p by -1.0027379 sidereal days/solar day.] C REAL*8 gastd, dut1p(2,2), deqec(2) REAL*8 pi,twopi, halfpi, convd, convds, convhs, secday COMMON /CMATH/ pi,twopi,halfpi, convd, convds, . convhs, secday EXtERNAL CMATHB C Variables From: C Convhs = conversion factor (radian/sec.of time) C INCLUDE 'ccon.i' C Variables From: C KDIUC - Flow control flag for EQE routines C = 0, GAST defined by new equation of equinoxes definition C = 1, GAST defined by old equation of equinoxes definition C KDIUD - debug flag for EQE routines C C Data base access: C PUT4 GASTD = differ. in GAST (new - old) (radian) C PUT4 DEQEC(2) = contrib. of GASTD to delay and delay rate C (sec, sec/sec) C Subroutine interface: C Caller Subr.: DRIVR C Called Subr.: None C C Programmer: C 930901 Norbert Zacharias: first implementation C 94.01.04 David Gordon - Fixed up documentation, sign of contributions C flipped in EQEG. C 94.06.02 David Gordon - Changed PUTR to PUT4. C 95.02.09 David Gordon - Corrected PUT4's of 'EQE DIFF' and 'EQE CONT'. C Only 8 charcters allocated, changed to 14. Bug C found by Bill Petrachenko using a Sun workstation C at the Dominion Radio Astrophysical Obs. (DRAO). C 98.01.20 David Gordon - Corrected old (Calc 8.x) delay and rate C contribution (Deqec) by dividing by -1.0027379 C siderial days/solar day. Changed database Lcode to C 'NEWEQCON' (from 'EQE CONT') and made it the C non-default. Added 'OLDEQCON' Lcode contribution C and made it the default. C C Calculate contributions to delay and rate C C Default case, contribution to convert from New to Old EQE definition If (KDIUC .eq. 0) then deqec (1) = -gastd * dut1p(1,1) / (-1.0027379 * convhs) deqec (2) = -gastd * dut1p(2,1) / (-1.0027379 * convhs) CALL PUT4 ('OLDEQCON ', deqec, 2,1,1) Endif C C Non-default case, contribution to convert from Old to New EQE definition If (KDIUC .eq. 1) then C Correct Calc 8.2 error by dividing by -1.0027379 deqec (1) = gastd * dut1p(1,1) / (-1.0027379 * convhs) deqec (2) = gastd * dut1p(2,1) / (-1.0027379 * convhs) CALL PUT4 ('NEWEQCON ', deqec, 2,1,1) Endif C C Output EQE difference to data base CALL PUT4 ('EQE DIFF ', gastd, 1,1,1) C C check flag for debug output: IF (KDIUD.NE.0) THEN WRITE (6,'(1x,a)') 'debug output for EQEC' WRITE (6,'(1x,a,2d25.16)') . 'dut1p = ', dut1p(1,1), dut1p(2,1) WRITE (6,'(1x,a,2d25.16)') . 'gastd, convhs = ', gastd, convhs WRITE (6,'(1x,a,2d25.16)') 'deqec = ', deqec ENDIF C Return END