skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.4k stars 211 forks source link

How to set up observations from bodies other than Earth? #709

Open AJMansfield opened 2 years ago

AJMansfield commented 2 years ago

I'm interested in plotting some night sky animations of the skies as observed from points on the surface of Mars and from some of the Jovian moons and minor planets, however I'm having trouble figuring out how to set up geographic positions on bodies other than Earth.

Setting up an reference ellipsoid for an alternate body isn't a problem; one can just construct skyfield.toposlib.Geoid with the appropriate parameters. E.g. for Mars:

# from https://ode.rsl.wustl.edu/mars/pagehelp/Content/Frequently_Asked_Questions/Coordinate_System.htm
# equatorial radius = 3396.19 km
# polar radius = 3376.2 km
iau2000_mars = Geoid("IAU2000 Mars", 3396190.0, 169.89)

Setting up the planet's rotation rate is a little less clear, but it seems like I'd just need to copy skyfield.toposlib.GeographicPosition, and replace the parameters in lst_hours_at with appropriate parameters for another body.

It's not all that obvious what the different constants in that function represent though, making it difficult to search for the equivalent parameters for other bodies of interest. But at least as a rough first approximation for Mars, it seems like I could just scale the output by the ratio of the mars and earth sidereal day lengths, and get at least the right average rotation rate:

class MarsGeographicPosition(GeographicPosition):
    ...
    def lst_hours_at(self, t):
        mars_earth_ratio = 1.0287656 # https://www.wolframalpha.com/input?i=mars+rotation+period+%2F+earth+rotation+period
        sprime = -47.0e-6 * (t.whole - T0 + t.tdb_fraction) / 36525.0
        return (t.gast + self.longitude._hours + sprime / 54000.0) * mars_earth_ratio % 24.0
    ...

I'd definitely like to set it up more accurately than that though.

And it's not clear at all to me how to set up obliquity parameters. Mars at least, seems to have this included in its de421.bsp the same as Earth, but I can't seem to find a way to set it up manually for bodies I'm only defining with e.g. a Kepler orbit.

To summarize:

How do I set up the parameters in GeographicPosition.lst_hours_at for bodies other than Earth? How do I set up obliquity for e.g. Kepler Orbit bodies?

brandon-rhodes commented 2 years ago

Good question. Skyfield at the moment provides only very simple support for non-Earth bodies, though returning a simple angle from lst_hours_at() which — as you say — increases at roughly the right rate for Mars might get you started. You can ignore the adjustments made for Earth's position by those little obliquity values; they appear because we know the Earth's position to such very high precision.

The official way that folks at the JPL deal with planetary reference frames is with a file like pck00008.tpc which is mentioned here:

https://rhodesmill.org/skyfield/planetary.html

All I've implemented so far is enough support for Moon locations, but the pck00008.tpc file has descriptions of several other bodies too. For example, try searching it for the word "Jupiter":

https://naif.jpl.nasa.gov/pub/naif/JUNO/kernels/pck/pck00008.tpc

You'll soon find this section of constants:

        BODY599_POLE_RA        = (   268.05      -0.009       0. )
        BODY599_POLE_DEC       = (    64.49       0.003       0. )
        BODY599_PM             = (   284.95     870.5366420   0. )
        BODY599_LONG_AXIS      = (     0.                        )

        BODY5_NUT_PREC_ANGLES  = (    73.32   91472.9
                                      24.62   45137.2
                                     283.90    4850.7
                                     355.80    1191.3
                                     119.90     262.1
                                     229.80      64.3
                                     352.35    2382.6
                                     113.35    6070.0   
                                     146.64  182945.8
                                      49.24   90274.4  )

You can see the data on where its pole points; on how fast its prime meridian (PM) spins; and even some nutation and precession angles if you want to correctly model how its pole wanders as it spins.

It feels like these should be enough to provide a reference frame, but since all I've implemented so far is the Moon, which has its frames defined explicitly in a file moon_080317.tf where they have names like MOON_ME_DE421. But what if you just want to take the Jupiter information from pck00008.tpc and treat it as a standalone reference frame? I'll go re-read Skyfield's planetarylib.py and re-read the SPICE documentation (which describes those files) and see how far Skyfield might or might not be from being able to understand something like Jupiter's orientation.