cosinekitty / astronomy

Astronomy Engine: multi-language calculation of Sun, Moon, and planet positions. Predicts lunar phases, eclipses, transits, oppositions, conjunctions, equinoxes, solstices, rise/set times, and other events. Provides vector and angular coordinate transforms among equatorial, ecliptic, horizontal, and galactic orientations.
MIT License
422 stars 59 forks source link

Predict Jupiter moon shadow and eclipse events #140

Open cosinekitty opened 2 years ago

cosinekitty commented 2 years ago

Here is an example of the kind of events I want to predict. This text comes from the Sky and Telescope Jupiter's Moon calculator. That tool also produces a nice image of the planet and moons. I am thinking about creating a similar rendering tool as a web page.

Wednesday, December 8, 2021

02:02 UT, Ganymede begins transit of Jupiter.
05:40 UT, Ganymede ends transit of Jupiter.
07:12 UT, Ganymede's shadow begins to cross Jupiter.
10:46 UT, Ganymede's shadow leaves Jupiter's disk.
15:36 UT, Io enters occultation behind Jupiter.
19:10 UT, Io exits eclipse by Jupiter's shadow.

Thursday, December 9, 2021

01:32 UT, Europa enters occultation behind Jupiter.
06:58 UT, Europa exits eclipse by Jupiter's shadow.
12:56 UT, Io begins transit of Jupiter.
14:12 UT, Io's shadow begins to cross Jupiter.
15:14 UT, Io ends transit of Jupiter.
16:32 UT, Io's shadow leaves Jupiter's disk.

Friday, December 10, 2021

06:06 UT, Callisto begins transit of Jupiter.
10:06 UT, Io enters occultation behind Jupiter.
10:46 UT, Callisto ends transit of Jupiter.
13:40 UT, Io exits eclipse by Jupiter's shadow.
18:14 UT, Callisto's shadow begins to cross Jupiter.
19:44 UT, Europa begins transit of Jupiter.
22:14 UT, Europa's shadow begins to cross Jupiter.
22:30 UT, Callisto's shadow leaves Jupiter's disk.
22:36 UT, Europa ends transit of Jupiter.

Saturday, December 11, 2021

01:08 UT, Europa's shadow leaves Jupiter's disk.
07:26 UT, Io begins transit of Jupiter.
08:42 UT, Io's shadow begins to cross Jupiter.
09:44 UT, Io ends transit of Jupiter.
11:00 UT, Io's shadow leaves Jupiter's disk.
16:06 UT, Ganymede enters occultation behind Jupiter.
19:44 UT, Ganymede exits occultation behind Jupiter.
21:08 UT, Ganymede enters eclipse by Jupiter's shadow.
prideout commented 2 years ago

This might be especially interesting since the 2006 analytical model for Jupiter moons that you are using might be more precise than what the S&T app uses.

By the way, something I discovered with your library that I find interesting is that the orbital period of the Galilean moons seems to vary. I wonder if that's due to LaPlace resonance. I know that their eccentricity varies due to resonance, but I didn't realize their periods vary too.

cosinekitty commented 2 years ago

@prideout Yes, it's possible the IMCCE formulas could be more accurate. I'm not sure what Sky & Telescope uses for their predictions, but they have been doing it for a long time. On the other hand, I truncate the IMCCE series to make the code smaller and faster, so it might be less accurate! I omit as many terms as possible while maintaining relative errors (position error divided by the satellite's distance from the center of Jupiter, and the velocity error divided by the satellite's Jovicentric speed) within 3.5e-7. This is especially important for the JavaScript version, where I value small downloads much more than the tiny residual improvement in accuracy.

I'm not sure about the physical cause of the orbital periods changing, but it would make sense that it is a resonance effect of the moons on each other. It makes sense from the mathematical model's point of view, because the IMCCE series are representing the orbits as ever-changing ellipses. All 6 Keplerian elements are functions of time in their model. Especially when you allow the semi-major axis value to change with time, it implies a change in orbital period.

I believe their models are "semi-analytic", meaning they are based on state-of-the-art numerical integrations (a direct gravity simulation) that is then fed into a best-fit model to produce the series coefficients. So it is possible to construct such a model without even thinking about resonance effects; any such effects will be captured by the original integration. I know, for example, that the VSOP87 model of the planets was based on NASA JPL numerical integrations.

prideout commented 2 years ago

I'm using the Python variant of your library rather than the JS, does that have the truncation that you describe? I found a paper that has graphs of Io's changing eccentricity and it seems roughly consistent with what I'm seeing with your library. It's hard to find literature about this effect that's actually readable to amateurs like me.

