OWAO 2025, Практически тур — Задача 2. At the World’s Edge

A Analysis of Eclipse Visibility with AstroPy

"""
Analysis of Lunar Eclipse Visibility
====================================

This script reads a catalog of 21st century lunar eclipses from
‘eclipses.xlsx‘, extracts the times of totality, and determines
which eclipses are visible from:
    - Cape Fligely
    - Point Xi on the equator

Visibility is defined as the Moon’s upper limb being above the
geometric horizon, corrected for refraction and the Moon’s
radius.
"""

import pandas as pd
import numpy as np
import astropy.units as u
from astropy.time import Time, TimeDelta
from astropy.coordinates import EarthLocation, AltAz, get_body


# ----------------------------------------------------------------------
# Load and preprocess the eclipse catalog
# ----------------------------------------------------------------------
def load_eclipses(path: str = "eclipses.xlsx") -> pd.DataFrame:
    """
    Load the lunar eclipse catalog from Excel and standardize columns.

    Parameters
    ----------
    path : str
        Path to the Excel file.

    Returns
    -------
    df : pandas.DataFrame
        DataFrame with standardized columns and auxiliary fields.
    """
    df = pd.read_excel(path, skiprows=1)

    df = df.rename(columns={
        "Unnamed: 0": "num",
        "U1": "U1", "U2": "U2", "U3": "U3", "U4": "U4",
        "Unnamed: 10": "maxt"
    })

    months = {m: i for i, m in enumerate(
        ["January","February","March","April","May","June",
         "July","August","September","October","November","December"], 1)}
    df["Month_num"] = df["Month"].map(months)

    df["date"] = pd.to_datetime(
        dict(year=df["Year"], month=df["Month_num"], day=df["Day"]),
        errors="coerce"
    )

    def frac_to_time(frac):
        if pd.isna(frac):
            return np.nan
        seconds = float(frac) * 24 * 3600
        h = int(seconds // 3600)
        m = int((seconds % 3600) // 60)
        s = int(seconds % 60)
        return f"{h:02d}:{m:02d}:{s:02d}"

    df["maxt_str"] = df["maxt"].apply(frac_to_time)

    return df


# ----------------------------------------------------------------------
# Convert tabulated times into astropy.Time arrays
# ----------------------------------------------------------------------
def build_times(df: pd.DataFrame):
    """
    Construct astropy.Time objects for total eclipses only.

    Returns
    -------
    df_tot : pandas.DataFrame
        Subset containing only total eclipses (U2 and U3 defined).
    t2, t3, tmax : astropy.Time
        Times of second contact, third contact, and maximum eclipse.
    """
    df_tot = df[df["U2"].notna() & df["U3"].notna()].copy()

    def to_iso(dates, times):
        iso_list = []
        for d, t in zip(dates, times):
            if pd.isna(t) or pd.isna(d):
                iso_list.append("NaT")
            else:
                iso_list.append(f"{d.strftime(’%Y-%m-%d’)}T{t}")
        return np.array(iso_list, dtype=str)

    iso_u2  = to_iso(df_tot["date"], df_tot["U2"])
    iso_u3  = to_iso(df_tot["date"], df_tot["U3"])
    iso_max = to_iso(df_tot["date"], df_tot["maxt_str"])

    t2   = Time(iso_u2,  format="isot", scale="utc", out_subfmt="date_hms")
    t3   = Time(iso_u3,  format="isot", scale="utc", out_subfmt="date_hms")
    tmax = Time(iso_max, format="isot", scale="utc", out_subfmt="date_hms")

    # adjust for midnight crossings
    one_day = TimeDelta(1*u.day)
    t3[tmax > t3] = t3[tmax > t3] + one_day
    t2[tmax < t2] = t2[tmax < t2] - one_day

    return df_tot, t2, t3, tmax


# ----------------------------------------------------------------------
# Visibility criterion
# ----------------------------------------------------------------------
def visible_mask_for_site(t2, t3, lat_deg, lon_deg=59 + 14/60, height_m=0,
                          n_samples=11, h_thresh_deg=-(35+16)/60.0):
    """
    Check eclipse visibility over the full totality interval [U2, U3]
    Samples each eclipse at ‘n_samples‘ evenly spaced times
    and returns True if the Moons center altitude exceeds ‘h_thresh_deg‘

    Parameters
    ----------
    t2, t3 : astropy.Time
        Start and end of totality.
    lat_deg, lon_deg : float
        Geographic coordinates of the site (degrees).
    height_m : float
        Elevation above sea level.
    n_samples : int, optional
        Number of samples per eclipse (inclusive of endpoints), default 11.
    h_thresh_deg : float, optional
        Altitude threshold (degrees) for the Moons center.

    Returns
    -------
    mask : ndarray of bool
        True if any sample during totality is above the threshold
        (i.e., visible).
    """
    loc = EarthLocation(
        lat=lat_deg*u.deg,
        lon=lon_deg*u.deg,
        height=height_m*u.m
    )

    frac = np.linspace(0.0, 1.0, n_samples)
    dt = t3 - t2
    t_grid = t2[:, None] + frac[None, :] * dt[:, None]

    alt = get_body("moon", t_grid.reshape(-1), location=loc)\
        .transform_to(AltAz(obstime=t_grid.reshape(-1), location=loc))\
        .alt.deg.reshape(len(t2), n_samples)

    return (alt > h_thresh_deg).any(axis=1)


# ----------------------------------------------------------------------
# Main pipeline
# ----------------------------------------------------------------------
def analyze_visibility(path="eclipses.xlsx"):
    """
    Analyze visibility of total lunar eclipses from Cape Fligely
    and point Xi.
    """
    df = load_eclipses(path)
    df_tot, t2, t3, tmax = build_times(df)

    events = df_tot["num"].to_numpy()

    # two sites
    vis_XI = visible_mask_for_site(t2, t3, lat_deg=0.0)
    vis_FL = visible_mask_for_site(t2, t3, lat_deg=81 + 51/60.0)

    xi_list  = events[vis_XI].tolist()
    fl_list  = events[vis_FL].tolist()
    both     = sorted(set(xi_list) & set(fl_list))

    print("Total eclipses:", len(events))
    print("Visible from point Xi:", len(xi_list))
    print("Visible from Cape Fligely:", len(fl_list))
    print("Visible from both:", len(both))

    return xi_list, fl_list, both


# ----------------------------------------------------------------------
# Entry point
# ----------------------------------------------------------------------
if __name__ == "__main__":
    xi_list, fl_list, both = analyze_visibility("eclipses.xlsx")