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.
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.
This is a thin wrapper function for the internal
function CalcMoon, which has already been extensively
validated. It will enable outside users to search
for ascending and descending nodes of the Moon,
or to calculate ecliptic spherical coordinates for the Moon
for any other useful purpose.
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.
Added data for testing Lagrange point calculations,
which has yet to be implemented. JPL Horizons
provides L1, L2, L4, and L5 (but *not* L3) calculations
for:
Sun / EMB
Earth / Moon
where EMB = Earth/Moon Barycenter.
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.)
Reduce the number of redundant Earth nutation calculations
by passing astro_time_t values as pointers in more functions.
Nutation values can then be cached in the time parameter
and passed to other functions that can then avoid calculating
the same nutation again.
Nutation is an expensive calculation, so reducing this overhead
can dramatically speed up certain use cases.
This was only needed in C, because this is the only language
in which times are passed by value. In Python, C#, and JavaScript,
times are objects that are already passed by reference, and
they already benefit from this nutation recyling approach.
The following functions have had their parameters changed.
This is a breaking change, but in every case, the caller
usually just needs to change `time` to `&time`.
Astronomy_Rotation_EQD_EQJ
Astronomy_Rotation_EQD_ECL
Astronomy_Rotation_EQD_HOR
Astronomy_Rotation_EQJ_EQD
Astronomy_Rotation_EQJ_HOR
Astronomy_Rotation_ECL_EQD
Astronomy_Rotation_ECL_HOR
Astronomy_Rotation_HOR_EQD
Astronomy_Rotation_HOR_EQJ
Astronomy_Rotation_HOR_ECL
Astronomy_RotationAxis
Astronomy_VectorObserver
For developers wanting to contribute to Astronomy Engine,
the tooling instructions were vague about xsltproc.
I added some useful hints.
I also removed things that don't need to be manually
installed: jsdoc2md, graphviz, etc.
In Windows builds, I was checking for the existence
of md5sum.exe, and if present, I used it for verifying
downloaded files. However, this is not a standard
utility that comes built into Windows 10.
I found there is a standard utility certutil.exe
that can calculate md5, sha256, etc, checksums.
However, it does not verify files created by the
Linux utilities md5sum, sha256sum, etc.
So I created a batch file checksum.bat that invokes
certutil.exe to process one of those listing files.
Reworked run.bat to call checksum.bat instead of
using md5sum.exe.
Weirdly, the python program for generating constellation
data did not seem to run, but that failure did not break
the build directly. I am adding an explicity 'py' command
to run each Python program. Also added a check for missing
output constellation test data.
The version of md5sum.exe that runs in GitHub Actions
gets confused when its input file ends with CRLF,
which is weird because this is supposed to be a Windows
environment. So I added a .gitattributes to force
these checksum files to have LF line endings.
In the GitHub Actions CI environment, the program md5sum.exe
resides in "C:\Program Files\Git\usr\bin\md5sum.exe".
Therefore, to execute it, I need quotes around the executable.
The 'windows' directory is mainly useful for
maintainers, not end users. So I moved it out of
the root to reduce distraction for a first-time
visitor.
While I was fixing up resulting breakage in
Visual Studio project files, I noticed I still had
some hard-coded absolute paths that would only work
on my own Windows laptop (e.g. "c:\don\github\astronomy").
I replaced those with relative paths that will work
regardless of what directory the repo is cloned into.