By the way I've authored some poor man's eclipse detection code. My strategy was to do point-line distance, where the point is Io's center, and the line goes from the Sun to Jupiter. However this approach approximates the shadow with a cylinder instead of a cone, so maybe you'll think of something smarter. Here it is for reference:

def is_totally_eclipsed():
    distance_from_line = io_position.cross(jupiter_vector).Length()
    within_shadow_cylinder = distance_from_line < (JUPITER_RADIUS - IO_RADIUS)
    far_side = io_position.Length() > jupiter_position.Length()
    return within_shadow_cylinder and far_side
cosinekitty commented 2 years ago

Yes, all 4 languages do the same exact numeric calculation, including the truncation. As much as possible, I keep them all numerically identical. The truncation also makes the math faster, which is especially important for the Python version.

I think your approximation of the shadow as a cylinder is perfectly fine, because of how far away Jupiter is from the Sun, and how close the moons are to Jupiter. The distance ratios make the cone angles negligible.

If you test this against S&T, you will probably find a plus or minus 8 minute error however, because of a plus or minus 1AU change the in the light arrival time as seen from the Earth. That is, unless you already thought of correcting for light travel time! (I can't tell from this snippet of code.)

You can take a look at what I did in an experimental image raytracer I wrote in C++: https://github.com/cosinekitty/astronomy/blob/master/demo/c/raytrace/main.cpp The function PlanetImage at line 316 calls Astronomy_GeoVector, which knows to correct for light travel time to find the apparent position of Jupiter as seen from the Earth. Then at line 332 it figures out how far into the past we are seeing the Jupiter system, and calculates the positions of Jupiter's moons at that time in the past.

That way you get an accurate prediction of what an observer on the Earth will see through a telescope.

offsky commented 2 years ago

Just one idea about truncating precision to optimize for download size. Would it be possible to give users of your library a choice? Maybe add an optional build product that is 2x the size and includes increased precision. For my uses, I dont need the extra accuracy, but maybe someone does.

cosinekitty commented 2 years ago

Hi @offsky, if anyone needs higher precision, I would help them adjust things so they could rebuild the code generator to do less truncation. However, truncated series are only part of the work. There are many places where I took shortcuts, sacrificing an arcsecond or two of accuracy for much smaller code. Just one example: I ignore relativistic corrections for the bending of light rays as they pass near the Sun and planets, and time dilation effects. There are other such things I leave out that would otherwise bloat the code, but only add fractional improvements to accuracy. The big burden would be going through and adjusting the unit tests to validate that the accuracy has actually improved.

I believe applications that really need extreme accuracy should use something like NOVAS C 3.1, and bear the burden of storing and accessing the large ephemeris files. Astronomy Engine is trying to fill a niche where you can actually do arcminute-precision calculations in a browser with a very small download of minified JS, and with very little CPU overhead. My mission is to serve the middle ground that I saw was missing: everything I found for astronomy calculation was either (a) too simple and woefully inaccurate or (b) hyper-aggressively accurate at the expense of being unusable in a browser or microcontroller environment.

prideout commented 2 years ago

This is a little GIF I made using this library combined with skia-python for rendering. It shows the path of Jupiter and Io through space, while staying fixed on Jupiter. Io's radius is exaggerated but otherwise this is drawn to scale. The portion of Io's orbit that is eclipsed is drawn in black. To detect the eclipse, I used something similar to the code snippet I posted earlier.

challenge-800k

cosinekitty commented 2 years ago

@prideout That is really cool!

prideout commented 2 years ago

Sorry for the spam but I just want to mention that the Search function that you provide seems to be working with my distance-from-line method. For the signed function that needs to cross 0, I simply return JUPITER_RADIUS - IO_RADIUS - distance_from_line. It's ok to ignore the far_side part of my function because I only call Search after I've established rough time bounds.

I'm not familiar with the algorithm that Search uses, is it a convex optimization technique?

cosinekitty commented 2 years ago

Not spam at all! I really enjoy the social/intellectual interaction with like-minded people. That is great to hear that Search works for your use case. This algorithm is something I mashed together, and I tuned it to work well with the kind of functions I wanted to find ascending roots for. It fits the search region with a best-fit quadratic curve and solves for the root. If the root of the quadratic is close enough to a tighter-fit range that contains the actual root, it shrinks the search range dramatically. Otherwise, it uses simple bisection for the current iteration, picking whichever half of the range contains the root (by checking for positive/negative endpoints and midpoint). It keeps looping until the caller-specified time tolerance is met. Most of the astronomy functions are very smooth, and a parabola is a very good approximation, so it works well enough for my needs.