Changed the documentation for the GeoMoon and GeoMoonState
functions to make it explicit that they calculate coordinates
oriented with respect to the Earth's J2000 equator (EQJ).
This is because I will soon add ecliptic (ECL) counterparts
for the GeoMoon function, to more directly search for ascending
and descending nodes of the Moon.
See this discussion:
https://github.com/cosinekitty/astronomy/issues/150
For the case of calculating a map, where each pixel
on the map represents a different location on the Earth,
it is more efficient to factor out expensive calculation
of sidereal times, assuming the entire map represents
some phenomenon at a single moment in time.
For example, to determine whether the Moon is visible
at different places on the Earth, the following
functions can be calculated across thousands of
different (lat, lon) geographic coordinates around
the world:
ObserverVector
Rotation_EQD_HOR
Before iterating over the map pixels, a program
can call GeoMoon, then convert EQJ coordinates to EQD.
Then by passing the same time value in a loop to
ObserverVector and Rotation_EQD_HOR, the program
can calculate a vector from the observer to the Moon
in EQD coordinates, then convert EQD to HOR.
The z-coordinate of the horizontal coordinates
determines whether the Moon is above or below the
observer's horizon at that point on the Earth.
This calculation pattern performed redundant
sidereal time calculations for each pixel on the map.
I changed the code for all 4 languages to cache
sidereal time so that it only needs to be calculated
once.
In the C version of Astronomy Engine, this resulted
in a speedup factor of about 2.3 in the above use case.
(See the function MapPerformanceTest in generate/ctest.c.)
The code generator was creating slightly different numeric
values for the Pluto state tables and the Jupiter rotation matrix.
I decreased the output precision by one decimal digit.
This should allow the code generator to produce identical
source code on both Linux and macOS.
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.
The documentation for the Python function `SearchAltitude`
was missing a mention of the `altitude` parameter.
I searched for similar mistakes in Python, C#, and C,
having just completed the same exercise in the JavaScript code.
I also found several places where extraneous newlines
between the parameter documentation caused the Markdown
to be rendered incorrectly.
Now the Python version of Astronomy Engine supports calculating
the Earth/Moon Barycenter (EMB) state vector (position and velocity)
relative to the Earth's center (geocentric) or relative
to the Solar System Barycenter (SSB).
This completes support for this feature across C, C#, JavaScript, and Python.
The BaryState function did not support Pluto before.
Refactored the code so that the internal CalcPluto function
returns both the position and velocity, and its caller
can select from heliocentric or barycentric coordinates.
HelioVector asks for heliocentric coordinates and keeps
only the position vector. BaryState asks for barycentric
coordinates and returns both position and velocity.
I added test data for Pluto generated by JPL Horizons.
It turns out the Pluto system barycenter is the best fit
for TOP2013, presumably because Charon causes Pluto to
wobble quite a bit.
I also generated JPL Horizons test data for the Moon
and the Earth/Moon barycenter, anticipating that I will
support calculating their barycentric state vectors soon.
I had to increase the enforced size limit for minified
JavaScript from 100000 bytes to 120000 bytes.
I guess this is like raising the "debt ceiling".
Fixed a bug in Python unit tests: if "-v" verbose option
was specified, it was printing a summary line for every
single line of input, instead of a single summary after
processing the whole file, as was intended. This is one
of those Python whitespace indentation bugs!
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.
The Pluto gravity simulator constants now come from
a single source: pluto_gravsim.h. This will allow me
to experiment with the Pluto state table to get a better
compromise between size and accuracy.
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.
There is now a Libration function in all 4 supported languages.
The returned structure contains libration angles in
ecliptic latitude and ecliptic longitude, along with
the Moon's ecliptic position and distance.
Also included is the Moon's apparent angular diameter.
All 4 languages have added a `diam_deg` field to the
structure returned by the Libration function.
It is the apparent angular diameter of the Moon as
seen from the center of the Earth, expressed in degrees.
In JavaScript and Python, throw an exception if provided
an invalid refraction option. Especially in JavaScript,
it was too easy to pass in a value like 'true', which did
not calculate refraction as expected.
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.
Added more exhaustive testing of VectorObserver.
I found a few cases where the height calculation
was off by more than 5 millimeters.
In the VectorObserver function, require the latitude solver
to keep iterating until the error is less than one billionth
of a degree. Now the height error is always within 1 mm.
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.
Ported conversion to/from galactic coordinates to Python.
Added unit test for new Python code.
Updated documentation for all 4 supported languages.
Fixed mistakes in JavaScript function documentation.
I wrote a quick Python program based on an original reference
paper defining the galactic orientation system.
It generates a rotation matrix from first principles
that matches one inside the NOVAS function equ2gal(),
within the expected 2.3 arcsecond difference between
ICRS and EQJ.
NOVAS equ2gal matrix:
double ag[3][3] = {
{-0.0548755604, +0.4941094279, -0.8676661490},
{-0.8734370902, -0.4448296300, -0.1980763734},
{-0.4838350155, +0.7469822445, +0.4559837762}};
This program's generated matrix:
B1950 = 1949-12-31T22:09:21.346Z
-0.0548624779711344 0.4941095946388765 -0.8676668813529025
-0.8734572784246782 -0.4447938112296831 -0.1980677870294097
-0.4838000529948520 0.7470034631630423 0.4559861124470794
Also added some JPL Horizons test data to confirm
conversion back and forth between EQJ and GAL, which
I will use for future tests.
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.
Started work on a Python demo for finding when the moon
reaches relative longitudes with other solar system bodies
that are multiples of 30 degrees. It is not finished yet,
but getting close.
Added operator overloads for the Python Time class so
that times can be compared against each other.
This makes it easier to sort a list of times, for example.
This script will help me figure out how to tune the search algorithm
I'm about to write for searching for times that a pair of bodies
reach a given relative apparent ecliptic longitude.
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 change has no effect on client-facing behavior.
It just makes the internal data tables for the array of
constellation appear more compact in C, C#, and Python.
This is what the TypeScript/JavaScript code was already doing.
