Handle changing ecliptic plane

[cosinekitty]

I have a better understanding now of the obliquity between the equatorial plane and the ecliptic plane. From page 15 of Astronomy on the Personal Computer by Montenbruck and Pfleger:

[The] obliquity of the ecliptic ε ... is not constant, but is decreasing by about 47" per century. This primarily arises from slow alterations in the Earth's orbit as a result of perturbations caused by the other planets.

Before now, I thought the ecliptic plane was essentially constant over human timescales, and that only the Earth's equatorial plane was changing. The equatorial plane is precessing and nutating, but I just now learned these have only a tiny effect on obliquity. The major long-term change in the obliquity angle ε is caused by changes in the Earth's orbital plane.

The implications for Astronomy Engine:

1. EclipticGeoMoon returns ecliptic coordinates relative to the moving ecliptic plane. That means it is very accurate for predicting eclipses.
2. Ecliptic returns coordinates for the ecliptic plane of the year 2000. Any place I use it for predicting lunar nodes, transits, etc, will lose accuracy for years long before or long after 2000. These need to be addressed.
3. I need to assess whether to keep a fixed year 2000 ecliptic, replace these coordinates with ecliptic-of-date coordinates, or add yet another coordinate system. Perhaps I need to replace what I call "ECL" with two mnemonics, "ECD" and "ECJ"?
4. I don't want to break backward compatibility too much, but maybe this isn't really breakage as much as fixing what was inaccurate.

This will take careful study to plan an approach. I will use this issue as a running journal of ideas, and to track progress.

[cosinekitty]

Here is the development plan:

