Eclipse 2024
For the complete solar eclipse in North America on April 8, 2024:
Generate eclipse statistic for any given location
Compute and plot the centerline of totalitys
Compute & plot center line of eclipse totality
[1]:
import numpy as np
import satkit as sk
import math as m
# Get an array of times over the day
theday = sk.time(2024, 4, 8)
timearr = np.array([theday + sk.duration.from_seconds(i) for i in range(0, 86400, 10)])
# Get the sun position at each time
sun_gcrf = np.array([sk.jplephem.geocentric_pos(sk.solarsystem.Sun, t) for t in timearr])
moon_gcrf = np.array([sk.jplephem.geocentric_pos(sk.solarsystem.Moon, t) for t in timearr])
# Rotate to Earth-fixed ITRF coordinates
qarray = sk.frametransform.qgcrf2itrf(timearr)
sun_itrf = np.array([q*x for q,x in zip(qarray, sun_gcrf)])
moon_itrf = np.array([q*x for q,x in zip(qarray, moon_gcrf)])
# Project the moon to the surface of the Earth along the Sun vector
# This is projecting onto the surface of a sphere. However, the Earth is actually
# an oblate spheroid. We can account for this by scaling the z axis by the flattening
# of the Earth.
scalefac = 1.0 / (1.0 - sk.consts.wgs84_f)
moon_itrf_scaled = moon_itrf
moon_itrf_scaled[:,2] = moon_itrf_scaled[:,2] * scalefac
sun_itrf_scaled = sun_itrf
sun_itrf_scaled[:,2] = sun_itrf_scaled[:,2] * scalefac
sun_itrf_scaled_hat = sun_itrf_scaled / np.linalg.norm(sun_itrf_scaled, axis=1)[:,None]
# Compute the distance to the surface of the Earth. This can be done
# via the law of Cosines and the quadratic equation
lcostheta = np.sum(sun_itrf_scaled_hat * moon_itrf_scaled, axis=1)
sqrtterm = lcostheta**2 - np.sum(moon_itrf_scaled**2, axis=1) + sk.consts.earth_radius**2
# Valid indices (where the projection hits the earth) are where the term under
# the square root is positive
vidx = np.argwhere(sqrtterm > 0).flatten()
# Distance to surface of the Earth along from the moon along the Sun vector
dist = lcostheta[vidx] - np.sqrt(sqrtterm[vidx])
# Now, get the positions on the surface of the Earth
pgnd_itrf_scaled = moon_itrf_scaled[vidx] - dist[:,None] * sun_itrf_scaled_hat[vidx]
# Undo the scaling of the zaxis to account for Earth flattening
pgnd_itrf = pgnd_itrf_scaled
pgnd_itrf[:,2] = pgnd_itrf[:,2] / scalefac
# Get the latitude and longitude of the points on the surface of the Earth
# that follow the center line of totality
coords = [sk.itrfcoord(x) for x in pgnd_itrf]
Downloading JPL Ephemeris file. File size is approx. 100MB
[2]:
import plotly.graph_objects as go
lat, lon = zip(*[(c.latitude_deg, c.longitude_deg) for c in coords])
fig = go.Figure()
fig.add_trace(
go.Scattergeo(
locationmode = 'USA-states',
lon = lon,
lat = lat,
mode = 'lines',
line = dict(width = 2,color = 'red'),
)
)
fig.update_layout(
title_text = 'Centerline of 2024 Solar Eclipse',
showlegend = False,
geo = dict(
projection_type = 'mercator',
resolution=50,
lataxis_range=[20,55],
lonaxis_range=[-140,-40],
showcountries=True,
countrycolor='black',
),
width=650,
height=500,
)
Eclipse statistics for arbitrary locations on Earth
[3]:
import satkit as sk
import numpy as np
import math as m
# Eclipse happens on April 8, 2024
time0 = sk.time(2024, 4, 8, 12, 0, 0)
timearr = np.array(
time0 + [sk.duration.from_days(x) for x in np.linspace(0, 0.5, 43200)]
)
# Get exact JPL ephemeris for sun & moon
sun_light_travel_time = sk.duration.from_seconds(sk.consts.au / sk.consts.c)
sun_gcrf = sk.jplephem.geocentric_pos(
sk.solarsystem.Sun, timearr - sun_light_travel_time
)
moon_gcrf = sk.jplephem.geocentric_pos(sk.solarsystem.Moon, timearr)
# Rotation to Earth-fixed frame
qitrf2gcrf = sk.frametransform.qitrf2gcrf(timearr)
def eclipse_stats(loc: sk.itrfcoord):
qitrf2ned = loc.qned2itrf.conj
# Location in GCRF
loc_gcrf = np.array([x * loc.vector for x in qitrf2gcrf])
# Compute angle between sun and moon at location
sun_diff = sun_gcrf - loc_gcrf
moon_diff = moon_gcrf - loc_gcrf
sun_norm = np.sqrt(np.sum(sun_diff**2, axis=1))
moon_norm = np.sqrt(np.sum(moon_diff**2, axis=1))
theta = np.arccos(np.sum(sun_diff * moon_diff, axis=1) / sun_norm / moon_norm)
# Compute angular extent of sun & moon
moon_dist = np.mean(moon_norm)
moon_extent_rad = sk.consts.moon_radius / moon_dist
sun_extent_rad = sk.consts.sun_radius / sk.consts.au
# How far off can they be while still having total eclipse?
max_eclipse_offset_rad = moon_extent_rad - sun_extent_rad
idx = np.argwhere(theta == np.min(theta))[0][0]
# Look for times where there is total eclipse
eidx = np.argwhere(theta < max_eclipse_offset_rad)
# Look for times of partial eclipse
pidx = np.argwhere(theta < (sun_extent_rad + moon_extent_rad))
data = {"latitude": loc.latitude_deg, "longitude": loc.longitude_deg}
if len(eidx) > 0:
data["total"] = {
"start": timearr[eidx[0][0]].datetime(),
"stop": timearr[eidx[-1][0]].datetime(),
"duration_seconds": (timearr[eidx[-1][0]] - timearr[eidx[0][0]]).seconds(),
}
data["partial"] = {
"start": timearr[pidx[0][0]].datetime(),
"stop": timearr[pidx[-1][0]].datetime(),
"duration_seconds": (timearr[pidx[-1][0]] - timearr[pidx[0][0]]).seconds(),
"peak": None,
"minangle_deg": None,
"max_area_occlusion": None,
"max_diam_occlusion": None,
}
elif np.min(theta) < (sun_extent_rad + moon_extent_rad):
durp = timearr[pidx[-1]][0] - timearr[pidx[0]][0]
mintheta = np.min(theta)
# Derived via traingles & law of cosines
theta_a = m.acos(
(sun_extent_rad**2 + mintheta**2 - moon_extent_rad**2)
/ (2 * sun_extent_rad * mintheta)
)
theta_b = m.acos(
(moon_extent_rad**2 + mintheta**2 - sun_extent_rad**2)
/ (2 * moon_extent_rad * mintheta)
)
h = sun_extent_rad * m.sin(theta_a)
Lb = h / m.tan(theta_b)
La = h / m.tan(theta_a)
# Area of right side of overlapping "lens"
aright = m.pi * moon_extent_rad**2 * theta_b / m.pi - Lb * h
# Area of left side of overlapping "lens"
aleft = m.pi * sun_extent_rad**2 * theta_a / m.pi - La * h
ashown = m.pi * sun_extent_rad**2 - aright - aleft
max_frac_area_occluded = 1 - ashown / (m.pi * sun_extent_rad**2)
max_frac_diam_occluded = 1 - (sun_extent_rad + mintheta - moon_extent_rad) / (
2 * sun_extent_rad
)
data["partial"] = {
"start": timearr[pidx[0][0]].datetime(),
"stop": timearr[pidx[-1][0]].datetime(),
"peak": timearr[idx].datetime(),
"duration_seconds": (timearr[pidx[-1][0]] - timearr[pidx[0][0]]).seconds(),
"minangle_deg": np.min(theta) * 180.0 / m.pi,
"max_area_occlusion": max_frac_area_occluded,
"max_diam_occlusion": max_frac_diam_occluded,
}
data["total"] = None
else:
data["total"] = None
data["partial"] = None
return data
[4]:
# A list of eclipse locations
locations = [
{ "name": "Mexico City", "latitude_deg": 19.4326, "longitude_deg": -99.1332, },
{ "name": "Austin, Tx", "latitude_deg": 30.2672, "longitude_deg": -97.7431, },
{ "name": "Dallas, Tx", "latitude_deg": 32.7767, "longitude_deg": -96.7970, },
{ "name": "St. Louis, Mo", "latitude_deg": 38.6270, "longitude_deg": -90.1994, },
{ "name": "New York City", "latitude_deg": 40.7128, "longitude_deg": -74.0060, },
{ "name": "Boston, MA", "latitude_deg": 42.3601, "longitude_deg": -71.0589, },
{ "name": "Burlington, Vt", "latitude_deg": 44.4759, "longitude_deg": -73.2121, },
{ "name": "Montreal", "latitude_deg": 45.5017, "longitude_deg": -73.5673, },
{ "name": "Quebec City", "latitude_deg": 46.8139, "longitude_deg": -71.2080, },
{ "name": "Halifax", "latitude_deg": 44.6488, "longitude_deg": -63.5752, },
{ "name": "Cleveland, Oh", "latitude_deg": 41.4993, "longitude_deg": -81.6944, },
]
for loc in locations:
loc["coord"] = sk.itrfcoord(latitude_deg=loc["latitude_deg"], longitude_deg=loc["longitude_deg"])
loc["stats"] = eclipse_stats(loc["coord"])
[5]:
# Display stats in a nice table
import pandas as pd
df = pd.DataFrame(
[
{
"Location": loc["name"],
"Total Eclipse Start": loc["stats"]["total"]["start"] if loc["stats"]["total"] else None,
"Total Eclipse Stop": loc["stats"]["total"]["stop"] if loc["stats"]["total"] else None,
"Total Eclipse Duration (s)": loc["stats"]["total"]["duration_seconds"] if loc["stats"]["total"] else None,
"Partial Eclipse Start": loc["stats"]["partial"]["start"] if loc["stats"]["partial"] else None,
"Partial Eclipse Stop": loc["stats"]["partial"]["stop"] if loc["stats"]["partial"] else None,
"Partial Eclipse Peak": loc["stats"]["partial"]["peak"] if loc["stats"]["partial"] else None,
"Partial Eclipse Duration (s)": loc["stats"]["partial"]["duration_seconds"] if loc["stats"]["partial"] else None,
"Min Seperation (deg)": loc["stats"]["partial"]["minangle_deg"] if loc["stats"]["partial"] else None,
"Max Area Occlusion": loc["stats"]["partial"]["max_area_occlusion"] if loc["stats"]["partial"] else None,
"Max Diameter Occlusion": loc["stats"]["partial"]["max_diam_occlusion"] if loc["stats"]["partial"] else None,
}
for loc in locations
]
)
df.style.format({"Total Eclipse Start": lambda x: x.strftime("%H:%M:%S") if not pd.isnull(x) else 'N/A',
"Total Eclipse Stop": lambda x: x.strftime("%H:%M:%S") if not pd.isnull(x) else 'N/A',
"Partial Eclipse Start": lambda x: x.strftime("%H:%M:%S") if not pd.isnull(x) else 'N/A',
"Partial Eclipse Stop": lambda x: x.strftime("%H:%M:%S") if not pd.isnull(x) else 'N/A',
"Partial Eclipse Peak": lambda x: x.strftime("%H:%M:%S") if not pd.isnull(x) else 'N/A',
"Total Eclipse Duration (s)": lambda x: f"{x:.0f}" if not pd.isnull(x) else 'N/A',
"Partial Eclipse Duration (s)": lambda x: f"{x:.0f}" if not pd.isnull(x) else 'N/A',
"Min Seperation (deg)": lambda x: f"{x:.2f}" if not pd.isnull(x) else "N/A",
"Max Area Occlusion": lambda x: f"{x:.2f}" if not pd.isnull(x) else '1.00',
"Max Diameter Occlusion": lambda x: f"{x:.2f}" if not pd.isnull(x) else '1.00',
})
[5]:
| Location | Total Eclipse Start | Total Eclipse Stop | Total Eclipse Duration (s) | Partial Eclipse Start | Partial Eclipse Stop | Partial Eclipse Peak | Partial Eclipse Duration (s) | Min Seperation (deg) | Max Area Occlusion | Max Diameter Occlusion | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | Mexico City | N/A | N/A | N/A | 16:55:32 | 19:36:12 | 18:14:15 | 9640 | 0.13 | 0.74 | 0.79 |
| 1 | Austin, Tx | N/A | N/A | N/A | 17:17:23 | 19:57:55 | 18:37:00 | 9632 | 0.01 | 1.00 | 1.00 |
| 2 | Dallas, Tx | 18:40:57 | 18:44:16 | 199 | 17:23:28 | 20:02:30 | N/A | 9542 | N/A | 1.00 | 1.00 |
| 3 | St. Louis, Mo | N/A | N/A | N/A | 17:43:08 | 20:17:15 | 19:00:48 | 9247 | 0.02 | 0.99 | 0.99 |
| 4 | New York City | N/A | N/A | N/A | 18:10:44 | 20:36:25 | 19:25:34 | 8741 | 0.06 | 0.90 | 0.91 |
| 5 | Boston, MA | N/A | N/A | N/A | 18:16:14 | 20:39:10 | 19:29:47 | 8576 | 0.05 | 0.92 | 0.93 |
| 6 | Burlington, Vt | 19:26:12 | 19:29:15 | 183 | 18:14:22 | 20:37:21 | N/A | 8579 | N/A | 1.00 | 1.00 |
| 7 | Montreal | 19:27:03 | 19:27:59 | 56 | 18:14:35 | 20:36:52 | N/A | 8537 | N/A | 1.00 | 1.00 |
| 8 | Quebec City | N/A | N/A | N/A | 18:18:27 | 20:38:24 | 19:30:15 | 8397 | 0.02 | 0.99 | 0.98 |
| 9 | Halifax | N/A | N/A | N/A | 18:27:26 | 20:44:13 | 19:38:02 | 8207 | 0.04 | 0.94 | 0.95 |
| 10 | Cleveland, Oh | 19:13:52 | 19:17:24 | 212 | 17:59:31 | 20:28:57 | N/A | 8966 | N/A | 1.00 | 1.00 |