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.
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.
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.
Use a private enumerated type to select which direction
the precession and nutation is to be done:
- from date to J2000
- from J2000 to date
Normalize the order of parameters to be consistent
between precession() and nutation(), and across languages.
Pass in AstroTime instead of a pair of floating point TT
values (one of which had to be 0).
Added TypeScript version of ObserverVector(),
but it has not yet been documented or tested.
The C functions that took a parameter of a pointer type
'astro_time_t *' were causing incorrect Markdown to be generated.
Now my custom Markdown translator (hydrogen.js) handles this case.
This function calculates the position of an observer on or
near the surface of the Earth (the geoid) in one of two
equatorial coordinate systems: J2000 or equator-of-date.
Moved the following constant definitions from astronomy.c
to astronomy.h, so external code can use them:
DEG2RAD
RAD2DEG
KM_PER_AU
My custom doxygen-to-markdown translator (hydrogen.js)
now emits markdown for the above constants.
Eliminated the obsolete constants MIN_YEAR and MAX_YEAR.
Astronomy Engine is no longer limited to calculating planets
within that range of years.
Fixed a couple of minor documentation issues in the C code.
Started work on a new function Astronomy_ObserverVector,
but it is just a stub for now.
Astronomy Engine used to use USNO historical and predictive tables,
along with linear interpolation, to calculate Delta-T values.
The problem with the USNO tables is, they did not work well outside
a few centuries around present day.
Later I replaced with Espenak & Meeus piecewise polynomials
that work over a much larger time span (thousands of years).
I just discovered there were still comments in the code referring
to the USNO models. I updated the ones I could find to reflect
the current truth about how the code works today.
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.
This caused me to discover I had forgotten to finish
making the necessary changes to astronomy.ts for saving
the cartesian vector inside the EquatorialCoordinates class.
I also realized I had made a mistake in the documentation
for the y-coordinate of the vector: it is the June solstice;
there is no such thing as a September solstice!
Also fixed some mistakes in demo tests: if something failed,
I was printing out the wrong filename (camera.c instead of camera.cs).
Added C function Astronomy_Pivot to transform a rotation matrix
by rotating it around one of its coordinate axes by a given angle.
Added C function Astronomy_IdentityMatrix that just returns
an identity matrix that can be used as the starting point in
a series of transforms.
C function Astronomy_Equator now also returns the topocentric
equatorial location in the form of a cartesian vector.
This is in a new member of the astro_equatorial_t struct
called 'vec'.
The unit test in ctest.c "Rotation_Pivot()" could be improved
with more and better tests.
Created a demo program camera.c that illustrates using
Astronomy_Pivot() to help calculate the tilt of the sunlit
side of the Moon, as seen by a camera pointing right at it.
The resulting tilt angle is not yet verified.
I need to have some confirmation that it is correct before
porting to the other languages.
I'm about to start working on adding a new output
from the Horizon functions. It was a good time to better
document the ideas behind these calculations, before
adding anything new. These are internal comments only
and do not affect generated documentation.
While I was in there, I noticed extra code that was
checking for impossible return values from atan2().
I eliminated these.
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".
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.
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.
Also reduced stack usage by sharing bary[5] more.
Eliminated a redundant call to MajorBodyBary() when
calculating one-way outside the bounds of PlutoStateTable.
There were two places where I was searching for table entries.
I realized I can just directly calculate the same indexes.
This makes the code smaller and faster.
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.
This is a happy compromise between speed and accuracy.
Solidified this as a #define PLUTO_DT, which will allow
me to do compile-time math to establish an array size for
an internal buffer of lazy-computed dt increments between
most recently requested segment endpoints.
Now we are getting somewhere!
Using the mean acceleration over the interval seems so much
better than the complicated calculus thing I did.
Now I am using a very large time step: 500 days,
with an error less than 0.2 arcminutes.
C PlutoCheck: dt = 500.000000
C PlutoCheck: 2049-12-19T12:00:00.000Z = 18250.000000 UT = 18250.001076 TT
C PlutoCheck: calc pos = [ 37.4373396436339121, -10.2445033319559062, -14.4764587025278448]
C PlutoCheck: ref pos = [ 37.4377303529113306, -10.2466292445590774, -14.4773101309873091]
C PlutoCheck: del pos = [ -0.0003907092774185, 0.0021259126031712, 0.0008514284594643]
C PlutoCheck: diff = 2.323163e-03 AU, 0.198 arcmin
Allow passing in the dt value as an environment variable PLUTO_DT:
$ PLUTO_DT=100 ./ctest pluto
C PlutoCheck: dt = 100.000000
C PlutoCheck: 2049-12-19T12:00:00.000Z = 18250.000000 UT = 18250.001076 TT
C PlutoCheck: calc pos = [ 37.4506659151144774, -10.2638146514621269, -14.4865003148403595]
C PlutoCheck: ref pos = [ 37.4377303529113306, -10.2466292445590774, -14.4773101309873091]
C PlutoCheck: del pos = [ 0.0129355622031468, -0.0171854069030495, -0.0091901838530504]
C PlutoCheck: error = 2.339073e-02 AU, 1.989 arcmin
Default to dt = 50 days.
I had a bug in the calculation of the Solar System Barycenter.
The position and velocity vectors were both backwards.
This caused the barycentric Sun to be in the wrong place.
C PlutoCheck: 2049-12-19T12:00:00.000Z = 18250.000000 UT = 18250.001076 TT
C PlutoCheck: calc pos = [ 37.4386594210222654, -10.2474739785185811, -14.4777836541558340]
C PlutoCheck: ref pos = [ 37.4377303529113306, -10.2466292445590774, -14.4773101309873091]
C PlutoCheck: del pos = [ 0.0009290681109348, -0.0008447339595037, -0.0004735231685249]
C PlutoCheck: error = 1.342001e-03
I didn't implement the math I had derived on paper exactly right.
There was one place where I meant to have (dt^3 / 6), but I had (dt^2 / 6).
This has only a tiny effect on the error. I know this can work better,
so I'm still searching for the real problem.
C PlutoCheck: 2049-12-19T12:00:00.000Z = 18250.000000 UT = 18250.001076 TT
C PlutoCheck: calc pos = [ 37.3682426267669996, -10.0216432448322834, -14.3724661626442582]
C PlutoCheck: ref pos = [ 37.4377303529113306, -10.2466292445590774, -14.4773101309873091]
C PlutoCheck: del pos = [ -0.0694877261443310, 0.2249859997267940, 0.1048439683430509]
C PlutoCheck: error = 2.577586e-01
I picked a worst-case starting time that is halfway between
two entries in the PlutoStateTable[]. This requires the most
iterations to calculate. Now I have a better metric for
position accuracy:
C PlutoCheck: 2049-12-19T12:00:00.000Z = 18250.000000 UT = 18250.001076 TT
C PlutoCheck: calc pos = [ 37.3682426469732860, -10.0216428456503461, -14.3724660445065933]
C PlutoCheck: ref pos = [ 37.4377303529113306, -10.2466292445590774, -14.4773101309873091]
C PlutoCheck: del pos = [ -0.0694877059380445, 0.2249863989087313, 0.1048440864807159]
C PlutoCheck: error = 2.577590e-01
Clearly I still have a lot of work to do!
Got 'ctest pluto' error down to 5.675183e-05 AU.
I was adding vectors when I should have been subtracting,
to find heliocentric Pluto from barycentric Pluto.
When finalizing the position of Pluto, I need to convert
barycentric coordinates to heliocentric coordinates.
I do this by adding the barycentric location of the Sun
to the barycentric location of Pluto to obtain the
heliocentric location of Pluto. But I was using the Sun
coordinates from the simulation starting point, not
at the final time.
This decreases the AU error from 1.369e-02 to 1.347e-02.
I'm still looking for the rest of the error.
I wasn't initialzing the acceleration vector in the value returned
by GravSim(). I'm not sure how any of this ever worked!
Found it by using valgrind.
There is no apparent advantage to looping for convergence in GravSim().
Also, it looks like I can use a much larger dt = 10 days.
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.
Fixed some build warnings that occur on various gcc
optimization levels, and only on this version of gcc.
For now, build ctest with fewer optimizations: -O1
instead of -O3. This is because -O3 and -O2 cause
excessive errors in 'ctest diff' of the order
1.0e-9, where I usually get 1.0e-12. I will have to
come back and figure out exactly which optimization(s)
are causing the problem and turn them off specifically.
This also means I need to document the dangerous optimizations
for people who are using the C version of Astronomy Engine.
When 'ctest diff' fails because of excessive numeric error,
print out the two lines of input text that had the worst
numeric error. This really helps on the Raspberry Pi
where memory is at a premium, and it's hard to open the
full output files using vi.
Added Astronomy_FormatTime to the topic index.
Reworded text to avoid "as explained above", because it turns
out the generated documentation does not always put things
in the same order they appear in the source code comments.
It's surprisingly tricky to print a time rounded to the
nearest millisecond, second, or minute using the C code.
I saw a case where positions.c printed '2020-07-09 04:29:60'.
Because printing a date/time is a basic need of an astronomy program,
I added the new function Astronomy_FormatTime to do this.
All the demo programs use this new function, which required
me to update the correct reference output for the unit tests.
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.
The high-frequency wobble in the Pluto position function was bothering me.
Decreased the arcminute error threshold from 1.0 to 0.5, resulting
in a much larger model, but a lot less ripple:
547 [ 78 140 94 115 21 99]