1. Create a new branch called ecliptic for the following work items.
2. We already report results using three ecliptic orientations: mean equinox of 2000 (the function Ecliptic), mean equinox of date (EclipticGeoMoon), and the true equinox of date (SunPosition).
3. Because ecliptic coordinates often imply some physical significance (e.g. eclipses and lunar phases), it makes sense to enable reporting ecliptic coordinates using the true equinox of date. Further, the API should lead people to using true equinox of date coordinates as the easiest/default option.
4. I want to minimize breakage of existing code that depends on Astronomy Engine. Instead of renaming existing functions, it will be more tolerable to cause existing functions to return slightly different numbers, especially if the newer numbers are more accurate for physical interpretations.
5. Existing functions like Rotation_EQJ_ECL do not accept a date parameter. Either changing its name or adding a parameter would break existing code. Therefore, it is better to leave "ECL" with its current meaning: ecliptic coordinates using the mean equinox of 2000.
6. I will add a new orientation code "ECT" to indicate ecliptic for true equinox of date.
7. I might add another code "ECM" to indicate ecliptic for mean equinox of date, although I would prefer not to. This is the current orientation of EclipticGeoMoon. I still need to decide whether EclipticGeoMoon should include a nutation correction (to produce ECT) or should remain ECM.
8. Research: existing JPL Horizons test data file generate/moonphase/moon_ecm.txt appears to be using ECT orientation: "Observer-centered IAU76/80 ecliptic-of-date longitude and latitude of the target centers' apparent position, with light-time, gravitational deflection of light, and stellar aberrations." (The "IAU76/80" refers to a nutation model.) Verify that converting the Moon's coordinates for nutation improves the accuracy of my unit test MoonVector().
9. Measure the difference in efficiency between the existing EclipticGeoMoon and the same code with a nutation correction. If the difference is less than 10% (as I suspect), go ahead and add it to EclipticGeoMoon.
10. Another research project: can I use an even simpler nutation model than IAU2000b? Calculating 77 terms every single time can be expensive. Can I get better than arcsecond precision much faster?
11. Otherwise, audit any internal use of EclipticGeoMoon to see if the nutation correction is really needed. I believe there are places where I'm just trying to prune out eclipse candidates, and that level of accuracy would be unhelpful if made too much slower.
12. Add new functions that generate rotation matrices to convert back and forth between ECT and EQJ. Hopefully (assuming I don't maintain any ECM coordinates in the public interface) I will not need to mention ECM in the documentation or add functions for them also. But that is an option if needed.

The summary of the plan: add support for ECT = ecliptic coordinates in true equinox of date: rotation matrix functions and new version of Ecliptic function. Convert EclipticGeoMoon to return ECT, just like SunPosition already does. Consider optimizing nutation calculations.

[cosinekitty]

Correcting EclipticGeoMoon from mean ecliptic to true ecliptic (by applying nutation) turns out to greatly decrease efficiency for very little gain of accuracy. With this correction, EclipticGeoMoon takes 3 times as long. In fact, I'm wondering if it's time to consider replacing the iau2000b nutation model with something much faster, even if it sacrifices an arcsecond or two of accuracy. But for now, I'm leaning toward treating all ecliptic coordinates of date as equivalent, whether actually referenced to the mean obliquity or the true obliquity.

[cosinekitty]

I found there is an even faster nutation that uses only the largest 13 terms of the IAU2000 model instead of 77. This comes from an "alternative" FORTRAN listing from NOVAS 3.1:

      SUBROUTINE NOD (T,DPSI,DEPS)
*
*     SUBROUTINE NOD VERSION 2.
*     IN LOW-ACCURACY MODE, THIS SUBROUTINE EVALUATES A SHORT
*     NUTATION SERIES AND RETURNS APPROXIMATE VALUES FOR NUTATION IN
*     LONGITUDE AND NUTATION IN OBLIQUITY FOR A GIVEN TDB JULIAN DATE.
*     IN THIS MODE, ONLY THE LARGEST 13 TERMS OF THE IAU 2000A NUTATION
*     SERIES ARE EVALUATED.  IN HIGH-ACCURACY MODE, THE STANDARD IERS
*     SUBROUINE IS CALLED TO EVALUATE THE FULL IAU 2000A NUTATION
*     SERIES.
*
*          T    = TDB TIME IN JULIAN CENTURIES SINCE J2000.0 (IN)
*          DPSI = NUTATION IN LONGITUDE IN ARCSECONDS (OUT)
*          DEPS = NUTATION IN OBLIQUITY IN ARCSECONDS (OUT)
*
*     NOTE:  IN LOW-ACCURACY MODE, MAX ERROR IN DPSI < 0.05 ARCSEC,
*     MAX ERROR IN DEPS < 0.02 ARCSEC, AVERAGE ERROR ABOUT 1/4 OF MAX.
*
*
      DOUBLE PRECISION T,DPSI,DEPS,PI,SECCON,T0,T1,DP,DE,
     .     X,EL,ELP,F,D,OM,ARG,DSIN,DCOS
      DIMENSION X(9,13)
      SAVE

      PARAMETER ( PI     = 3.14159265358979324D0 )
      PARAMETER ( SECCON = 180.D0 * 3600.D0 / PI )

*     T0 = TDB JULIAN DATE OF EPOCH J2000.0 (TT)
      DATA T0 / 2451545.00000000D0 /

*     LARGEST 13 TERMS OF IAU 2000A NUTATION SERIES, WITH PRECISION
*     OF COEFFICIENTS TRUNCATED
      DATA X / 0., 0., 0., 0., 1., -17.2064,-0.01747, 9.2052, 0.00091,
     .         0., 0., 2.,-2., 2.,  -1.3171,-0.00017, 0.5730,-0.00030,
     .         0., 0., 2., 0., 2.,  -0.2276,-0.00002, 0.0978,-0.00005,
     .         0., 0., 0., 0., 2.,   0.2075, 0.00002,-0.0897, 0.00005,
     .         0., 1., 0., 0., 0.,   0.1476,-0.00036, 0.0074,-0.00002,
     .         0., 1., 2.,-2., 2.,  -0.0517, 0.00012, 0.0224,-0.00007,
     .         1., 0., 0., 0., 0.,   0.0711, 0.00001,-0.0007, 0.00000,
     .         0., 0., 2., 0., 1.,  -0.0387,-0.00004, 0.0201, 0.00000,
     .         1., 0., 2., 0., 2.,  -0.0301, 0.00000, 0.0129,-0.00001,
     .         0.,-1., 2.,-2., 2.,   0.0216,-0.00005,-0.0096, 0.00003,
     .         0., 0., 2.,-2., 1.,   0.0128, 0.00001,-0.0069,-0.00000,
     .        -1., 0., 2., 0., 2.,   0.0123, 0.00000,-0.0053, 0.00000,
     .        -1., 0., 0., 2., 0.,   0.0157, 0.00000,-0.0001, 0.00000 /
*     REMAINING TERMS ALL HAVE AMPLITUDES < 0.01 ARCSECOND

*     GET METHOD/ACCURACY MODE
      CALL GETMOD ( MODE )

      IF ( MOD ( MODE, 2 ) .EQ. 0 ) THEN

*         HIGH ACCURACY MODE -- CALL IERS SUBROUTINE

          T1 = T * 36525.D0
          CALL NU2000A ( T0, T1, DP, DE )
          DPSI = DP * SECCON
          DEPS = DE * SECCON

      ELSE

*         LOW ACCURACY MODE -- EVALUATE SHORT NUTATION SERIES ABOVE

*         COMPUTATION OF FUNDAMENTAL ARGUMENTS
          CALL FUNARG ( T,   EL, ELP, F, D, OM )

          DPSI = 0.D0
          DEPS = 0.D0

*         SUM NUTATION SERIES TERMS
          DO 10 I = 13, 1, -1
              ARG = X(1,I) * EL
     .            + X(2,I) * ELP
     .            + X(3,I) * F
     .            + X(4,I) * D
     .            + X(5,I) * OM
              DPSI = ( X(6,I) + X(7,I) * T ) * DSIN ( ARG ) + DPSI
              DEPS = ( X(8,I) + X(9,I) * T ) * DCOS ( ARG ) + DEPS
  10      CONTINUE

*         ADD IN OUT-OF-PHASE COMPONENT OF PRINCIPAL (18.6-YEAR) TERM
*         (TO AVOID SMALL BUT LONG-TERM BIAS IN RESULTS)
          DPSI = DPSI + 0.0033D0 * DCOS ( OM )
          DEPS = DEPS + 0.0015D0 * DSIN ( OM )

      END IF

      RETURN

      END

[cosinekitty]

As of 576eea2245894421aa25209bdf12253daeba22ea, I have a 5-term truncated version of the nutation formula. It is much faster than the original 77-term IAU2000B formula, but with minimal loss of accuracy. Next I will try again to adapt EclipticGeoMoon to return coordinates in the true ecliptic and equator of date, rather than using the mean equinox. Hopefully this time the slowdown of EclipticGeoMoon will be minimal.

[cosinekitty]

As of 12e68a1931024a6e1db36fb8c48cf7b34d9e5c13, the EclipticGeoMoon functions now correct for nutation of the Earth's axis. This brings their worst case error to less than 1.9 arcseconds in ecliptic latitude, 6.0 arcseconds in ecliptic longitude, over the years 1900 to 2100, when compared with JPL Horizons data. The automated unit tests for the C version of Astronomy Engine now verify this. The tests for the other 4 languages compare output with the C version's output to verify numbers are consistent to 15 decimal digits.

[cosinekitty]

Fixed in v2.1.13 at commit 1a4f842764ca71e1aff49a58597a4fd0ddd7e23c.