4

I'm trying to determine the latitude and longitude of say the Sun, the Moon and Mars. I need the result relative to the Earth's equator and the Prime Meridian in order to produce a result similar to this map.

I believe that's also what the author of this question wanted, however the answer there doesn't add up for me (comparing with values from the first link).

Expected result, obtained from the page linked to earlier:

On Thursday, 1 January 2015, 00:00:00 UTC the Sun is at its zenith at Latitude: 23° 02' South, Longitude: 179° 29' West

>>> import ephem; from math import degrees
>>> b = ephem.Sun(epoch='date'); b.compute('2015/1/1 00:00:00')
>>> print("{},{}".format(degrees(b.dec), degrees(b.ra)))
-23.040580418272267,281.12827017399906

So the latitude/declination seems about right, but no 180° wraparound will fix that right ascension, probably because it starts at the Vernal Equinox.

I have also unsuccessfully tried to use an observer at 0,0.

Can this be done using PyEphem, Skyfield or astropy? It seems odd that artificial satellites in PyEphem have the convenient sublat and sublong attributes, but it's so hard for celestial bodies.

4

1 回答 1

2

我终于弄明白了。有点。实际上，我只是将 libastro 的相关部分移植到 Python。请注意，此代码针对 Skyfield 的当前 git 版本（be6c7296）运行。

这里是（要点版本）：

#!/usr/bin/env python3

from datetime import datetime, timezone
from math import atan, atan2, degrees, floor, pi, radians, sin, sqrt

from skyfield.api import earth, JulianDate, now, sun


def earth_latlon(x, y, z, time):
    """
    For an object at the given XYZ coordinates relative to the center of
    the Earth at the given datetime, returns the latitude and longitude
    as it would appear on a world map.

    Units for XYZ don't matter.
    """
    julian_date = JulianDate(utc=time).tt
    # see https://en.wikipedia.org/wiki/Julian_date#Variants
    # libastro calls this "mjd", but the "Modified Julian Date" is
    # something entirely different
    dublin_julian_date = julian_date - 2415020

    # the following block closely mirrors libastro, so don't blame me
    # if you have no clue what the variables mean or what the magic
    # numbers are because I don't either
    sidereal_solar = 1.0027379093
    sid_day = floor(dublin_julian_date)
    t = (sid_day - 0.5) / 36525
    sid_reference = (6.6460656 + (2400.051262 * t) + (0.00002581 * (t**2))) / 24
    sid_reference -= floor(sid_reference)
    lon = 2 * pi * ((dublin_julian_date - sid_day) *
                    sidereal_solar + sid_reference) - atan2(y, x)
    lon = lon % (2 * pi)
    lon -= pi
    lat = atan(z / sqrt(x**2 + y**2))

    return degrees(lat), degrees(-lon)


if __name__ == '__main__':
    print("2015-01-01 00:00:00:")
    time = datetime(2015, 1, 1, tzinfo=timezone.utc)
    x, y, z = earth(JulianDate(utc=time)).observe(sun).apparent().position.au
    print(earth_latlon(x, y, z, time))

    print("now:")
    time = datetime.now(timezone.utc)
    x, y, z = earth(JulianDate(utc=time)).observe(sun).apparent().position.au
    print(earth_latlon(x, y, z, time))

输出：

2015-01-01 00:00:00:
(-23.05923949080624, -179.2173856294249)
now:
(-8.384551051991025, -47.12917634395421)

如您所见， 2015-01-01 00:00:00 的值与问题中的参考值相匹配。不完全是，但对我来说已经足够了。据我所知，我的价值观可能会更好。

由于我对 libastro 代码中使用的无证幻数一无所知，我无法为地球以外的物体进行这项工作。

@BrandonRhodes：如果您有兴趣在 Skyfield 中拥有此功能，请告诉我，然后我会尝试提出一个拉取请求。

于 2015-02-27T15:52:48.867 回答