8 апреля 2024 года автор статьи, основатель и СЕО компании Modal Labs, Эрик Бернхардссон планировал посмотреть свое первое полное солнечное затмение. За день до этого ему пришла в голову идея — что, если попробовать рассчитать периодичность этого явления в прошлом и будущем, используя Python? Несмотря на незначительные сложности с системой координат, автору удалось создать работоспособное решение всего за несколько часов.

Под катом читайте, как с помощью ~100 строк кода удалось вычислить и проследить путь каждого солнечного затмения в период с 2020 по 2030 год.

*Обращаем ваше внимание, что позиция автора может не всегда совпадать с мнением МойОфис.

Когда я собрался посмотреть солнечное затмение, мне стало любопытно, насколько сложно вычислить его периодичность с помощью Python. Я не хотел углубляться в астрономию для расчетов, поэтому решил использовать замечательную экосистему Python для всего. Оказалось, что в пакете Astropy есть около 80 %  того, что мне нужно. В частности, пакет позволяет довольно просто вычислять положение Солнца и Луны на небе.

После нескольких минут гугления у меня появилось код, вычисляющий перекрытие Солнца Луной для определенной точки Земли:

from astropy.coordinates import EarthLocation, get_body
from astropy.time import Time
from astropy.units import deg, m

def sun_moon_separation(lat: float, lon: float, t: float) -> float:
    loc = EarthLocation(lat=lat * deg, lon=lon * deg, height=0 * m)
    time = Time(t, format="unix")
    moon = get_body("moon", time, loc)
    sun = get_body("sun", time, loc)

    sep = moon.separation(sun)
    return sep.deg

Для этого код берет пару широта/долгота, а также временную метку unix и рассчитывает угловое расстояние между Солнцем и Луной. В основном это просто означает расстояние между центрами объектов на небе, видимых с Земли. Если угловое расстояние очень близко к нулю, получается солнечное затмение.

Но! Моей целью были не вычисления для заданной координаты. Я хотел рассчитать местоположение полного затмения по временной метке (если она есть).

В идеале у нас должны быть 3D-координаты Земли, Солнца и Луны. Затем нам нужно спроецировать линию между Солнцем и Луной, чтобы посмотреть, проходит ли эта линия через Землю, и если да, то необходимо найти широту и долготу этого пересечения. Возможно, это «правильный» способ, и если бы у меня было время, я бы обратился к своим подзабытым знаниям по геометрии и сделал именно так.

Однако времени у меня не было. Затмение уже завтра, и я просто хотел рассчитать координаты наименее заковыристым способом. У нас уже есть инструмент, вычисляющий близкие параметры, но нам нужно немного изменить ситуацию. Мы сделаем это с помощью штуки, к которой я в подобных случаях всегда прибегаю — методу прямого поиска (оптимизация «черного ящика»).

Решение для координат с помощью метода прямого поиска

У нас есть функция, которая принимает три значения (временная метка, широта и долгота) и выдает угловое расстояние между Солнцем и Луной. Но давайте вместо этого попробуем решить смежную задачу: по заданной временной метке найдем широту и долготу, которые минимизируют расстояние между Солнцем и Луной.

Если минимальное расстояние практически равно нулю, это означает, что мы нашли солнечное затмение. В этом случае координата, минимизирующая функцию — центр лунной тени на Земле.

Минимизировать такую произвольную функцию довольно просто. Для этого отлично подойдет пакет scipy.optimize,  где есть куча проверенных функций, которые, вероятно, были реализованы на Фортране 77 – при желании, можете проверить это сами. У нас даже нет градиента функции, но это не страшно — Нелдер-Мид вам в помощь.

Самое приятное, что мы можем рассматривать эту функцию как абсолютный «черный ящик» и оптимизировать ее извне. Правда, это требует определенных вычислительных затрат, но лично я не стал бы из-за этого переживать.

Код для использования scipy.optimize.minimize, чтобы найти местоположение затмения в итоге будет таким:

def find_eclipse_location(dt: datetime) -> tuple[float, float] | None
    """Возвращает координаты полного затмения или `None`"""
    t = datetime.timestamp(dt)
    fun = lambda x: sun_moon_separation(x[0], x[1], t)

    ret = minimize(fun, bounds=[(-90, 90), (-180, 180)], x0=(0, 0))
    return ret.x if ret.fun < 1e-3 else None

По сути, мы привязываем время к sun_moon_separation и строим новую функцию с двумя переменными: широтой и долготой. А затем перебираем эту функцию (с заданными ограничениями), чтобы найти минимум.

Все почти работает! Но часть проблемы заключалась в том, что я потратил два часа впустую из-за глупой ошибки со знаками в широтах и долготах. Но даже после исправления я получил странные искаженные координаты.

Я думаю, что это произошло из-за ложных минимумов, так как считаю, что антипод одного решения – это другое решение. Очевидно, мы должны отбрасывать решения, когда Солнца не видно. Две простые модификации — и все заработает:

1. Если Солнце или Луна ниже горизонта — вернем какое-нибудь большое число.

