# Eclipse prediction

# 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.

- To compute
`MC_EP`

, between- the sun-earth axis and
- the moon’s plane

- 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

- the line $(rt)$, originating from the perimeter of the sun’s disk in the direction of
- 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

- the line $(rp)$, originating from the perimeter of the sun’s disk in the opposite direction of

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.