I'm getting much better accuracy sticking with my original
gravity simulator, just with smaller time increments, than
I was with the Runge-Kutta 4 method. The PlutoStateTable
gets a bit larger (51 state vectors instead of 41), but the
accuracy is so much higher.
Removed the Runge-Kutta code because I won't be going back to it.
I realize some use cases require adjustments for
stellar aberration. The existing aberration adjustments
are only supplied for calculating planet positions.
Some users will benefit from being able to add/subtract
aberration corrections to arbitrary vectors, including
for star positions.
I have added some JPL Horizons test data to help
validate the aberration functionality I'm about to add.
I created the beginning of a unit test in ctest.c,
but currently there is no aberration correction
implemented, so the test has no error threshold.
Allow this demo to use the current date and time by
default if the user does not specify one on the command
line. This required changing the order of the command line
parameters.
Given the right ascension and declination of a star,
expressed in J2000 coordinates, this demo converts those coordinates
to right ascension and declination expressed in the Earth's
equator at any given date and time. This example illustrates
how to use rotation matrices to convert one coordinate system
to another.
This example was prompted by the question at:
https://github.com/cosinekitty/astronomy/discussions/114
Before making these changes, I had the following discrepancies
between the calculations made by the different programming
language implementations of Astronomy Engine:
C vs C#: 5.55112e-17, worst line number = 6
C vs JS: 2.78533e-12, worst line number = 196936
C vs PY: 1.52767e-12, worst line number = 159834
Now the results are:
Diffing calculations: C vs C#
ctest(Diff): Maximum numeric difference = 5.55112e-17, worst line number = 5
Diffing calculations: C vs JS
ctest(Diff): Maximum numeric difference = 1.02318e-12, worst line number = 133677
Diffing calculations: C vs PY
ctest(Diff): Maximum numeric difference = 5.68434e-14, worst line number = 49066
Diffing calculations: JS vs PY
ctest(Diff): Maximum numeric difference = 1.02318e-12, worst line number = 133677
Here is how I did this:
1. Use new constants HOUR2RAD, RAD2HOUR that directly convert between radians and sidereal hours.
This reduces tiny roundoff errors in the conversions.
2. In VSOP longitude calculations, keep clamping the angular sum to
the range [-2pi, +2pi], to prevent it from accumulating thousands
of radians. This reduces the accumulated error in the final result
before it is fed into trig functions.
The remaining discrepancies are largely because of an "azimuth amplification" effect:
When converting equatorial coordinates to horizontal coordinates, an object near
the zenith (or nadir) has an azimuth that is highly sensitive to the input
equatorial coordinates. A tiny change in right ascension (RA) can cause a much
larger change in azimuth.
I tracked down the RA discrepancy, and it is due to a different behavior
of the atan2 function in C and JavaScript. There are cases where the least
significant decimal digit is off by 1, as if due to a difference of opinion
about rounding policy.
My best thought is to go back and have a more nuanced diffcalc that
applies less strict tests for azimuth values than the other calculated values.
It seems like every other computed quantity is less sensitive, because solar
system bodies tend to stay away from "poles" of other angular coordinate
systems: their ecliptic latitudes and equatorial declinations are usually
reasonably close to zero. Therefore, right ascensions and ecliptic longitudes
are usually insensitive to changes in the cartesian coordinates they
are calculated from.
This caused me to discover I had forgotten to finish
making the necessary changes to astronomy.ts for saving
the cartesian vector inside the EquatorialCoordinates class.
I also realized I had made a mistake in the documentation
for the y-coordinate of the vector: it is the June solstice;
there is no such thing as a September solstice!
Also fixed some mistakes in demo tests: if something failed,
I was printing out the wrong filename (camera.c instead of camera.cs).
In all four versions of Astronomy Engine (C, C#, JavaScript, and Python),
starting a search for a full moon near December 19, 2020 would fail.
I added a unit test to all four languages and it failed consistently
across them all.
The root cause: I was too optimistic about how narrow I could make
the window around the approximate moon phase time in the
SearchMoonPhase functions. Finding the exact moon phase time failed
because it was outside this excessively small window around the approximate
time. I increased the window from 1.8 days to 3.0 days.
This should handle all cases with minimal impact on performance.
Now all four of the new unit tests pass.
I believe this wraps up the Python integrator.
It now works in all 4 languages and passes all tests.
Fixed up demo tests to match new output.
Turned on Travis CI checking in this branch again.
To be consistent, when calculating the geocentric position of the Sun,
we do need to correct for light travel time just like we would for any
other object. This reduces the maximum time error for predicting transits
from 25 minutes to 11 minutes.
Also had to disable aberration when calculating moon phases
(longitude from Sun) in order to keep a good fit with test data.
In all 4 supported languages, use consistent constant names for
Earth and Moon radii.
Use Moon's equatorial radius for rise/set timing.
Use Moon's mean radius for calculating Moon's umbra radius for
detecting solar eclipses.
Also use Moon's mean radius for determining whether the Earth's shadow
touches the Moon, for finding lunar eclipses.
Use the Moon's polar radius for distinguishing between total
and annular eclipses, with a 14 meter bias (instead of 1420 meters!)
to match Espenak data.
Use consistent unit test error threshold of 0.57 minutes for rise/set.
Updated demo test data for slight changes to rise/set prediction times.
Updated doxygen options to issue an error on any warnings.
Fixed the incorrect function name link that doxygen was warning me about.
