Added the following new functions to all 4 languages:
MassProduct: find the GM product for all Solar System bodies.
LagrangePoint: calculate L1..L5 state vectors for a pair of bodies.
LagrangePointFast: calculate L1..L5 state vectors given
state vectors and GM products of a pair of bodies.
Reworked the tests that use JPL Horizons output files containing
state vectors so that they generalize to different parameters.
Specifically, soon I will need to pass in (major_body, minor_body,
point) to support Lagrange point tests.
I want to be able to re-use the code for loading state
vectors from a JPL Horizons text file, instead of copy-n-paste
like I did in C. So I reworked it as an iterator.
There were a few places inside the unit test function LoadStateVectors
where I had error messages that printed the wrong function name
(VerifyStateBody instead of LoadStateVectors) if an error was detected.
This is because of a copy-n-paste oversight. They are fixed.
In languages that support it, using hypot(x,y) is a little
easier to read than sqrt(x*x + y*y). Some documentation
(e.g. the man page for the C function) leads me to believe
hypot might also be better behaved than sqrt in some cases.
The JavaScript Math.hypot() is especially nice because it works
for any number of dimensions, so I can use it in 2D and 3D cases.
C only allows 2D usage, as does Python 3.7. Python 3.8 added
support for any number of dimensions, but I don't want to break
compatibility with Python 3.7 just yet. Therefore, in C and Python,
I am only using hypot for 2D cases.
C# does not appear to have any kind of hypot function,
so no changes were made to the C# code.
Thanks to https://github.com/ebraminio for this suggestion.
There is no function double.IsFinite() in .NET Framework.
Reworked the sanity check in Astronomy.Pivot so the C# code
builds in these older .NET platforms.
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.
My custom Markdown documentation generator for C had
a bug when emitting the listing of a #define.
It is not valid to try to hyperlink to other symbols,
because the Markdown syntax gets listed literally inside
the context of a C code block.
In most cases, people calculating Lagrange points will just
want to pass in the bodies and not have to worry about calculating
their state vectors and masses.
Renamed Astronomy_LagrangePoint to Astronomy_LagrangePointFast.
Added new function Astronomy_LagrangePoint that accepts body enum
values instead of state vectors and masses. It knows to optimize
the precision of the calculation by calling GeoMoonState for the
Earth/Moon case.
I confirmed that the Mac version of GitHub Actions does not
flush stderr on exit, so that is why I couldn't see my
diagnostic error messages. Now I can tell VerifyEquilateral
is failing due to a slightly out of bounds length ratio
in the Earth/Moon/M4 equilateral triangle. Made the tolerance
window a little larger, and trying again.
I think I have finally tracked down where the 5.4 arcminute
discrepancy comparing my L4 with JPL Horizons L4 is coming from.
When I feed JPL's geocentric Moon state vector through my
L4/L5 calculator, the result is in a slightly different plane
than it should be. It looks like a mistake in JPL Horizons!
I also fixed a bug where ctest's LoadStateVectors() was not
initializing the state.status before appending to the array.
The result was uninitialized random garbage in the status.
It is conceptually simpler to take cross products to
generate 3 coordinate axes (essentially a rotation matrix)
that represent radial, tangential, and normal directions
with respect to the major and minor bodies.
Before comparing my Lagrange point calculations to JPL Horizons,
I check to see if my Lagrange point calculations form an
equilateral triangle from (major body, minor body, L4/L5).
When ctest.c detects that a state vector error is
unacceptably large, it now prints extra diagnostics
about the two vector values, their magnitudes, and
how much of the error is angular and how much is
a magnitude discrepancy.
Calculated the angles between JPL Horizons velocity vectors
for the geocentric Moon and the geocentric L4/L5.
They are always very close to 60 degrees apart, within 0.15 arcsec.
This is not large enough to explain the velocity vector
errors my code calculates.
The Lagrange test was using Solar System barycentric
state vectors for the pair of bodies. This involved
a lot of unnecessary calculation.
For the Sun/EMB test, use heliocentric coordinates.
For the Earth/Moon test, use geocentric coordinates.
Fail the lagrange_jpl test if we see the major/minor/L4,L5
triangle deviate more than a tiny fraction from equilateral,
after adding the velocity components.
This convinces me that the JPL Horizons velocity vectors
make sense for L4/L5.
I'm having problems confirming formulas for L4/L5 velocity vectors.
So I wanted to test the assumption that these Lagrange points would
have velocity vectors that would leave the major body, minor body,
and Lagrange point in an equilateral triangle after all 3 bodies
continued in a straight line at their current relative velocities.
This does turn out to be the case, which means there is just
a bug in how I'm calculating the velocity vectors.
I just need to find the bug.
Now correctly calculating L4 and L5 positions, but
there is a large error in their velocity vectors.
Refactored ctest.c LagrangeTest() to be a lot easier
to understand and modify. A new function VerifyStateLagrange()
allows passing test parameters in a more function-oriented way.
Confirmed that L4 and L5 always lie in the same plane with
the position vector and velocity vector.
By taking the cross product of the Moon's position with
its velocity, I confirmed this is the same normal vector
as taking the Moon's position crossed with the L4/L5 position.
This means I should be able to calculate L4 and L5
using position and velocity vectors to define the
instantaenous co-orbital plane.
Use the formulas I already had to calculate first
approximations for L1, L2, L3 distances.
Then use Newton's Method to home in on the positions
where centrifugal acceleration balances with net
gravitational acceleration.
I realized there was a small mistake in how I was
calculating the distance scaling factor for L1 and L2.
It was relative to the distance between the minor body
and the barycenter, not the minor body and the major body.
This significantly improves the accuracy for Earth/Moon
Lagrange points, but still has more error compared
to JPL Horizons than I currently understand.
Analysis of Lagrange point calculations by JPL Horizons
now includes standard deviations of position/velocity
magnitude ratios. They confirm the ratios are very consistent.
The Microsoft C compiler is oddly picky about declaring a const variable.
Apparently it cannot do math with other const variables in its
initializer expression, unlike other C compilers.
So I had to change MOON_GM from a const to a #define.
Wrote code in ctest.c to load state vector data from
JPL Horizons output into a dynamically-allocated array.
This makes it easier to detangle the logic for loading
the data from the logic for doing statistical analysis.
The Lagrange point calculation is still not finished,
but L1 and L2 are working. L3 is probably correct, but there
is no test data for it.
I replaced the test data with new JPL Horizons output that
is centered on the primary body instead of the Solar System Barycenter.
This allows Astronomy_LagrangePoint() to be agnostic about
the coordinate systems of the state vectors handed to it.
I still need to get L4 and L5 calculations to match JPL Horizons
data, but it is not yet clear how to do that.
Python and npm package version: 2.0.11.
Finished implementing new functions across all
supported languages:
EclipticGeoMoon
Calculate the Moon's ecliptic geocentric position
in angular coordinates. The ecliptic longitude is
measured with respect to the mean equinox of date.
SearchMoonNode
NextMoonNode
A pair of functions to search for consecutive occurrences
of the Moon's center passing through the ecliptic plane.
Implemented a pair of C functions for finding a series of
Moon nodes:
Astronomy_SearchMoonNode
Astronomy_NextMoonNode
Finished the C unit test "moon_nodes" that verifies
my calculations against Fred Espenak's test data.
Fixed another bug in my parsing of the original Espenak
data for moon nodes. Added verification that my own
Moon calculations match his:
The Moon is always within 0.182 arcminutes ecliptic
longitude of the node when he says it is crossing the node.
The Moon is always within 1.54 arcminutes of the equatorial
coordinates he says.
