Add routine and test for the solar eclipse calculation.

skyfield/eclipselib.py

@@ -3,10 +3,10 @@
 
 from __future__ import division
 
-from numpy import arcsin, byte
+from numpy import arcsin, byte, arctan2, sin, pi, arccos, sqrt, dot, array, empty
 from .constants import AU_KM, C_AUDAY, ERAD
 from .functions import angle_between, length_of
-from .searchlib import find_maxima
+from .searchlib import find_maxima, find_minima
 from .relativity import add_aberration
 
 LUNAR_ECLIPSES = [
@@ -15,6 +15,35 @@
     'Total',
 ]
 
+SOLAR_ECLIPSES = [
+    'Partial',
+    'Total/Hybrid',
+    'Annular'
+]
+
+
+def _compute_angle(t, earth_barycenter, earth, moon, sun):
+    jd, fr = t.whole, t.tdb_fraction
+    b, velocity = earth_barycenter.compute_and_differentiate(jd, fr)
+    e = earth.compute(jd, fr)
+    m = moon.compute(jd, fr)
+    s = sun.compute(jd, fr)
+
+    earth_to_sun = s - b - e
+    earth_to_moon = m - e
+
+    # The aberration routine requires specific units.  (We can leave
+    # the `earth_to_moon` vector unconverted because we only need
+    # its direction.)  We approximate the Earth’s velocity as being
+    # that of the Earth-Moon barycenter.
+
+    earth_to_sun /= AU_KM
+    velocity /= AU_KM
+    light_travel_time = length_of(earth_to_sun) / C_AUDAY
+    add_aberration(earth_to_sun, velocity, light_travel_time)
+
+    return angle_between(earth_to_sun, earth_to_moon)
+
 def lunar_eclipses(start_time, end_time, eph):
     """Return the lunar eclipses between ``start_time`` and ``end_time``.
 
@@ -46,30 +75,10 @@ def lunar_eclipses(start_time, end_time, eph):
     earth = sdict[3,399]
     moon = sdict[3,301]
 
-    def f(t):
-        jd, fr = t.whole, t.tdb_fraction
-        b, velocity = earth_barycenter.compute_and_differentiate(jd, fr)
-        e = earth.compute(jd, fr)
-        m = moon.compute(jd, fr)
-        s = sun.compute(jd, fr)
-
-        earth_to_sun = s - b - e
-        earth_to_moon = m - e
-
-        # The aberration routine requires specific units.  (We can leave
-        # the `earth_to_moon` vector unconverted because we only need
-        # its direction.)  We approximate the Earth’s velocity as being
-        # that of the Earth-Moon barycenter.
-
-        earth_to_sun /= AU_KM
-        velocity /= AU_KM
-        light_travel_time = length_of(earth_to_sun) / C_AUDAY
-        add_aberration(earth_to_sun, velocity, light_travel_time)
-
-        return angle_between(earth_to_sun, earth_to_moon)
-
+    f = lambda x: _compute_angle(x, earth_barycenter, earth, moon, sun)
     f.step_days = 5.0
-    t, y = find_maxima(start_time, end_time, f, num=4)
+
+    t, y = find_maxima(start_time, end_time, f , num=4)
 
     jd, fr = t.whole, t.tdb_fraction
     b = earth_barycenter.compute(jd, fr)
@@ -132,3 +141,173 @@ def f(t):
     }
 
     return t, code, details
+
+def _cross_area(d, r, R):
+    """
+    Computes common area of the circles with radius ``r`` and ``R``
+    that are separated by distance ``d``.
+    """
+
+    return (
+        r * r * arccos((d * d + r * r - R * R) / (2 * d * r))
+        + R * R * arccos((d * d + R * R - r * r) / (2 * d * R))
+        - 1 / 2 * sqrt((-d + r + R) * (d + r - R) * (d - r + R) * (d + r + R))
+    )
+
+
+def _line_sphere_intersection(unit_vec, sphere_center, sphere_radius):
+    """
+    Computes closest point of line intersecting the sphere.
+    The line starts at origin (0, 0, 0) and has unit vector ``unit_vec``.
+    If line does intersect the sphere, returns closest point of sphere
+    surface to the line.
+    """
+
+    center_sq = dot(sphere_center, sphere_center)
+    prod1 = dot(unit_vec, sphere_center)
+    delta = prod1 * prod1 - center_sq + sphere_radius * sphere_radius
+    if delta < 0:
+        s_to_point = unit_vec * prod1
+        unit_e_to_point = (s_to_point - sphere_center) / length_of(
+            s_to_point - sphere_center
+        )
+        return sphere_center + unit_e_to_point * sphere_radius
+    else:
+        d = prod1 - sqrt(delta)
+        return unit_vec * d
+
+
+def _magnitude_obscurity(e_to_s, e_to_m, s_radius, m_radius, e_radius):
+    """
+    Calculates magnitude and obscurity of the eclipse.
+    """
+
+    s_to_m = e_to_m - e_to_s
+    s_to_ecl = _line_sphere_intersection(s_to_m / length_of(s_to_m),
+                                         -e_to_s, e_radius)
+    m_to_ecl = s_to_ecl - s_to_m
+    scale = length_of(s_to_ecl) / length_of(m_to_ecl)
+    s_radius_res = s_radius / scale
+    s_to_ecl_res = s_to_ecl / scale
+    dist = length_of(m_to_ecl - s_to_ecl_res)
+    if dist >= s_radius_res + m_radius:
+        return 0.0, 0.0
+    mag = m_radius / s_radius_res
+    if m_radius + dist <= s_radius_res:
+        obs = m_radius * m_radius / (s_radius_res * s_radius_res)
+    elif s_radius_res + dist <= m_radius:
+        obs = 1.0
+    else:
+        mag = (s_radius_res + m_radius - dist) / (2 * s_radius_res)
+        cross_a = _cross_area(dist, s_radius_res, m_radius)
+        obs = cross_a / (pi * s_radius_res * s_radius_res)
+    return mag, obs
+
+
+def solar_eclipses(start_time, end_time, eph):
+    """Return the solar eclipses between ``start_time`` and ``end_time``.
+
+    Returns a 4-item tuple:
+
+    * A :class:`~skyfield.timelib.Time` giving the dates of each eclipse.
+    * An integer array of codes identifying what type eclipse is.
+    * A float indicating the magnitude of the eclipse.
+    * A float indicating the obscurity of the eclipse.
+
+    Comments regarding approximations for the method lunar_eclipses()
+    are applicable for this routine as well.
+    """
+
+    sdict = dict(((s.center, s.target), s.spk_segment) for s in eph.segments)
+    sun = sdict[0,10]
+    earth_barycenter = sdict[0,3]
+    earth = sdict[3,399]
+    moon = sdict[3,301]
+
+    f = lambda x: _compute_angle(x, earth_barycenter, earth, moon, sun)
+    f.step_days = 5.0
+
+    t, y = find_minima(start_time, end_time, f, num=4)
+
+    jd, fr = t.whole, t.tdb_fraction
+    b, velocity = earth_barycenter.compute_and_differentiate(jd, fr)
+    e = earth.compute(jd, fr)
+    m = moon.compute(jd, fr)
+    s = sun.compute(jd, fr)
+
+    earth_to_sun = s - b - e
+    moon_to_earth = e - m
+    moon_to_sun = s - b - m
+
+    earth_to_sun_AU = earth_to_sun / AU_KM
+    add_aberration(earth_to_sun_AU, velocity / AU_KM,
+                   length_of(earth_to_sun_AU) / C_AUDAY)
+    earth_to_sun = earth_to_sun_AU * AU_KM
+
+    solar_radius_km = 696340.0
+    moon_radius_km = 1737.1
+    earth_radius_km = ERAD / 1000
+
+    angle_penumbra = arctan2(moon_radius_km + solar_radius_km, length_of(moon_to_sun))
+    origin_pen_from_sun = (length_of(moon_to_sun) * solar_radius_km /
+                           (moon_radius_km + solar_radius_km))
+    origin_pen = s - origin_pen_from_sun * moon_to_sun / length_of(moon_to_sun)
+    pen_to_earth = b + e - origin_pen
+    angle_earth_from_pen = angle_between(-moon_to_sun, pen_to_earth)
+
+    earth_from_pen_cone = (sin(angle_earth_from_pen - angle_penumbra) *
+                           length_of(pen_to_earth))
+    is_eclipse = ((angle_earth_from_pen < angle_penumbra) +
+                  (earth_from_pen_cone <= earth_radius_km))
+
+    origin_umbra_from_sun = (solar_radius_km * length_of(moon_to_sun) /
+                             (solar_radius_km - moon_radius_km))
+    angle_umbra = arctan2(solar_radius_km, origin_umbra_from_sun)
+    origin_umbra = s - origin_umbra_from_sun * moon_to_sun / length_of(moon_to_sun)
+    umbra_to_earth = b + e - origin_umbra
+
+    angle_earth_from_umbra = angle_between(moon_to_sun, umbra_to_earth)
+
+    is_angle_inside_umbra = angle_umbra - angle_earth_from_umbra > 0
+    is_umbra_inside_earth_radius = length_of(umbra_to_earth) <= earth_radius_km
+
+    angle_antumbra = pi - angle_umbra
+    earth_from_antumbra_cone = sin(angle_antumbra - angle_earth_from_umbra) * length_of(
+        umbra_to_earth
+    )
+    is_eclipse_annular = (angle_earth_from_umbra - angle_antumbra > 0) + (
+            earth_from_antumbra_cone <= earth_radius_km
+    )
+    is_eclipse_annular *= [not x for x in is_umbra_inside_earth_radius]
+
+    earth_from_umbra_cone = sin(angle_earth_from_umbra - angle_umbra) * length_of(
+        umbra_to_earth
+    )
+    is_earth_radius_inside_umbra = earth_radius_km >= earth_from_umbra_cone
+
+    is_eclipse_total = is_eclipse * (
+            is_angle_inside_umbra
+            + ((is_earth_radius_inside_umbra) *
+               (angle_earth_from_umbra - angle_umbra < pi / 2))
+            + is_umbra_inside_earth_radius
+    )
+
+    times = t[is_eclipse]
+    code = is_eclipse_total.astype(byte)
+    code += 2 * (is_eclipse_annular.astype(byte) -
+                 (is_eclipse_annular * is_eclipse_total).astype(byte))
+    code = code[is_eclipse]
+
+    mag = empty(sum(is_eclipse))
+    obs = empty(sum(is_eclipse))
+
+    for i in range(len(times)):
+        mag[i], obs[i] = _magnitude_obscurity(
+            array([array(x[is_eclipse][i]) for x in earth_to_sun]),
+            -array([array(x[is_eclipse][i]) for x in moon_to_earth]),
+            solar_radius_km,
+            moon_radius_km,
+            earth_radius_km,
+        )
+
+    return times, code, mag, obs

skyfield/tests/test_eclipselib.py

@@ -14,3 +14,19 @@ def test_lunar_eclipses():
     assert len(t) == len(y) == 2
     for name, item in details.items():
         assert len(item) == len(t)
+
+def test_solar_eclipses():
+    ts = load.timescale()
+    eph = load('de421.bsp')
+
+    t0 = ts.utc(2019, 1, 1)
+    t1 = ts.utc(2020, 1, 1)
+    times, codes, mag, obs = eclipselib.solar_eclipses(t0, t1, eph)
+
+    result = []
+    for t, c, m, o in zip(times, codes, mag, obs):
+        result += [f"{t.tt_strftime()} {eclipselib.SOLAR_ECLIPSES[c]} {m:.4} {o:.2}"]
+
+    assert ','.join(result) == ("2019-01-06 01:42:40 TT Partial 0.7192 0.63,"
+                                "2019-07-02 19:24:04 TT Total/Hybrid 1.046 1.0,"
+                                "2019-12-26 05:18:52 TT Annular 0.9699 0.94")