Made the internal data structure for the C version of barystate
use named fields in a struct for Sun...Neptune, instead of
an array. This makes the C code look more like the other 3
language implementations. I am going to experiment with adding
more bodies to see if it helps accuracy, and this makes the
code easier to modify for that experiment.
I am starting to work on a function to find the position
and velocity vectors for an observer on the surface of the Earth.
I created the C function Astronomy_ObserverState(), but I don't
yet have a unit test for it.
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.
Now the C 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).
Astronomy_BaryState() now supports body == BODY_EMB.
Added a new function Astronomy_GeoEmbState() to calculate
the geocentric state for the EMB.
Both have been verified using test data generated by JPL Horizons.
I wanted to test the geocentric Moon state vector I
calculate for the sake of the Moon state relative to
the Solar System Barycenter (SSB). Because the geocentric
portion has such a small magnitude, I decided to go ahead
and expose GeoMoonState as part of the API, and create
a test for it specifically. I used JPL Horizons to generate
the test file GeoMoon.txt.
The C function Astronomy_BaryState() now supports BODY_MOON.
Because of the complexity of the CalcMoon() function, I ended
up calculating two positions close together in time, and
using dr/dt to estimate the velocity vector.
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.
Reworked the Pluto gravity sim constants so they are defined
in one place: a new header file gravsim/pluto_gravsim.h.
Then the code generator writes the #defines to the C code, instead
of having two independent versions of the same constants.
I will continue down the road of having a single-source-of-truth
for these constants across all 4 supported languages.
Also, confusingly, I had one constant called PLUTO_DT in codegen.c
that was called PLUTO_TIME_STEP in astronomy.c. Also, astronomy.c
had a different constant PLUTO_DT that didn't mean the same thing.
I reworked the naming to be consistent in all places.
I already had a TopPosition() function that knows how to calculate
exact equatorial coordinates, so I eliminated the redundant logic
from gravsim_test.c
Significantly decreased the calculation error:
0.20 arcmin to 0.12 arcmin in my test metric.
However, the amount of extra work may not be
worth the accuracy, compared to just stepping more
increments between the segments, or simply making
more segments in the first place.
As they say in government-funded academia,
"more research is needed."
I'm interested in understanding drift error calculation
between the known-correct state waypoints for Pluto's orbit.
I'm trying to figure out whether there is some unintended
asymmetry between the forward calculation and the reverse
calculation. I will likely have to compare against TOP2013
data for the major planets, because I am using truncated
VSOP87, which also introduces position errors.
I should also validate TOP2013 versus VSOP87, both
untruncated, for Jupiter..Neptune, which they both cover.
I'm trying to get a better feel for the amount of error
in my gravity simulator calculations for the movement of Pluto.
Added conditionally-compiled code to log state vectors calculated
in the forward and reverse time directions, along with the
exact endpoints that frame the interpolated values.
Also log errors measured between both directions.
There is a curious asymmetry in the first case I tried
(roughly the years 2000..2100), where the forward calculation
seems less accurate than the reverse calculation.
For bodies other than the Sun, Moon, or Saturn, the
C function Astronomy_Illumination calls an internal function
VisualMagnitude. If VisualMagnitude is passed an invalid body,
it returns an error code. Astronomy_Illumination was not checking
for an error code. In the case of being passed a pseudo-body
like BODY_SSB (the Solar System Barycenter) or BODY_EMB
(the Earth/Moon Barycenter), VisualMagnitude is called and
returns the error ASTRO_INVALID_BODY. Astronomy_Illumination was
ignoring the error and returning ASTRO_SUCCESS to the caller,
even though the magnitude was NAN (the "not a number" value).
Note that other invalid bodies than EMB and SSB would not
cause this problem, because the earlier call to HelioVector
would fail and be noticed.
I added unit tests that confirmed this bug, then made fixes
to the code so that the unit tests pass.
I confirmed this same problem does NOT exist in the Python,
JavaScript, or C# versions of Astronomy Engine. In all
the other languages, this case causes VisualMagnitude to
throw an exception, so no error checking is needed in the
Illumination functions.
The C and C# Illumination functions now return
a `phase_fraction` result to complement `phase_angle`.
This makes them consistent with the Python and JavaScript
versions.
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.
Because I have to perform the expensive calculation to find
the Moon's ecliptic coordinates, I might as well return them
to the caller. This could help reduce calculation overhead
for some uses, and doesn't add any significant cost.
I could not measure a significant difference in calculation
accuracy from doing the expensive earth-tilt step.
I removed it to significantly speed up the calculation.
I discovered that when I tried to build astronomy.c as C++ code,
I got several errors and warnings. So I fixed those issues and
added the C++ build-check to the unit tests.
I'm doing Astronomy Engine development from different
Debian versions (Buster and Bullseye). Buster installs
doxygen version 1.8.13, but Bullseye installs version 1.9.1.
These two versions of doxygen generate slightly different output
for function pointer typedefs: the older version adds an extra
space between the '*' and the defined type name.
I need the output to be exactly the same so that
the continuous integration tests don't see any changed
files in git after they finish running.
So I added an extra step in hydrogen.js (the code I wrote
that converts the doxygen output into markdown) to squash
multiple contiguous spaces into a single space in the
typedef output.
Implemented the C function Astronomy_ObserverGravity.
It implements the WGS 84 Ellipsoidal Gravity Formula,
yielding the effective observed gravitational acceleration
at a location on or above the Earth's surface.
Wrote a demo program that also serves as a unit test.
I verified a few of the calculations, so the file
demo/c/test/gravity_correct.txt also serves as correct
unit test output.
Ported the C version of BaryState to JavaScript.
Fixed an issue in both the C and JS unit tests:
the JPL Horizons data is given in terms of TT, not UT.
I updated the C aberration unit test to use the barycentric
velocity of the Earth to adjust the apparent position of
a star. This brought the error compared to JPL Horizons
data down from 20.5+ arcseconds to less than 0.4 arcseconds.
Success!
Implemented the C function Astronomy_BaryState().
Used JPL Horizons to generate some test data.
Started work on the C unit test for BaryState,
but it is not yet finished. This is just a good
checkpoint for this work in progress.
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.
Added the following C functions:
Astronomy_Rotation_EQJ_GAL
Astronomy_Rotation_GAL_EQJ
These return rotation matrices to convert between
the galactic and J2000 equatorial orientation systems.
This function is a generalization of Astronomy_LongitudeFromSun,
which it replaces. It calculates the relative ecliptic longitude of one body
with respect to another body, as seen from the Earth.
After implementing the same function in C#, JavaScript,
and Python, I will come back and create a generalized
search algorithm to find the next time two bodies are
at a given apparent relative longitude. Even though this
is a generalization of SearchRelativeLongitude, I will have
to figure out a more general way of tuning the search.
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.
Now there are constants for the mean radii of Jupiter's
four major moons available in the C, C#, Python, and JavaScript
versions of Astronomy Engine.
Clarified that these are all mean radii.
Fixed some lingering "//" comments in the C code
(I want to keep ANSI C code as portable as possible.)
Now callers can create time objects from either UT (UT1/UTC civil time)
or ephemeris/dynamical Terrestrial Time (TT). The new TT functions
numerically solve to find the UT that produces the given TT based
on the Delta-T value at that UT. This is always a very fast
numerical convergence, because TT and UT are almost perfectly
linear over brief time windows.
I increased the error tolerance slightly for the Jupiter moons model.
This shrank the model tables significantly, giving me some more
breathing room to stay under 100K download size.
