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.
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.
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.
