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.
The aberration unit test that relies on barycentric
velocity calculation for the Earth's geocenter
has been ported from C to JS and shows identical results.
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.
I realize some use cases require adjustments for
stellar aberration. The existing aberration adjustments
are only supplied for calculating planet positions.
Some users will benefit from being able to add/subtract
aberration corrections to arbitrary vectors, including
for star positions.
I have added some JPL Horizons test data to help
validate the aberration functionality I'm about to add.
I created the beginning of a unit test in ctest.c,
but currently there is no aberration correction
implemented, so the test has no error threshold.
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.
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 added this test, but unfortunately I could not figure
out how to make JPL Horizons generate equatorial and galactic
coordinates using the same aberration model. This appears
to introduce an extra 22 arcseconds of error.
Added a sanity check in the unit tests that the functions
Astronomy_Rotation_EQJ_GAL and Astronomy_Rotation_GAL_EQJ
return matrices that really are inverses of each other.
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.
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.
Starting work on support for galatic coordinates.
Generate a test data file using calculations made
by the NOVAS function equ2gal(). Later I will use
this data to verify the conversion functions I
write for Astronomy Engine.
Decreased the minified browser code from 94918 bytes to 94221 bytes.
Did this by using a more efficient encoding of the IAU2000B nutation model:
instead of making {nals:[_], cls:[_]} objects, make lists of lists [[_], [_]].
Ran 'npm audit fix' to resolve some security vulnerabilities
in the developer tools in the 'generate' directory.
None of the vulnerabilities affect the npm package
astronomy-engine, because it has no external dependencies.
The risk was only to developers who run the code generation
tools, not end users. Even then, the risk is minimal because
these tools run with well-defined inputs that are not subject
to external tampering.
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 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.
The test build failed because diffcalc reported a small
discrepancy between the C and C# output.
So I made the threshold more lenient for now.
I want to come back later and figure out if I can get back
to exact agreement between C and C# code.
Told wget not to output rediculous progress bar stuff
that eats thousands of lines of log output.
I ran into a problem recently that was confusing to debug.
It turned out that I was calling fgets() providing a line buffer
that was not long enough for all of the lines in the input file.
This caused the unread portion of the long line to appear as if
it were the beginning of another line, failing the test in
a weird way.
So I replaced all calls to fgets() in ctest.c with a new
wrapper function ReadLine(). It checks for this issue
and immediately aborts with a helpful diagnostic.
It turns out different Node.js versions do math differently,
which caused a Travis CI build failure.
Scale topocentric distance the same way I scale heliocentric distance.
Adjusted diffcalc bash script and diffcalc.bat Windows batch file accordingly.
The differ now prints the final "score" so I'm less likely to make
a mistake spotting the correct maximum difference.
Removed unused variable in ctest.c DiffLine(): maxdiff.
When comparing calculations of body vectors, scale
the size of the difference by the minimum orbital
radius (or typical radius in the case of the Solar
System Barycenter).
This concludes my investigations of discrepancies between
the various language calculations. I have done as much
as I can without implementing my own trig functions,
which is not worth the effort (or the loss of efficiency
in JavaScript).
Scaling the errors relative the measurement units reveals
that the discrepancies are reasonable for the 16-digit
precision one expects from 64-bit floating point numbers.
The worst case is C vs JavaScript, with a scaled error
of about 7.2e-15. I can live with that.
A given amount of error in an angle measured in
sidereal hours is 15 times more important than the
same numeric error in an angle measured in degrees.
Scale angular errors by the range of values they
could take on. Longitude-like angles in degrees
have a range of 360, while latitude-like angles
range over 180 degrees (-90 to +90).
Split out separate Windows batch file diffcalc.bat,
just like I already split out bash script diffcalc.
