Определение лунного затмения в небесном поле

Мне дан список дат в формате UTC, все часы переведены на 00:00.

Я хотел бы определить, произошло ли (лунное) затмение в данный день (то есть за последние 24 часа)

Учитывая фрагмент Python

from sykfield.api import load
eph = load('de421.bsp')
def eclipticangle(t):

    moon, earth = eph['moon'], eph['earth']
    e = earth.at(t)
    x, y, _ = e.observe(moon).apparent().ecliptic_latlon()

    return x.degrees

Я предполагаю, что можно определить, произошло ли затмение в течение 24 часов времени t на

  1. Проверка того, что первый угол достаточно близок к 180 (легко)
  2. Проверка, достаточно ли вторая степень близка к 0 (не так уж и сложно?)

Теперь, что касается ответа в комментарии, не так уж и тривиально решить вторую проблему, просто проверив, близок ли угол к 0.

Поэтому у меня вопрос

Может ли кто-нибудь предоставить функцию для определения того, произошло ли лунное затмение в данный день t?

Изменить. Этот вопрос был отредактирован с учетом отзывов Брэндона Роудса, оставленных в комментариях ниже.


person Jernej    schedule 03.11.2020    source источник
comment
Будет ли первая проблема - которая потребует от вас знания скорости, с которой Луна движется по небу, которая меняется - исчезнет, ​​если вы спросите положение Луны в «t минус один день» и проверите, было ли это положение на другом? сторона 180 °?   -  person Brandon Rhodes    schedule 03.11.2020
comment
@BrandonRhodes Спасибо, это отлично работает. Можно ли применить аналогичную логику и к затмениям? Монотонно ли ведет себя угол между двумя плоскостями? Если нет, то какой разумный способ определять затмения?   -  person Jernej    schedule 03.11.2020
comment
- Увы, я никогда не писал никаких программ для поиска затмений (поэтому я добавил просто комментарий и не пытался ответить на ваш вопрос). Я знаю, что он включает в себя трехмерную форму тени Земли в космосе, но, поскольку тень может быть аппроксимирована парой конусов (тени и полутени) с точки зрения Луны, диаграммы всегда, кажется, аппроксимируют тень как пару кругов, через которые проходит Луна - плоское сечение этих двух конусов. Полагаю, они измеряют положение Луны относительно тех кругов тени на расстоянии?   -  person Brandon Rhodes    schedule 05.11.2020
comment
@BrandonRhodes Спасибо за общую идею. Это говорит мне, что у меня недостаточно опыта, чтобы самому его кодировать. Поэтому я просто рефакторирую этот вопрос и увеличу небольшую награду.   -  person Jernej    schedule 05.11.2020
comment
Поскольку дни идут в обратном порядке, я спрошу здесь в комментарии, приближается ли мой ответ к условиям вашей награды, или он оставляет какие-то важные вопросы по-прежнему открытыми. Спасибо!   -  person Brandon Rhodes    schedule 10.11.2020
comment
@BrandonRhodes Да, пока это выглядит именно так, как я хотел. Я приму ответ до истечения срока действия награды, так как перед этим хочу провести один обширный тест.   -  person Jernej    schedule 10.11.2020
comment
Мне будет интересно услышать, как выйдет тест!   -  person Brandon Rhodes    schedule 13.11.2020


Ответы (1)


Я только что просмотрел раздел 11.2.3 Пояснительного дополнения к астрономическому альманаху и попытался превратить его в код Skyfield Python. Вот что я придумал:

import numpy as np

from skyfield.api import load
from skyfield.constants import ERAD
from skyfield.functions import angle_between, length_of
from skyfield.searchlib import find_maxima

eph = load('de421.bsp')
earth = eph['earth']
moon = eph['moon']
sun = eph['sun']

def f(t):
    e = earth.at(t).position.au
    s = sun.at(t).position.au
    m = moon.at(t).position.au
    return angle_between(s - e, m - e)

f.step_days = 5.0

ts = load.timescale()
start_time = ts.utc(2019, 1, 1)
end_time = ts.utc(2020, 1, 1)

t, y = find_maxima(start_time, end_time, f)

e = earth.at(t).position.m
m = moon.at(t).position.m
s = sun.at(t).position.m

solar_radius_m = 696340e3
moon_radius_m = 1.7371e6

pi_m = np.arcsin(ERAD / length_of(m - e))
pi_s = np.arcsin(ERAD / length_of(s - e))
s_s = np.arcsin(solar_radius_m / length_of(s - e))

pi_1 = 0.998340 * pi_m

sigma = angle_between(s - e, e - m)
s_m = np.arcsin(moon_radius_m / length_of(e - m))

penumbral = sigma < 1.02 * (pi_1 + pi_s + s_s) + s_m
partial = sigma < 1.02 * (pi_1 + pi_s - s_s) + s_m
total = sigma < 1.02 * (pi_1 + pi_s - s_s) - s_m

mask = penumbral | partial | total

t = t[mask]
penumbral = penumbral[mask]
partial = partial[mask]
total = total[mask]

print(t.utc_strftime())
print(0 + penumbral + partial + total)

Он создает вектор времени, когда происходили лунные затмения, а затем оценку того, насколько полное затмение:

['2019-01-21 05:12:51 UTC', '2019-07-16 21:31:27 UTC']
[3 2]

Время его затмений находится в пределах 3 секунд от времени, указанного в огромной таблице лунных эфемерид НАСА:

https://eclipse.gsfc.nasa.gov/5MCLE/5MKLEcatalog.txt

person Brandon Rhodes    schedule 07.11.2020