I suddenly realized today that I could eliminate lots of red wavy
lines in the source editor when working on generate/template/astronomy.*
files, by tweaking the $ASTRO_... tokens to look like comments.
For some reason, I didn't think of doing this before now.
I was motivated by the Go language support, where the tokens were
breaking `go fmt`.
On the Raspberry Pi 4, using latest versions of cppcheck
and pylint, a few more minor fixes were needed for eliminating
warnings.
Also had to soften a tolerance for the Kotlin unit tests.
The generated code for the Pluto state table in Python
now uses a class `_pstate` for better type checking.
It also makes the code easier to understand.
Moved class _TerseVector higher in the source file to
reduce the need for quoted forward type declarations.
I added the mypy option `--disallow-untyped-defs` to fail
any function lacking complete type hints.
Then I fixed all the resulting errors.
I ended up changing the Python code generator to create
some tuple types instead of list, because it is possible
to write stricter type checks that way. This was in
the Pluto and Jupiter Moon tables.
I still should come back and do the same thing for the VSOP tables.
The type checking revealed a couple of places where I wasn't
checking for a search failure. I fixed those too.
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.
Modified the code generator and source templates to allow
using fewer than 77 terms from the IAU2000B nutation model.
The number of terms causes calculations to be far slower than
I would like, plus most of these terms provide far less than
one arcsecond of difference in the output.
I want to experiment with truncated versions of the series.
https://apps.dtic.mil/sti/pdfs/AD1112517.pdf
Page 4 in the above document references a shorter series
"NOD version 2" that has only 13 terms instead of 77 as used in IAU2000B.
I found that code and confirmed the 13 terms are the first 13
consecutive terms from the original 77.
This commit does not change the number of terms calculated;
it only enables doing so by changing a single #define in
generate/codegen.c.
Once I change the nutation model, I'm sure there will be multiple
tweaks to unit tests to get everything working again.
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 calculation of heliocentric and barycentric
state vectors for Pluto. This is done using
a gravity simulator that treats Pluto as a negligible
mass that is affected by the major masses of the
Solar System: Sun, Jupiter, Saturn, Uranus, Neptune.
Updated the code generator to write the Kotlin
version of the Pluto state table, a lookup table
of known-correct state vectors of Pluto at long
intervals, derived from the TOP2013 model.
The gravity simulator interpolates state vectors
of Pluto between these known-correct states.
Minor code style cleanup in the Kotlin source.
Fixed a possible thread-safety issue in the C# code.
Implemented the Montenbruck/Pfleger version of the NAO1954
geocentric Moon model in Kotlin.
In the process, code review helped simplify parts of the C# code.
I may want to go back and see if I can simplify the CalcMoon
code in the other languages too.
Added code generator for VSOP87 major planet models in Kotlin.
Added a lookup function for VSOP87 models to Kotlin template.
Moved IAU2000b delcaration outside class Astronomy.
Implemented the iau2000b nutation formula.
Implemented calculation of the precession
rotation matrix.
This is the very first use of a code generator
macro for the Kotlin code. I am going to try
keeping all these macros toward the bottom
of the source tepmlate, so that as I look at
line numbers for compiler errors, they will match
between the target code and the template code.
I may go back and rework the other languges to
do this also. I'm not sure why I didn't think
of this before!
The code generator was creating slightly different numeric
values for the Pluto state tables and the Jupiter rotation matrix.
I decreased the output precision by one decimal digit.
This should allow the code generator to produce identical
source code on both Linux and macOS.
The Pluto gravity simulator constants now come from
a single source: pluto_gravsim.h. This will allow me
to experiment with the Pluto state table to get a better
compromise between size and accuracy.
Reworked the Pluto gravity sim constants so they are defined
in one place: a new header file gravsim/pluto_gravsim.h.
Then the code generator writes the #defines to the C code, instead
of having two independent versions of the same constants.
I will continue down the road of having a single-source-of-truth
for these constants across all 4 supported languages.
Also, confusingly, I had one constant called PLUTO_DT in codegen.c
that was called PLUTO_TIME_STEP in astronomy.c. Also, astronomy.c
had a different constant PLUTO_DT that didn't mean the same thing.
I reworked the naming to be consistent in all places.
I already had a TopPosition() function that knows how to calculate
exact equatorial coordinates, so I eliminated the redundant logic
from gravsim_test.c
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 [[_], [_]].
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.
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.
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.
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.
Improved the type checking by using tsc --strict.
Nothing substatial changed in the generated JavaScript, and no
actual bugs were found, but I removed a lot of loose/sloppy
type signatures. This should make mistakes less likely
in the JavaScript code going forward.
I forgot that my build process automatically updates
copyright years when the current year changes.
My Travis CI unit tests verify that there are no local
changes after running all the tests.
That test failed because the update_copyrights.py changed
all the "2019-2020" to "2019-2021".
I believe this wraps up the Python integrator.
It now works in all 4 languages and passes all tests.
Fixed up demo tests to match new output.
Turned on Travis CI checking in this branch again.
Ported Pluto integrator to C#.
Along the way, I noticed that I had VSOP87 latitude and longitude
swapped in such a way that they worked, but were labeled wrong.
This confused me quite a bit as I tried to implement functions
to calculate the derivatives of the VSOP87 spherical coordinates.
Fixed this in the code generator and the C and C# template files.
The PlutoStateTable was slightly different when generated
in Linux and Windows, because there was so much precision
after the decimal point. Reduced precision until Linux
and Windows generated the exact same output.
Instead of remembering the most recent 3 segments of Pluto's orbits,
cache up to all 40 segments within the year span 0000..4000.
Use dynamic memory to allocate them instead of static memory.
Added Astronomy_Reset() to free memory before exit.
To make Pluto calculations have a low amortized time cost,
calculate up to 3 segments between adjacent pairs of
pre-calulcated states and recycle them. Do linear ramp
fade-mixing between them.
Currently, something is wrong with this because it fails
by inaccurately calculating horizontal coordinates of Pluto
in one test. I'm not sure what's going wrong there yet,
but likely related to multiple calls via search.
Added PlutoStateTable to C code generator.
This is a table of known correct [tt, pos, vel] tuples for Pluto,
calculated using TOP2013. These serve as seed points from which
to integrate Pluto's motion.
Added PlutoCheck() function to ctest, just to get going.
I have a lot more peformance work to do in order to make
the full blown unit test to finish in a reasonable amount of
time.
Changes in astronomy.c:
Added some generic "terse vector" support.
A terse vector contains 3 components and nothing else.
This is handy for implementing compact formulas for various
vector expressions.
Created enhanced VSOP87 calculations that provide
velocity vectors as well as position vectors for
the Sun, Jupiter, Saturn, Uranus, and Neptune.
These are the "major" bodies that have significant
effects on the motion of Pluto. Also used to convert
heliocentric coordinates to barycentric coordinates.
Implemented first version of the integrator logic.
It is not accurate enough yet, and it is far too slow.
I need to debug the accuracy first, then I will work on
making it faster.
I was experimenting with an alternative way to calculate Pluto's position.
That idea didn't work out, but in the process, I did make some changes
I want to keep:
1. A typecast in generate/codegen.c that eliminates a build warning
in Windows (Visual Studio 2015).
2. Expanded generate/vsop/vsop.c to calculate both position and velocity
vectors. This was tricky to get working, and could come in handy.
Version 9.3.0 of gcc has extra warnings that detect snprintf
truncations and use of possibly uninitialized variables.
Fixed some warnings of both types.
Also corrected code generator to output term coefficients
in scientific notation. In the C code, it was dropping signficant
digits by outputting in fixed point notation.
