Instead of stopping at the first failure, diffcalc
scripts (Linux and Windows) will now continue to run
all tests before reporting pass/fail. This allows
fixing multiple issues all in one GitHub Actions pass.
Also fixed the immediate issue where I needed to
increase tolerances slightly, as a follow-up fix
for nutation formula changes.
While trying to convert ecliptic coordinates from mean
equinox of date to true equinox of date, I ran into excessive
overhead from the IAU2000B nutation model. The fact that it
uses 77 trigonometric terms made the calculations a lot slower.
https://apps.dtic.mil/sti/pdfs/AD1112517.pdf
Page 4 in the above document mentions a shorter series
“NOD version 2” that has 13 terms instead of 77 as used in IAU2000B.
I had not noticed NOD2 before, because it appears only in
the FORTRAN version of NOVAS 3.x, not the C version.
After reading the FORTRAN code, I realized NOD2 is the same
as IAU2000B, only it keeps the first 13 of 77 terms.
The terms are already arranged in descending order of
significance, so it is easy to truncate the series.
Based on this discovery, I realized I could achieve all of
the required accuracy needed for Astronomy Engine by
keeping only the first 5 terms of the nutation series.
This tremendously speeds up nutation calculations while
sacrificing only a couple of arcseconds of accuracy.
It also makes the minified JavaScript code smaller:
Before: 119500 bytes.
After: 116653 bytes.
So that's what I did here. Most of the work was updating
unit tests for accepting slightly different calculation
results.
The nutation formula change did trigger detection of a
lurking bug in the inverse_terra functions, which convert
a geocentric vector into latitude, longitude, and elevation
(i.e. an Observer object). The Newton's Method loop in
this function was not always converging, resulting in
an infinite loop. I fixed that by increasing the
convergence threshold and throwing an exception
if the loop iterates more than 10 times.
I also fixed a couple of bugs in the `demotest` scripts.
The GitHub Actions test comparing C and Kotlin
output failed for the Mac OS platform.
The test failed locally on my Windows 10 system
for the same reason.
Relaxed the constraints on both so they will pass again.
The `horizon` function was taking the sine and cosine
of the right ascension, but treating it as degrees.
This was wrong because right ascension is in sidereal hours.
Fix: multiply hours by 15 to get degrees.
Keeping the unit test that confirmed this bug.
Added comparison of C and Kotlin output to `diffcalc`
and `diffcalc.bat`. This test now passes on Linux.
This is a big milestone!
Implemented Kotlin functions and code generator
for calculating the state vectors of Jupiter's
largest 4 moons.
Added cautionary comments about needing to correct
Jupiter's moons for light travel time.
This is the first pass to get everything needed
for the AstroCheck tests. I tried comparing
C output to Kotlin output, and there are some
serious problems to debug:
$ ./ctest diff 2.8e-16 temp/{c,k}_check.txt
First file: temp/c_check.txt
Second file: temp/k_check.txt
Tolerance = 2.800e-16
lnum a_value b_value factor diff name
FAIL 137746 4.2937184148112564e+01 4.2944101081740065e+01 0.03364 2.327e-04 helio_x
FAIL 373510 1.4197190315274938e+01 1.4193716564905307e+01 0.03364 1.168e-04 helio_y
FAIL 137746 -6.5897675150466091e+00 -6.5929481589493522e+00 0.03364 1.070e-04 helio_z
FAIL 59150 1.8035183339348251e+01 1.8035909197904104e+01 0.01730 1.255e-05 sky_j2000_ra
FAIL 137747 -8.1222057639092533e+00 -8.1250990689970894e+00 0.00556 1.607e-05 sky_j2000_dec
FAIL 137747 4.8436159305823310e+01 4.8441487614058218e+01 0.03481 1.855e-04 sky_j2000_dist
FAIL 322846 8.7596368704201495e+01 2.6760770774700188e+02 0.00278 4.995e-01 sky_hor_az
FAIL 405828 -6.5075824596574279e+01 5.6922941329250996e+01 0.00556 6.778e-01 sky_hor_alt
OK 92717 4.1268347083494783e-03 4.1268347083494774e-03 223.21429 1.936e-16 jm_x
OK 45091 -8.0149190392649894e-03 -8.0149190392649929e-03 79.42812 2.756e-16 jm_y
OK 135377 1.5470777280065808e-03 1.5470777280065804e-03 223.21429 9.680e-17 jm_z
OK 216836 4.5725777238332412e-03 4.5725777238332394e-03 126.58228 2.196e-16 jm_vx
OK 351647 5.1351566793199944e-03 5.1351566793199962e-03 126.58228 2.196e-16 jm_vy
OK 351647 2.5217607180929289e-03 2.5217607180929298e-03 126.58228 1.098e-16 jm_vz
Score = 6.778e-01
ctest(Diff): EXCEEDED ERROR TOLERANCE.
So I'm checking this in as work-in-progress.
Added the rotationAxis, which calculates dynamic orientation
of planet, Sun, and Moon rotation axes. Added the first
unit test that verifies against JPL Horizons data.
Eliminated more redundant time parameters in precession functions.
More cleanup of C# code: I realized the private function
vector2radec was redundant with the public function EquatorFromVector.
Deleted vector2radec.
Two different people are currently helping me get the
build process working on Mac OS. They both ran into different
amounts of comparison error in the calculations for different
langauges. I updated the 'diffcalc' bash script to have
slightly less strict tolerances, so the unit tests pass.
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.
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.
I was able to get the diffcalc test to confirm that
the C and C# algorithms are producing absolutely identical
output values. I just needed to make the C code print its
output in scientific notation with 18 significant figures.
For angular results that are longitude-like, the severity
of the error decreases as the location approaches the
poles of that coordinate system. So scale the longitude
difference by the cosine of the latitude.
For example, in horizontal coordinates, the closer an object
is to the zenith, the less significant an error in its azimuth
becomes. At the zenith, azimuth loses all meaning.
Also, now we pass in the tolerance on the command line, because
there are known issues for each language. I expect C# and C
to be in very close agreement, and I want to know if they start
drifting apart. However, JS has known issues with atan2(), so I'm more
lenient with C vs JS.
Now track and report max difference independently for each
type of calculated value, instead of merging them all
together in one big soup.
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.
