Generate longitude output files that are tested by
the `generate` program.
Discovered that Windows `run.bat` was not correctly testing
the longitude output test files generated by the JavaScript
code. Fixed that.
Implemented searchRelativeLongitude, which finds
planetary conjunctions and oppositions.
Discovered I can make all languages' unit tests
more strict: 6.8 minute error tolerance instead of 15.
Fixed documentation mistake in C# function SearchRelativeLongitude:
the function cannot return null. It either finds a solution time
or throws an exception.
Simplified Kotlin unit tests: use a more compact pattern of
scanning space-delimited tokens in lines.
For years before 1582 or years after 3668, the Seasons functions
were unable to find many equinoxes and/or solstices.
The problem was that over time, the Earth's axis precesses
enough that the calendar dates of these events drifts outside
the fixed search ranges I had provided for them.
I expanded the search ranges so all season changes can be found
for a much wider range of years, as verified by unit tests:
C/C++: -2000..9999
C#: 1..9999
JavaScript: -2000..9999
Python: 1..9999
Kotlin: 1..9999
Note: C#, Python, and Kotlin currently do not allow
years values below +1. In fact, I discovered we were not
noticing when an invalid year was passed into the Kotlin code.
I updated that code to throw an exception when the year does
not match what was expected. It is disturbing that the
GregorianCalendar class silently ignores invalid years!
Constricted the search tolerance from 1 second to 0.01
seconds for the seasons search, to ensure more consistent
behavior.
Fixed a bug in the Kotlin search() function's
quadratic interpolation that was causing the convergence
to be slower than it should have been.
Added Astronomy.searchAltitude, which enables more generic
searches for altitude events. The most common use for these
are finding civil, nautical, and astronomical twilight times.
Implemented the following Kotlin functions:
Astronomy.searchHourAngle
Astronomy.searchRiseSet
The Kotlin code can now search for rise/set
times for a given Earthbound observer for the
Sun, Moon, or any planet.
It can also search for times when a given body
reaches a desired hour angle. This has its own
value (for example, culmination), but is also
used to assist finding time brackets that bound
rise/set events.
Added special case exception type:
EarthNotAllowedException.
This follows the pattern of the other languages,
and makes diagnosing a violation easier than
the more generic InvalidBodyException.
Minor simplifications to the C# function
Astronomy.InternalSearchAltitude.
Improved the comments for the C# unit test
function RiseSetTest. They make the algorithm
easier to understand.
Added the following Kotlin functions:
equatorialToEcliptic
pairLongitude
moonPhase
searchMoonPhase
searchMoonQuarter
nextMoonQuarter
Discovered I could tighten the tolerance for the moon phase
unit tests from 120 seconds to 90 seconds and they still pass.
Implemented Astronomy.siderealTime() in Kotlin.
Updated all languages' unit tests for sidereal time
to verify exact conformity between them, rather than
to an externally derived value. I wanted to make
sure all languages, including Kotlin, are calculating
the exact same value.
I don't need an external authoritative test for
sidereal time, because it will be indirectly tested
through its involvement in thousands of other calculations
that depend on it. I just need a quick sanity check
before implementing those other things that depend on it.
There was already an internal function for calculating
Greenwich Apparent Sidereal Time (GAST). By request,
I have exposed this function for outside users.
Added a minimal unit test to verify the function is
callable and returns the correct result for one case.
This function is already exhaustively tested by unit
tests that verify other functions that already called
this function when it was internal, so minimal testing
is sufficient in this case.
Added the following new functions to all 4 languages:
MassProduct: find the GM product for all Solar System bodies.
LagrangePoint: calculate L1..L5 state vectors for a pair of bodies.
LagrangePointFast: calculate L1..L5 state vectors given
state vectors and GM products of a pair of bodies.
I want to be able to re-use the code for loading state
vectors from a JPL Horizons text file, instead of copy-n-paste
like I did in C. So I reworked it as an iterator.
Now that Microsoft has officially released .NET 6,
I have upgraded the C# version of Astronomy Engine to use it.
No source code changes were needed. I just bumped the
version number in the project files, and targeted .NET 6
in the GitHub Actions continuous integration tests.
Fixed some obsolete wording in generate/README.md.
The reason SSB vector errors were larger than other bodies
is because the Sun/Barycenter relationship does not have
position and velocity vectors of a "typical" size.
The distance of a planet from the SSB is fairly constant,
and the speed a planet travels is fairly constant.
Therefore, comparing errors by dividing by the magnitude
of the correct vector usually makes sense for scaling.
But for the barycentric Sun (or the heliocentric SSB),
the magnitude of the vectors can become arbitrarily small,
nearly zero in fact, resulting in surprisingly large ratios.
I compensated for this in all the tests by adding a new rule.
When the error thresholds r_thresh and v_thresh are negative,
it is a flag that indicates they are absolute, not relative.
In other words, r_thresh < 0 indicates that abs(r_thresh) is
the maximum number of astronomical units allowed in position
errors, and v_thresh < 0 specifies maximum AU/day.
This results in more consistent numbers that give confidence
the errors are indeed very small and not worth worrying about.
In the C unit test for barystate, heliostate, and topostate,
I had switched from checking absolute differences to relative
differences. I forgot to do that in the other 3 languages
until now. They are all working consistently in how they
measure calculation errors.
Implemented the C# version of the ObserverState function.
This returns the geocentric position and velocity for
a point on the Earth's surface at a given time.
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.
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.
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.
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
----------------------------------------------------------------------------------------------------
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 is technically a breaking change, but only for clients
that use the cartesian coordinates in an ecliptic coordinate
return type. Before now, the coordinates were just separate
floating-point members ex, ey, ez. Now they are a standard
vector type.
The purpose is to allow seamless interfacing with vector
rotation functions, and to be consistent with the equatorial
coordinate types.
Now that equatorial coordinates include both angles
and cartesian coordinates, there is no need for the
VectorFromEquator function. It has been removed
from all four supported languages.
The expression "VectorFromEquator(equ, time)" can be
replaced with "equ.vec" in any calling code.
Created new rotation matrix functions for the C# version.
IdentityMatrix creates a new instance of the 3x3 identity matrix
1 0 0
0 1 0
0 0 1
Pivot transforms a rotation matrix by pivoting it about
one of its coordinate axes by a specified angle.
Still need to port the C version of the "camera" demo.
When built from a system with a European (or similar) culture setting
where a comma is used as a decimal marker instead of a period,
the C# unit tests and demos would fail.
Now explicitly specify InvariantCulture to resolve these problems.
In all four versions of Astronomy Engine (C, C#, JavaScript, and Python),
starting a search for a full moon near December 19, 2020 would fail.
I added a unit test to all four languages and it failed consistently
across them all.
The root cause: I was too optimistic about how narrow I could make
the window around the approximate moon phase time in the
SearchMoonPhase functions. Finding the exact moon phase time failed
because it was outside this excessively small window around the approximate
time. I increased the window from 1.8 days to 3.0 days.
This should handle all cases with minimal impact on performance.
Now all four of the new unit tests pass.
