Plotting Galilean moon positions per time?

[slegeza]

Hi,

I was wondering if there exists a system with Skyfield that plots the positions of Jupiter's moons as a function of time? It would be very useful for people intending to photograph Jupiter, and it makes it a lot easier to identify which moons are which.

[brandon-rhodes]

Not currently. You can yourself compute the ecliptic positions of Jupiter and one of its moons and subtract to find how far off of Jupiter's center that moon is, but there is no helper function. I'll keep this issue open, though, with the idea of someday adding one!

[brandon-rhodes]

It turns out that Jovian moons are a big topic because users will probably not only want their position, but data on which moons are casting a shadow on Jupiter, are in Jupiter's shadow themselves, are behind Jupiter from Earth's point of view, and are in front of Jupiter from Earth's point of view. Until I can tackle all of those circumstances (for which the Explanatory Supplement provides no approach — so I'll have to make one up on my own), here's a simple routine I've written that will give you an x and y relative to Jupiter for each of the four popular satellites:

def jovian_moons(ephemeris, center):
    """Compute the relative positions to Jupiter of its 4 major moons.

    Given both an ephemeris and the Solar System position of a
    ``center`` from which to observe Jupiter (usually the result of
    calling ``earth.at(t)``), return a list of ``(name, x, y)`` tuples
    each giving the name of a moon and its angular separation from
    Jupiter in arcseconds.  Positive arcseconds are east of Jupiter
    while negative values indicate a moon that is trailing to the west.

    """
    t = center.t

    j = ephemeris['JUPITER BARYCENTER']
    lat, lon, distance = center.observe(j).ecliptic_latlon(t)
    jlon_arcseconds = lon.arcseconds

    names = ('Io', 'Europa', 'Ganymede', 'Callisto')
    result = []
    for name in names:
        moon = ephemeris[name]
        lat, lon, distance = center.observe(moon).ecliptic_latlon(t)
        result.append((name, lon.arcseconds - jlon_arcseconds))
    return result

Note that you'll need the ephemeris jup310.bsp.

I'll keep this issue open until I have time to do more work on it!

[slegeza]

It should also factor in the apparent radius of Jupiter to tell the user if one of the moons is eclipsing, which also means it'll have to calculate the variance in declination for the moons' positions (looks like you already have that). The shadowing will probably be a little more complicated, having to compensate for Jupiter's current orbital position around the Sun. If you'd like, I can submit a pull request and take a stab at it when I'm available!

[ghost]

You might want to look at CSPICE for existing code. I haven't looked deeply, but if you visit:

https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/cspice/index.html

and search for the word "occult", you'll see several helpful functions.

I'm almost sure you can use CSPICE to compute Jupiter's moons passing behind or in front of Jupiter.

[psbaltar]

I took a stab at calculating moon shadow positions on Jupiter. What I came up with is close, but not exact. I couldn't find any authoritative sources for the positions of the shadows, so I eyeballed it in Stellarium.

I'm projecting the shadow onto a plane that's perpendicular to the shadow and intersects Jupiter's center. This gives me the shadow's position relative to Jupiter's center (and assumes that the planet was flattened or cut in half). I then convert this back to ICRF and calculate where the shadow is when viewed from Earth.

The obvious problem is that Jupiter isn't actually flat. Though, because of the way the angles work out, I expect that this method is always fairly close for Jupiter and its moons. It should work best close to opposition since we're viewing it along the same direction as the shadow.

from skyfield.api import load
from skyfield.functions import mxm, mxv, length_of
from skyfield.positionlib import ICRF

import numpy as np

# calculate a reference frame that has a plane perpendicular to shadow vector
# - x-axis is parallel to moon->sun vector (sv)
# - shadow is projected onto YZ plane
def calcR(sv):
  x_ = sv # x-axis is parallel to shadow vector
  x = x_/length_of(x_) # get a unit vector
  z_ = np.array([0,0,1]) # dummy z to calc y-axis that's orthogonal to x-axis
  y_ = np.cross(z_, x)
  y = y_/length_of(y_)
  z = np.cross(x, y)

  R = np.array([x, y, z]) # ICRF to "projection plane" rotation matrix
  return R

# need to load both main and planet-specific ephemerides
# https://github.com/skyfielders/python-skyfield/issues/145
main_ephemeris = load('de421.bsp')
ephemeris = load('jup365.bsp')
ephemeris.segments.extend(main_ephemeris.segments)

earth = ephemeris['EARTH']
sun = ephemeris['SUN']
jup = ephemeris['JUPITER']
io = ephemeris['IO']

ts = load.timescale()
t = ts.utc(2022,9,27,16,4,0) # near opposition
t = ts.utc(2023,5,31,16,4,0) # some other time with a shadow

obs = earth.at(t)
target = jup
moon = io

obs_target = earth.at(t).observe(target)
t_target = t - obs_target.light_time
target_moon = target.at(t_target).observe(moon)
target_sun = target.at(t_target).observe(sun)

# want to project the shadow onto a plane
# calculate a reference frame that has a plane perpendicular to the shadow
shadow_vector = (target_sun-target_moon).position.au
R = calcR(shadow_vector)
print(R)

# rotation matrix will be orthogonal if built correctly
# i.e. R*R.T will be identity matrix
print('this should be an identity matrix:\n', mxm(R, R.T))

target_moon_pp = mxv(R, target_moon.position.au) # moon pos in projection plane coords
x, y, z = target_moon_pp # x is distance to moon, [y, z] is position of shadow
position_on_plane = [0, y, z]
shadow_icrf = mxv(R.T, position_on_plane) # convert back to ICRF

obs_shadow_vec = obs_target.position.au + shadow_icrf
obs_shadow = ICRF(obs_shadow_vec, t=t, center=target.target)

print('target: ', obs_target.radec(epoch='date')[0:2])
print('shadow:', obs_shadow.radec(epoch='date')[0:2])

I realized that what I'm doing is calculating solar eclipses on Jupiter. While researching that, I came across Besselian elements. It looks like eclipses on Earth are calculated the same way, just a much more rigorous version of it. Then, after calculating the shadow's position on the plane, they project it back up to the surface. I haven't gone through all of it yet, but it should be possible to adapt the technique for other planets. https://en.wikipedia.org/wiki/Besselian_elements http://ytliu.epizy.com/eclipse/Besselian_solar.html https://www.eclipsewise.com/solar/SEhelp/SEbeselm.html

[davidmikolas]

Orbits of the Galilean moons and occultations thereof will no doubt be interesting to folks wanting to calculate and/or watch the mutual occultations (where one moon occults another) which happen roughly every six years.

So as soon as one implements an approximate ephemeris based on orbital elements, one can be sure there will be complaints in a few years that mutual occultation predictions are failing.

Thus it would be best if one could also provide a way where spice kernels for the Jupiter system could be imported.

• https://astronomy.stackexchange.com/questions/23477/when-will-the-next-series-of-mutual-eclipses-of-jupiters-moons-begin