The demo program ecliptic_of_date.py shows how to
calculate the true ecliptic of date (ECT) angular coordinates
of the Sun, Moon, and planets for an observer somewhere on the Earth.
It calculates the equatorial of date (EQD) coordinates, then uses
a rotation matrix to convert the vector to ECT, then converts
the vector to spherical coordinates: latitude, longitude, and distance.
I refactored the unit tests for all the demo programs
to follow a different pattern that makes it simpler
to add more demo tests in the future.
The main thing is that correct output and generated
output are now in separate directories `correct` and `test`.
I have moved the test scripts from `test/test` to `./demotest`
in all the langauge demo directories.
This makes it simpler to clean up any stale generated
files before each test run by `rm -f test/*.txt`.
I stumbled across this while making the Java demo tests,
and it was a better solution, so now all the other languages
are consistent with the Java demo tests.
In the C demo tests, I also decided to compile all the
binary executables into a subdirectory `bin` that can
be cleaned out before each run, to make sure there are
no stale executables from an earlier run.
For years before 1582 or years after 3668, the Seasons functions
were unable to find many equinoxes and/or solstices.
The problem was that over time, the Earth's axis precesses
enough that the calendar dates of these events drifts outside
the fixed search ranges I had provided for them.
I expanded the search ranges so all season changes can be found
for a much wider range of years, as verified by unit tests:
C/C++: -2000..9999
C#: 1..9999
JavaScript: -2000..9999
Python: 1..9999
Kotlin: 1..9999
Note: C#, Python, and Kotlin currently do not allow
years values below +1. In fact, I discovered we were not
noticing when an invalid year was passed into the Kotlin code.
I updated that code to throw an exception when the year does
not match what was expected. It is disturbing that the
GregorianCalendar class silently ignores invalid years!
Constricted the search tolerance from 1 second to 0.01
seconds for the seasons search, to ensure more consistent
behavior.
Fixed a bug in the Kotlin search() function's
quadratic interpolation that was causing the convergence
to be slower than it should have been.
The phrase "Moon phase" is ambiguous, because sometimes
it means relative ecliptic longitude, other times it means
illuminated fraction. The "moonphase" demos were only
calculating the relative ecliptic longitude, which was
confusing. Now they calculate both.
More work getting MacOS build process to work.
Avoid excessive number of floating point digits of
output in the demo tests, so that insignificant
floating point variations don't cause unit test failures.
I found a mistake in the raytracer's Spheroid class,
thanks to a warning about an unused member variable.
I don't believe it had any effect on the currently
generated images, but it was important to fix it before
I ever do any set operations on Spheroids.
On macOS, there is no 'realpath' command by default.
So I eliminated some more attempts to use 'realpath'
in the demo test scripts.
Renamed the GitHub Actions tests to be consistent:
Astronomy-Engine-Linux
Astronomy-Engine-Macos
The demo tests on Mac OS failed because of very tiny
floating point discrepancies that don't matter.
Changed the output of the "Moon check" so that slight
differences in vector residue no longer fail the unit tests.
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.
Tonight as I was walking outside, I saw a fairly bright
star about half a degree away from the edge of the Moon.
I wondered what it was, so I decided to write a quick
program to find out.
This Python demo program scans the HYG Database
(https://github.com/astronexus/HYG-Database)
to find which bright stars are within a small angular
distance of the Moon, as seen at a given time, latitude, and longitude.
It turns out the star I saw was Nunki (Sigma Sagittarii).
It was handy to do vector subtraction to implement this program,
and it was trivial to do in the Python code's Vector class,
so I went ahead and added that.
Just like the Python version, this program calculates
the best-fit intersection point for two lines of sight
as seen by two observers. It demonstrates converting
back and forth between geographic coordinates and
geocentric vectors.
Instead of the hack call to Search(), the latitude
solver now uses Newton's Method directly. This
significantly speeds up the code, and is more elegant.
The program triangulate.py finds the point in space
where two vantage lines come closest to each other.
It is the midpoint between the closest points on both lines.
Now I print the distance between those two points also,
as a measure of how much uncertainty there is in the
estimation of the target object.
I already had the function ObserverVector that converts geographic
coordinates (latitude, longitude, elevation) to an equatorial-of-date
(EQD) vector.
Now I'm in the process of adding the inverse function VectorObserver
that calculates geographic coordinates from an EQD vector.
This commit implements VectorObserver in Python.
The other languages will follow in future commits.
The motivation was from the following request:
https://github.com/cosinekitty/geocalc/issues/1
The goal is to find the near-intersection between two different lines
of sight from two different observers on the Earth's surface.
Added a demo program triangulate.py that solves this problem.
Finished the script demos/python/lunar_angles.py
that shows how to search for times when the Moon and other
solar system bodies reach apparent ecliptic longitude separations
as seen from the Earth.
This is also a good demo of how to perform a custom search
for events using Astronomy Engine. This is the same technique
used internally by Astronomy Engine to search for lunar phases,
eclipses, solstices, etc.
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.
The demo shows how to correct for light travel
time to render Jupiter's moons as they appear
from the Earth.
Created an addition operator for the Vector
class in the Python code, because it is handy.
Corrected a bug in the string representation
of the Python StateVector class.
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.
Created skeleton test harness for validating the demo programs.
Created stub moonphase.py.
Copied correct demo program outputs from nodejs; will tweak as needed.
Call the Python demo test harness from the 'run' script.
