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.
All along, ctest(diff) has been treating the "v" lines as
v x y z
instead of
v tt x y z
and ignoring the z values. Fixed that. It was a relief
it didn't reveal any lurking problems in the z values.
Also, dramatically decreased the disagreement measured between the
C and C# code just by printing out more decimal places in both.
First file: temp/c_check.txt
Second file: dotnet/csharp_test/csharp_check.txt
Tolerance = 5.600000e-17
lnum a_value b_value factor diff name
OK 9269 4.3776373521942023e-05 4.3776373521942017e-05 1.00000 6.776e-21 helio_x
OK 5589 -4.3623354680948618e-05 -4.3623354680948625e-05 1.00000 6.776e-21 helio_y
OK 606 -4.3091817463880202e-05 -4.3091817463880195e-05 1.00000 6.776e-21 helio_z
OK 199595 -2.5752073258047489e-05 -2.5752073258047493e-05 1.00000 3.388e-21 sky_j2000_dec
----------------------------------------------------------------------------------------------------
First file: temp/c_check.txt
Second file: temp/js_check.txt
Tolerance = 1.200000e-12
lnum a_value b_value factor diff name
OK 388830 1.0342029871624860e+01 1.0342029871624890e+01 1.00000 3.020e-14 helio_x
OK 67200 -2.9057622338509663e+00 -2.9057622338509423e+00 1.00000 2.398e-14 helio_y
OK 67200 -4.3874295048838763e-01 -4.3874295048837719e-01 1.00000 1.044e-14 helio_z
OK 282864 4.2562907944628900e+00 4.2562907944628821e+00 0.94021 7.516e-15 sky_j2000_ra
OK 333852 -9.3967798183388087e+00 -9.3967798183388513e+00 1.00000 4.263e-14 sky_j2000_dec
OK 151220 3.0052362065237119e+01 3.0052362065237109e+01 1.00000 1.066e-14 sky_j2000_dist
OK 123271 8.3934732121148897e+01 8.3934732121147533e+01 0.83370 1.137e-12 sky_hor_az
OK 123271 -3.3519630572770836e+01 -3.3519630572771831e+01 1.00000 9.948e-13 sky_hor_alt
----------------------------------------------------------------------------------------------------
First file: temp/c_check.txt
Second file: temp/py_check.txt
Tolerance = 5.100000e-14
lnum a_value b_value factor diff name
OK 22 -2.5294053992874876e-01 -2.5294053992874882e-01 1.00000 5.551e-17 helio_x
OK 29 -3.3385318243368106e-01 -3.3385318243368112e-01 1.00000 5.551e-17 helio_y
OK 22 3.7850398890096215e-01 3.7850398890096221e-01 1.00000 5.551e-17 helio_z
OK 135350 2.3433972116724325e-01 2.3433972116724319e-01 1.00000 5.551e-17 sky_j2000_ra
OK 603 4.5414554731189877e-01 4.5414554731189882e-01 1.00000 5.551e-17 sky_j2000_dec
OK 490 4.2751511640162104e-01 4.2751511640162099e-01 1.00000 5.551e-17 sky_j2000_dist
OK 49066 3.2035956679701377e+02 3.2035956679701371e+02 0.88694 5.042e-14 sky_hor_az
OK 49066 -2.7508172411136329e+01 -2.7508172411136300e+01 1.00000 2.842e-14 sky_hor_alt
----------------------------------------------------------------------------------------------------
First file: temp/js_check.txt
Second file: temp/py_check.txt
Tolerance = 1.200000e-12
lnum a_value b_value factor diff name
OK 388830 1.0342029871624890e+01 1.0342029871624860e+01 1.00000 3.020e-14 helio_x
OK 67200 -2.9057622338509423e+00 -2.9057622338509663e+00 1.00000 2.398e-14 helio_y
OK 67200 -4.3874295048837719e-01 -4.3874295048838757e-01 1.00000 1.038e-14 helio_z
OK 282864 4.2562907944628821e+00 4.2562907944628900e+00 0.94021 7.516e-15 sky_j2000_ra
OK 333852 -9.3967798183388513e+00 -9.3967798183388087e+00 1.00000 4.263e-14 sky_j2000_dec
OK 151220 3.0052362065237109e+01 3.0052362065237119e+01 1.00000 1.066e-14 sky_j2000_dist
OK 123271 8.3934732121147533e+01 8.3934732121148897e+01 0.83370 1.137e-12 sky_hor_az
OK 123271 -3.3519630572771831e+01 -3.3519630572770836e+01 1.00000 9.948e-13 sky_hor_alt
----------------------------------------------------------------------------------------------------
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.
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.
The demo shows how to correct for light travel
time to render Jupiter's moons as they appear
from the Earth.
Created an addition operator for the Vector
class in the Python code, because it is handy.
Corrected a bug in the string representation
of the Python StateVector class.
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.)
To assist software that wants to depict Jupiter and its 4 major moons
as they would appear in a telescope, it is important to know their
physical sizes. I already had constants for Jupiter's equatorial
and polar radii. Here I add constants for the radii of the moons
Io, Europa, Ganymede, and Callisto. They are all nearly spherical,
so a single mean radius value is sufficient.
My pydown.py custom Markdown generator was printing bogus
warnings about unknown symbol types, when it was actually
generating correct documentation for those symbols.
Eliminated the warnings, and improved the output format
for global constant documentation: no more extraneous spaces.
If there really is an undocumented symbol detected, fail the build!
Don't just print a warning that slides up the screen unnoticed.
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.
I don't like how close I am to my 100K target size, now
that I'm calculating Jupiter's moons.
Simplified the spin() function so its minified code is smaller.
I will look for other things I can shrink too.
The optimizer makes the Jupiter moons series as short as
possible while keeping error within an acceptable limit.
This should help produce much smaller code, especially
for JavaScript where it really matters.
I want to experiment with truncating the L1.2 series to
sacrifice some accuracy for smaller generated code.
To that end, I implemented the ability to save the
Jupiter moons model after loading it. I added a 'jmopt'
command to the 'generate' program that will do this
optimization. For now, it just loads the model and
saves it back to a different file. Then the code generator
loads from the saved file instead of the original.
This commit verifies that everything is still working,
before I start truncating the series.
I want to make sure I have thorough coverage of exercising
the Jupiter moon calculations before I start tinkering
with truncating the L1.2 models; the eventual goal is
to decrease the code size, especially later for JavaScript.
Only 101 test cases seemed far too small. Now there are 5001.
Increased the time range to cover the years 1931..2068.
Decreased the time step from 100 days to 10 days.
Eliminated all but the geometric data, because I'm testing
code without any light-time or aberration correction.
Updated the README.txt instructions for generating the
JPL Horizons test data.
Output the Jupiter moon model data tables in a tidier format.
Format the amplitudes as fixed-point instead of exponential,
so that the JavaScript minifier will have an easier time
shrinking the data (later, when I get to the JavaScript version).
I translated the L1.2 FORTRAN code into C, and verified
that the calculations match the Stellarium code I modified
to produce EQJ coordinates. I still need to compare against
JPL Horizons data.
This is the beginning of a unit test for the new C function
Astronomy_JupiterMoons(). It reads the stellarium file and
can iteratively solve for the ut corresponding to a given tt.
The test "passes" because it doesn't actually check the value
returned by Astronomy_JupiterMoons() yet.
Instead of calculating ECL coordinates that I will later
have to convert to EQJ, I re-ran all the JPL Horizons test
data sets for EQJ, and updated the rotation matrix in the
Stellarium sample code to generate EQJ output. I'm still
getting reasonably good fit between the two: the max
error is about 1 part in 5000. I was hoping for better,
so I still wonder if I'm just a tiny bit off in some respect.
Stellarium uses the same L1.2 model for Jupiter's moons that
I am implementing. I wanted to confirm that I have valid test
data that matches something authoritative.
Work in progress.
Generating the data tables for Jupiter's moons, but not using them yet.
Created a stub function Astronomy_JupiterMoons(), but it just
returns invalid vectors. The formulas have not yet been implemented.
I found some FORTRAN code for calculating the positions of Jupiter's
moons. I recorded its origin in README.md. I'm going to experiment
with using it as a basis for doing these calculations in Astronomy Engine.
I am starting the process of implementing calculation
of Jupiter's four largest moons: Io, Europa, Ganymede, Callisto.
This commit just contains constant declarations for the
equatorial, polar, and volumetric mean radii of Jupiter.
The positions of the moons will be related to the center
of Jupiter and be expressed in Jupiter equatorial radius units,
so I felt it would be good to give users a way to convert to
kilometers, which can in turn be converted to AU.
Automatically update the front page README.md to include the current
byte size of astronomy.browser.min.js. Fail the build process if
this file ever grows to 100000 bytes or larger.
Python docstrings don't work for variables, so I hacked
a special comment format for helping pydown generate Markdown
text for the README.md for the exported constant KM_PER_AU,
or any other constants I may want to expose in the future.
