Partial vs Total eclipse

A total eclipse is what someone located in the umbra sees, the moon blocks all the incoming light of the sun. A partial eclipse is what someone in the penumbra sees, the moon only blocks some the sunlight.

Is there an eclipse right now ?

Approximations

• The sun, the earth and the moon are considered spherical.
• The other celestial bodies are ignored.

The interesting lines

Without further loss of precision, we can consider the sun, the moon and the earth as disks, each sitting on one of the 3 planes normal to the sun-earth axis and containing a body’s center.

To detect whether a total or partial eclipse is occurring, all we have to do is to find out if either the penumbra or the umbra intersects the earth’s disk. In concrete terms, we need to compute the bit of umbra and penumbra that is the closest to the sun-center-to-earth-center axis and see whether they fall on the earth’s disk or not.

The (normalized) vector moon-center to projection of the earth’s center on the moon plane, MC_EP gives the direction on the sun-facing planes in which to seek the points that define the lines leading to the interesting shadow bits, as shown in the following schemas.

In a given position, two lines (on for total, and one for partial) are enough determine if an eclipse is occurring. Each of these lines is defined by two points: one on the sun’s perimeter and another on the moon’s perimeter as demonstrated in the first schema.

If they intersect the earth’s disk then there’s an eclipse.

Note that in some circumstances this is not exactly correct. In such case the rules are slightly different and there is no total eclipse.

Computing these points

For different reasons, we’ll need to compute 3 intersections between a line and a plane.

1. To compute MC_EP, between
   • the sun-earth axis and
   • the moon’s plane
2. To find if there’s a total eclipse, between
   • the line $(rt)$, originating from the perimeter of the sun’s disk in the direction of MC_EP passing by the point of the moon’s perimeter that is the closest to the earth’s center (i.e. in the direction of MC_EP) and
   • the earth’s plane
3. To find if there’s a partial eclipse, between
   • the line $(rp)$, originating from the perimeter of the sun’s disk in the opposite direction of MC_EP passing by the point of the moon’s perimeter that is the closest to the earth’s center and
   • the earth’s plane

Let $A$, $B$ be two points, the intersection $I$ of the line directed by $\vec{AB}$ with the plane of normal $\vec{n}$ passing by $O$ is given by:

In python:

def intersectionPoint(O, n, A, B):
    """retruns, if it exists, the intersection point between:
       -the plane of normal vector n containing O
       -the line (AB)"""
    AB = B - A
    if abs(np.dot(AB, n)) < espilon:
        raise ValueError("No intersection")

    d = -(np.dot(A - O, n) / np.dot(AB, n))

    X = A + d * AB

    return X

That’s enough to create an isEclipse function that can tell for a given position of the bodies whether an eclipse is occurring or not.

def isEclipse(S, E, M, r_s, r_e, r_m):
    """returns a tuple of boolean, (partial, total)
       S, E, M: np array [x,y,z] of the position of the Sun, the Earth and the Moon
       r_s, r_e, r_m their radius"""
    partial, total = False, False

    earth_sun = S - E
    moon_sun = S - M

    # The moon must be between the earth and the sun
    if norm(moon_sun) > norm(earth_sun):
        return partial, total

    # The point of the moon's sun-facing plane which is on the sun-earth line
    earth_on_moon_plane = intersectionPoint(
        M, moon_sun, S, E
    )

    # Moon to Earth's Projection direction
    MC_EP = earth_on_moon_plane - M
    MC_EP = MC_EP / norm(MC_EP)

    # The lines pass by the same point of the moon
    moon_point = M + r_m * MC_EP


    # The lines pass by opposite points of the sun
    sun_point_partial = S - r_s * MC_EP
    sun_point_total = S + r_s * MC_EP

    # moon_point is the point of the moon that is the closest to the sun-earth axis
    # hence the bit of shadow that is the closest to the earth's center is on the line that
    # passes by sun_point_partial (for partial eclipse) and moon_point
    # or by sun_point_total (for total eclipse) and moon_point

    shadow_for_partial = intersectionPoint(
        E, earth_sun, sun_point_partial, moon_point
    )

    shadow_for_total = intersectionPoint(
        E, earth_sun, sun_point_total, moon_point
    )

    if distance(shadow_for_partial, E) < r_e:
        partial = True

    if distance(shadow_for_total, E) < r_e:
        total = True

    return partial, total

Approximating the solar system

Now that we can tell if there’s an eclipse in a given position, we need to compute all future positions given the current (observable) one.

The odeint python function can solve differential equations looking like $y'(t) = f(y(t), t)$ for a set of instants $t$ given initial conditions $y(0) = y_0$ .

Let’s use $p_i: t \mapsto \begin{pmatrix} x_i(t) \\ y_i(t) \\ z_i(t) \end{pmatrix}$ as the position function for the body $i$.

Using the fundamental principle of dynamic, for all body $i$ we have:

That gives us a way to compute the derivative of the vector $\begin{pmatrix}p_i(t) \\ p_i'(t)\end{pmatrix}$:

Knowing the derivative of the vector at a given time and its value for $t = 0$, that is the current measurable positions and speed of the celestial bodies, we can use odeint to compute their future positions.

In practice, our $f$ function will take a flat vector of all bodies positions and speeds and will output it’s derivative.

def f(y, t):
    """y: flat array of bodies [x1,y1,z1,dx1,dy1,dz1, x2, y2, ...]
       returns a flat array
       to be used in odeint"""
    s = np.zeros((nb_bodies, 6))
    y = y.reshape((nb_bodies, 6))
    for i in range(nb_bodies):
        for j in range(0, nb_bodies):
            if i != j:
                a = y[j][0:3] - y[i][0:3] # p_j - p_i
                d3 = sum(a**2)**1.5  # squarred then cubed distance
                s[i][3:6] += a / d3 * masses[j]
        s[i][3:6] *= G
        s[i][0:3] = y[i][3:6]
    y = y.reshape(nb_bodies * 6)
    return s.reshape(nb_bodies * 6)

# an instant every seconds for 3 days
instants = np.linspace(0, 3, 3600 * 24 * 3 + 1)
positions = sint.odeint(f, init_bodies, intstants)

Now just check if each position is a partial or a total eclipse ! My result is:

Partial Eclipse
Beginning:  2016-03-08T23:20:06
End:        2016-03-09T04:35:30

Total Eclipse
Beginning:  2016-03-09T00:16:41
End:        2016-03-09T03:39:00

which is very close (< 1min error !) to the actual timestamps from the french’s wikipedia of the eclipse.

The full code with the necessary init values is available here.

Since this method gives the position of the eclipse on the earth’s disk, it can be extended to infer where the eclipse will be visible from.