Получение широты и долготы Солнца на карте мира с помощью PyEphem

Question

Я пытаюсь определить широту и долготу, скажем, Солнца, Луны и Марса. Мне нужен результат относительно экватора Земли и Первичного меридиана, чтобы получить результат, подобный этой карте.

Я полагаю, что это также то, что хотел автор этого вопроса, однако ответ там не складывается для меня (сравнение со значениями из первой ссылки).

Ожидаемый результат, полученный со страницы, на которую ссылались ранее:

В четверг, 1 января 2015 года, 00:00:00 UTC Солнце находится в зените на широте: 23 ° 02 'южной широты, долгота: 179 ° 29' западной долготы

>>> 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

Таким образом, широта / склонение кажутся примерно правильными, но никакое вращение на 180° не исправит правильное восхождение, вероятно, потому что оно начинается в день весеннего равноденствия.

Я также безуспешно пытался использовать наблюдателя в 0,0.

Можно ли это сделать с помощью PyEphem, Skyfield или astropy? Кажется странным, что искусственные спутники в PyEphem имеют удобные атрибуты sublat и sublong, но это очень трудно для небесных тел.

Answer 1

Я наконец-то понял. Вроде, как бы, что-то вроде. На самом деле я просто портировал соответствующие фрагменты libastro на Python. Обратите внимание, что этот код работает с текущей версией git Skyfield (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, тогда я постараюсь собрать запрос на извлечение.