2. Вместо того, чтобы использовать (0, 0) в качестве исходной точки, выполните простой поиск по сетке в нескольких точках Земли и выберите ту, где расстояние между Солнцем и Луной будет наименьшим. Затем используйте эту точку в качестве исходной для оптимизации.

Мой финальный код для sun_moon_separation и find_eclipse_location в итоге получился лишь немного сложнее того, что я описал выше. С помощью этих модификаций мы получили функцию, которая надежно берет любую временную метку и определяет широту/долготу солнечного затмения (если оно есть).

Поиск всех затмений

Итак, давайте найдем больше затмений! В частности, выясним путь каждого затмения в промежутке между 2020 и 2030 годами. Для этого нам потребуется перебрать множество временных меток.

Увы, функция find_eclipse_location работает довольно медленно.

Как мы поступим? Добавим еще пару приемов:

1. Проводим грубый поиск по всему десятилетию, проверяя только каждый час. Если мы найдем затмение, выполняем более детальный поиск и прокладываем путь минута за минутой.

2. Распараллеливаем!!!

Я использовал фреймворк Modal, чьи инструменты позволяют взять код на Python и запустить его в облаке. Платформа отлично подходит для масштабирования функций, требующих больших вычислений.

Теперь мы можем найти все затмения 2020-30 годов. Добавляем простой декоратор к find_eclipse_location и затем выполняем процедуру маппинга. В итоге код выглядит так:

def run():
    dt_a = datetime(2020, 1, 1, 0, 0, 0, tzinfo=timezone.utc)
    dt_b = datetime(2030, 1, 1, 0, 0, 0, tzinfo=timezone.utc)

    # Compute evenly spaced datetimes
    dt = dt_a
    dts = []
    while dt < dt_b:
        dts.append(dt)
        dt = dt + timedelta(seconds=3600)

    # Map over it using Modal!!!
    for tup in find_eclipse_location.map(dts):
        if tup is not None:
            print("Found eclipse at", tup

Построение графика

Я опускаю некоторые детали в самом коде, но потерпите. Теперь, когда у нас есть все пути, можем начертить их. Я использовал Basemap и довольно быстро получил такое:

from matplotlib import pyplot
from mpl_toolkits.basemap import Basemap

def plot_path(dts: list[datetime], lats: list[float], lons: list[float]):
    # Set up a world map                                                                                                                                                                              
    pyplot.figure(figsize=(6, 6))
    lat_0, lon_0 = lats[len(lats) // 2], lons[len(lons) // 2]
    bm = Basemap(projection="ortho", lat_0=lat_0, lon_0=lon_0)
    bm.drawmapboundary(fill_color="navy")
    bm.fillcontinents(color="forestgreen", lake_color="blue")
    bm.drawcoastlines()

    # Plot eclipse path
    x, y = bm(lons, lats)
    bm.plot(x, y, color="red")

Я добавил еще несколько элементов в свой финальный скрипт, включая местное время, используя timezonefinder для поиска местных часовых поясов по парам широта/долгота.

Вот как будет выглядеть затмение 8 апреля 2024 года, если мы нанесем данные на карту с помощью нашего скрипта:

Великолепно!

На самом деле, за это вряд ли получишь награду в номинации «Лучший дизайнер», но это выглядит довольно прилично как отправная точка — смысл здесь не в красоте, а в поиске затмений в ~100 строках Python.

Что скрипт и делает! Фактически, он находит все затмения в период 2020-2030 гг.:

2020-06-21 над Африкой, Ближним Востоком и Азией

2020-12-14 над крошечным кусочком Южной Америки

2021-06-10 над северной Канадой и Гренландией

2021-12-04 над Антарктидой

2023-04-20 над Австралией и Папуа-Новой Гвинеей

2023-10-14 над США, Центральной Америкой и Южной Америкой

2024-04-08 над Мексикой, США и Канадой

2024-10-02 над небольшим кусочком Южной Америки (опять?)

2026-02-17 над крошечным кусочком Антарктиды (кто-нибудь увидит?)

2026-08-12 над Гренландией и Испанией

2027-02-06 над крошечным кусочком Южной Америки (и в третий раз??)

2027-08-02 над Северной Африкой и Ближним Востоком

2028-01-26 над Южной Америкой и Испанией

2028-07-22 над Австралией и Новой Зеландией

Мой список идентичен тем, что я нашел в сети, и это весьма обнадеживает.

Общее время расчетов – несколько минут.

Конечно, это грубый подход. Я уверен, что у NASA есть версия на C++, которая работает в 1000 раз быстрее. Однако с точки зрения продуктивности разработчиков это — победа, даже если не принимать во внимание тот факт, что мы также строили карты!

Заметки

• Счастливчики на юге Южной Америки ловят три затмения за десятилетие.

• Если хотите проверить код — он здесь.

• Я был вдохновлен этим постом, где для похожих на мои расчеты, использована система Mathematica.

• В качестве исходной точки в моем коде использован код из Stackoverflow.

• Я проигнорировал разницу между кольцевыми и полными затмениями в своем коде. Думаю, это несложно исправить при необходимости.

• Я также не вычислял ширину пути полного затмения, то есть ширину солнечной тени на Земле. Только путь центра этой тени.