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.
Astronomy Engine examples in Python
Camera
Suppose you want to photograph the Moon, and you want to know what it will look like in the photo. Given a location on the Earth, and a date/time, this program calculates the orientation of the sunlit side of the Moon with respect to the top of your photo image. It assumes the camera faces directly toward the Moon's azimuth and tilts upward to its altitude angle above the horizon.
Culmination
Finds when the Sun, Moon, and planets reach their highest position in the sky on a given date, as seen by an observer at a specified location on the Earth. Culmination is also the moment a body crosses the meridian, the imaginary semicircle in the sky that passes from due north on the horizon, through the zenith (straight up), and then toward due south on the horizon.
Horizon Intersection
This is a more advanced example. It shows how to use coordinate transforms to find where the ecliptic intersects with an observer's horizon at a given date and time.
Jupiter's Moons
Calculates the coordinates of Jupiter and its four major moons (Io, Europa, Ganymede, and Callisto) as seen from the Earth at a given date and time. This program illustrates how to correct for the delay caused by the time it takes for light to reach the Earth from the Jupiter system.
Lunar Eclipse
Calculates details about the first 10 partial/total lunar eclipses after the given date and time.
Moon Phase Calculator
This example shows how to determine the Moon's current phase, and how to predict when the next few quarter phases will occur.
Positions
Calculates equatorial and horizontal coordinates of the Sun, Moon, and planets.
Rise/Set
Shows how to calculate sunrise, sunset, moonrise, and moonset times.
Seasons
Calculates the equinoxes and solstices for a given calendar year.
API Reference
Complete documentation for all the functions and types available in the Python version of Astronomy Engine.