solareclipseproject/psestartendtime.py at main · i230010/solareclipseproject · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
"""
psestartendtime.py
-----------------
Compute start and end times for penumbral or umbral intersections
using cubic Besselian polynomials and numerical root-finding.

The zero crossings of the distance function correspond to the
first and last contact of the eclipse shadow with the Earth.
"""

from typing import List, Tuple
import math

from scipy.optimize import brentq


# ---------------------------------------------------------------------------
# Polynomial utilities


def poly(coeffs: List[float], t: float) -> float:
    """
    Evaluate a cubic polynomial:
        P(t) = c0 + c1*t + c2*t^2 + c3*t^3

    Parameters
    ----------
    coeffs : list of float
        Polynomial coefficients [c0, c1, c2, c3].
    t : float
        Input variable (time).

    Returns
    -------
    float
        Polynomial evaluated at t.
    """
    return coeffs[0] + coeffs[1] * t + coeffs[2] * t * t + coeffs[3] * t * t * t


# ---------------------------------------------------------------------------
# Eclipse distance function
# ---------------------------------------------------------------------------


def penumbra_distance(
    t: float,
    x_coeffs: List[float],
    y_coeffs: List[float],
    l_coeffs: List[float],
) -> float:
    """
    Compute the distance function whose zero defines shadow contact.

        D(t) = sqrt(x(t)^2 + y(t)^2) - (1 + L(t))

    A root of D(t) corresponds to the time when the penumbral or umbral
    boundary intersects the Earth's disk.

    Parameters
    ----------
    t : float
        Time parameter (same units as the Besselian polynomials).
    x_coeffs, y_coeffs, l_coeffs : sequence of float
        Cubic polynomial coefficients for X(t), Y(t), and L(t).

    Returns
    -------
    float
        Signed distance value.
    """
    x: float = poly(x_coeffs, t)
    y: float = poly(y_coeffs, t)
    radius: float = poly(l_coeffs, t)

    # hypot(x, y) computes sqrt(x^2 + y^2) in a numerically stable way
    return math.hypot(x, y) - (1.0 + radius)


# ---------------------------------------------------------------------------
# Start and end time solver
# ---------------------------------------------------------------------------


def startendtime(
    x_coeffs: List[float],
    y_coeffs: List[float],
    l_coeffs: List[float],
    t_start: float = -6.0,
    t_mid: float = 0.0,
    t_end: float = 6.0,
) -> Tuple[float, float]:
    """
    Solve for the start and end times of penumbral or umbral contact (TT).

    This function assumes the eclipse event is centered near t = 0 and
    that the distance function changes sign across the provided brackets.

    Parameters
    ----------
    x_coeffs, y_coeffs, l_coeffs : sequence of float
        Cubic polynomial coefficients for Besselian elements X(t), Y(t), L(t).
    t_start : float, optional
        Lower bound for the start time search (default: -6).
    t_mid : float, optional
        Midpoint separating start and end roots (default: 0).
    t_end : float, optional
        Upper bound for the end time search (default: 6).

    Returns
    -------
    tuple of float
        (start_time in TT, end_time in TT) in the same units as the input polynomials.

    Raises
    ------
    ValueError
        If the root is not bracketed within the provided intervals.
    """
    # Solve for first contact (ingress)
    start_time: float = brentq(
        penumbra_distance,
        t_start,
        t_mid,
        args=(x_coeffs, y_coeffs, l_coeffs),
    )

    # Solve for last contact (egress)
    end_time: float = brentq(
        penumbra_distance,
        t_mid,
        t_end,
        args=(x_coeffs, y_coeffs, l_coeffs),
    )

    return start_time, end_time