{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Eclipse 2024\n",
"\n",
"This tutorial analyzes the total solar eclipse that crossed North America on April 8, 2024. It demonstrates two capabilities:\n",
"\n",
"1. **Centerline computation**: Project the Moon's shadow onto the Earth's surface using JPL ephemerides for the Sun and Moon, accounting for Earth's oblateness\n",
"2. **Eclipse statistics**: For any ground location, compute the timing and extent of partial and total eclipse phases by comparing the angular separation of the Sun and Moon against their angular radii"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compute & Plot Centerline of Eclipse Totality\n",
"\n",
"The centerline is found by projecting the Moon onto the Earth's surface along the Sun direction. Since the Earth is an oblate spheroid, we scale the z-axis by $1/(1-f)$ where $f$ is the WGS-84 flattening factor, solve for the ray-sphere intersection using the law of cosines, then undo the scaling to recover true ITRF coordinates."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import satkit as sk\n",
"import math as m\n",
"\n",
"# Get an array of times over the day\n",
"theday = sk.time(2024, 4, 8)\n",
"timearr = np.array([theday + sk.duration.from_seconds(i) for i in range(0, 86400, 10)])\n",
"\n",
"# Get the sun position at each time\n",
"sun_gcrf = np.array([sk.jplephem.geocentric_pos(sk.solarsystem.Sun, t) for t in timearr])\n",
"moon_gcrf = np.array([sk.jplephem.geocentric_pos(sk.solarsystem.Moon, t) for t in timearr])\n",
"# Rotate to Earth-fixed ITRF coordinates\n",
"qarray = sk.frametransform.qgcrf2itrf(timearr)\n",
"sun_itrf = np.array([q*x for q,x in zip(qarray, sun_gcrf)])\n",
"moon_itrf = np.array([q*x for q,x in zip(qarray, moon_gcrf)])\n",
"\n",
"# Project the moon to the surface of the Earth along the Sun vector\n",
"# This is projecting onto the surface of a sphere. However, the Earth is actually\n",
"# an oblate spheroid. We can account for this by scaling the z axis by the flattening\n",
"# of the Earth.\n",
"scalefac = 1.0 / (1.0 - sk.consts.wgs84_f)\n",
"moon_itrf_scaled = moon_itrf\n",
"moon_itrf_scaled[:,2] = moon_itrf_scaled[:,2] * scalefac\n",
"sun_itrf_scaled = sun_itrf\n",
"sun_itrf_scaled[:,2] = sun_itrf_scaled[:,2] * scalefac\n",
"sun_itrf_scaled_hat = sun_itrf_scaled / np.linalg.norm(sun_itrf_scaled, axis=1)[:,None]\n",
"\n",
"# Compute the distance to the surface of the Earth. This can be done\n",
"# via the law of Cosines and the quadratic equation\n",
"lcostheta = np.sum(sun_itrf_scaled_hat * moon_itrf_scaled, axis=1)\n",
"sqrtterm = lcostheta**2 - np.sum(moon_itrf_scaled**2, axis=1) + sk.consts.earth_radius**2\n",
"# Valid indices (where the projection hits the earth) are where the term under\n",
"# the square root is positive\n",
"vidx = np.argwhere(sqrtterm > 0).flatten()\n",
"\n",
"# Distance to surface of the Earth along from the moon along the Sun vector\n",
"dist = lcostheta[vidx] - np.sqrt(sqrtterm[vidx])\n",
"# Now, get the positions on the surface of the Earth\n",
"pgnd_itrf_scaled = moon_itrf_scaled[vidx] - dist[:,None] * sun_itrf_scaled_hat[vidx]\n",
"# Undo the scaling of the zaxis to account for Earth flattening\n",
"pgnd_itrf = pgnd_itrf_scaled\n",
"pgnd_itrf[:,2] = pgnd_itrf[:,2] / scalefac\n",
"\n",
"# Get the latitude and longitude of the points on the surface of the Earth\n",
"# that follow the center line of totality\n",
"coords = [sk.itrfcoord(x) for x in pgnd_itrf]\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import ssl, certifi\n",
"ssl._create_default_https_context = lambda: ssl.create_default_context(cafile=certifi.where())\n",
"import matplotlib.pyplot as plt\n",
"import scienceplots # noqa: F401\n",
"plt.style.use([\"science\", \"no-latex\", \"../satkit.mplstyle\"])\n",
"%config InlineBackend.figure_formats = ['svg']\n",
"import warnings\n",
"warnings.filterwarnings(\"ignore\", \"Downloading\")\n",
"import cartopy.crs as ccrs\n",
"import cartopy.feature as cfeature\n",
"\n",
"lat, lon = zip(*[(c.latitude_deg, c.longitude_deg) for c in coords])\n",
"\n",
"fig, ax = plt.subplots(figsize=(10, 6), subplot_kw={'projection': ccrs.Mercator()})\n",
"ax.set_extent([-140, -40, 20, 55], crs=ccrs.PlateCarree())\n",
"ax.add_feature(cfeature.LAND, facecolor='lightgray')\n",
"ax.add_feature(cfeature.BORDERS, linewidth=0.5)\n",
"ax.add_feature(cfeature.COASTLINE, linewidth=0.5)\n",
"ax.plot(lon, lat, color='red', linewidth=2, transform=ccrs.PlateCarree())\n",
"ax.set_title(\"Centerline of 2024 Solar Eclipse\")\n",
"plt.tight_layout()\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Eclipse Statistics for Arbitrary Locations\n",
"\n",
"For any observer, we compute the angular separation between the Sun and Moon over time and compare it to their angular radii:\n",
"\n",
"- **Total eclipse** occurs when the separation is less than the difference of the angular radii (Moon fully covers the Sun)\n",
"- **Partial eclipse** occurs when the separation is less than the sum of the angular radii (discs overlap)\n",
"- For partial-only locations, the area and diameter fractions of the Sun that are occluded are computed using the geometry of two overlapping circles (lens area via law of cosines)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import satkit as sk\n",
"import numpy as np\n",
"import math as m\n",
"\n",
"# Eclipse happens on April 8, 2024\n",
"time0 = sk.time(2024, 4, 8, 12, 0, 0)\n",
"timearr = np.array(\n",
" time0 + [sk.duration.from_days(x) for x in np.linspace(0, 0.5, 43200)]\n",
")\n",
"\n",
"# Get exact JPL ephemeris for sun & moon\n",
"sun_light_travel_time = sk.duration.from_seconds(sk.consts.au / sk.consts.c)\n",
"sun_gcrf = sk.jplephem.geocentric_pos(\n",
" sk.solarsystem.Sun, timearr - sun_light_travel_time\n",
")\n",
"moon_gcrf = sk.jplephem.geocentric_pos(sk.solarsystem.Moon, timearr)\n",
"# Rotation to Earth-fixed frame\n",
"qitrf2gcrf = sk.frametransform.qitrf2gcrf(timearr)\n",
"\n",
"def eclipse_stats(loc: sk.itrfcoord):\n",
" qitrf2ned = loc.qned2itrf.conj\n",
"\n",
" # Location in GCRF\n",
" loc_gcrf = np.array([x * loc.vector for x in qitrf2gcrf])\n",
"\n",
" # Compute angle between sun and moon at location\n",
" sun_diff = sun_gcrf - loc_gcrf\n",
" moon_diff = moon_gcrf - loc_gcrf\n",
" sun_norm = np.sqrt(np.sum(sun_diff**2, axis=1))\n",
" moon_norm = np.sqrt(np.sum(moon_diff**2, axis=1))\n",
" theta = np.arccos(np.sum(sun_diff * moon_diff, axis=1) / sun_norm / moon_norm)\n",
"\n",
" # Compute angular extent of sun & moon\n",
" moon_dist = np.mean(moon_norm)\n",
" moon_extent_rad = sk.consts.moon_radius / moon_dist\n",
" sun_extent_rad = sk.consts.sun_radius / sk.consts.au\n",
" # How far off can they be while still having total eclipse?\n",
" max_eclipse_offset_rad = moon_extent_rad - sun_extent_rad\n",
"\n",
" idx = np.argwhere(theta == np.min(theta))[0][0]\n",
" # Look for times where there is total eclipse\n",
" eidx = np.argwhere(theta < max_eclipse_offset_rad)\n",
" # Look for times of partial eclipse\n",
" pidx = np.argwhere(theta < (sun_extent_rad + moon_extent_rad))\n",
"\n",
" data = {\"latitude\": loc.latitude_deg, \"longitude\": loc.longitude_deg}\n",
"\n",
" if len(eidx) > 0:\n",
" data[\"total\"] = {\n",
" \"begin\": timearr[eidx[0][0]].as_datetime(), # type: ignore\n",
" \"end\": timearr[eidx[-1][0]].as_datetime(), # type: ignore\n",
" \"duration_seconds\": (timearr[eidx[-1][0]] - timearr[eidx[0][0]]).seconds, # type: ignore\n",
" }\n",
" data[\"partial\"] = {\n",
" \"begin\": timearr[pidx[0][0]].as_datetime(), # type: ignore\n",
" \"end\": timearr[pidx[-1][0]].as_datetime(), # type: ignore\n",
" \"duration_seconds\": (timearr[pidx[-1][0]] - timearr[pidx[0][0]]).seconds, # type: ignore\n",
" \"peak\": None,\n",
" \"minangle_deg\": None,\n",
" \"max_area_occlusion\": None,\n",
" \"max_diam_occlusion\": None,\n",
" }\n",
" elif np.min(theta) < (sun_extent_rad + moon_extent_rad):\n",
" durp = timearr[pidx[-1]][0] - timearr[pidx[0]][0]\n",
" mintheta = np.min(theta)\n",
" # Derived via triangles & law of cosines\n",
" theta_a = m.acos(\n",
" (sun_extent_rad**2 + mintheta**2 - moon_extent_rad**2)\n",
" / (2 * sun_extent_rad * mintheta)\n",
" )\n",
" theta_b = m.acos(\n",
" (moon_extent_rad**2 + mintheta**2 - sun_extent_rad**2)\n",
" / (2 * moon_extent_rad * mintheta)\n",
" )\n",
" h = sun_extent_rad * m.sin(theta_a)\n",
" Lb = h / m.tan(theta_b)\n",
" La = h / m.tan(theta_a)\n",
" # Area of right side of overlapping \"lens\"\n",
" aright = m.pi * moon_extent_rad**2 * theta_b / m.pi - Lb * h\n",
" # Area of left side of overlapping \"lens\"\n",
" aleft = m.pi * sun_extent_rad**2 * theta_a / m.pi - La * h\n",
" ashown = m.pi * sun_extent_rad**2 - aright - aleft\n",
" max_frac_area_occluded = 1 - ashown / (m.pi * sun_extent_rad**2)\n",
" max_frac_diam_occluded = 1 - (sun_extent_rad + mintheta - moon_extent_rad) / (\n",
" 2 * sun_extent_rad\n",
" )\n",
"\n",
" data[\"partial\"] = {\n",
" \"begin\": timearr[pidx[0][0]].as_datetime(), # type: ignore\n",
" \"end\": timearr[pidx[-1][0]].as_datetime(), # type: ignore\n",
" \"peak\": timearr[idx].as_datetime(), # type: ignore\n",
" \"duration_seconds\": (timearr[pidx[-1][0]] - timearr[pidx[0][0]]).seconds, # type: ignore\n",
" \"minangle_deg\": np.min(theta) * 180.0 / m.pi,\n",
" \"max_area_occlusion\": max_frac_area_occluded,\n",
" \"max_diam_occlusion\": max_frac_diam_occluded,\n",
" }\n",
" data[\"total\"] = None\n",
" else:\n",
" data[\"total\"] = None\n",
" data[\"partial\"] = None\n",
"\n",
" return data"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# A list of eclipse locations\n",
"\n",
"locations = [\n",
" { \"name\": \"Mexico City\", \"latitude_deg\": 19.4326, \"longitude_deg\": -99.1332, },\n",
" { \"name\": \"Austin, Tx\", \"latitude_deg\": 30.2672, \"longitude_deg\": -97.7431, },\n",
" { \"name\": \"Dallas, Tx\", \"latitude_deg\": 32.7767, \"longitude_deg\": -96.7970, },\n",
" { \"name\": \"St. Louis, Mo\", \"latitude_deg\": 38.6270, \"longitude_deg\": -90.1994, },\n",
" { \"name\": \"New York City\", \"latitude_deg\": 40.7128, \"longitude_deg\": -74.0060, },\n",
" { \"name\": \"Boston, MA\", \"latitude_deg\": 42.3601, \"longitude_deg\": -71.0589, },\n",
" { \"name\": \"Burlington, Vt\", \"latitude_deg\": 44.4759, \"longitude_deg\": -73.2121, },\n",
" { \"name\": \"Montreal\", \"latitude_deg\": 45.5017, \"longitude_deg\": -73.5673, },\n",
" { \"name\": \"Quebec City\", \"latitude_deg\": 46.8139, \"longitude_deg\": -71.2080, },\n",
" { \"name\": \"Halifax\", \"latitude_deg\": 44.6488, \"longitude_deg\": -63.5752, },\n",
" { \"name\": \"Cleveland, Oh\", \"latitude_deg\": 41.4993, \"longitude_deg\": -81.6944, },\n",
"]\n",
"\n",
"for loc in locations:\n",
" loc[\"coord\"] = sk.itrfcoord(latitude_deg=loc[\"latitude_deg\"], longitude_deg=loc[\"longitude_deg\"])\n",
" loc[\"stats\"] = eclipse_stats(loc[\"coord\"])"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Display stats in a nice table\n",
"import pandas as pd\n",
"\n",
"df = pd.DataFrame(\n",
" [\n",
" {\n",
" \"Location\": loc[\"name\"],\n",
" \"Total Eclipse Begin\": loc[\"stats\"][\"total\"][\"begin\"] if loc[\"stats\"][\"total\"] else None,\n",
" \"Total Eclipse End\": loc[\"stats\"][\"total\"][\"end\"] if loc[\"stats\"][\"total\"] else None,\n",
" \"Total Eclipse Duration (s)\": loc[\"stats\"][\"total\"][\"duration_seconds\"] if loc[\"stats\"][\"total\"] else None,\n",
" \"Partial Eclipse Begin\": loc[\"stats\"][\"partial\"][\"begin\"] if loc[\"stats\"][\"partial\"] else None,\n",
" \"Partial Eclipse End\": loc[\"stats\"][\"partial\"][\"end\"] if loc[\"stats\"][\"partial\"] else None,\n",
" \"Partial Eclipse Peak\": loc[\"stats\"][\"partial\"][\"peak\"] if loc[\"stats\"][\"partial\"] else None,\n",
" \"Partial Eclipse Duration (s)\": loc[\"stats\"][\"partial\"][\"duration_seconds\"] if loc[\"stats\"][\"partial\"] else None,\n",
" \"Min Separation (deg)\": loc[\"stats\"][\"partial\"][\"minangle_deg\"] if loc[\"stats\"][\"partial\"] else None,\n",
" \"Max Area Occlusion\": loc[\"stats\"][\"partial\"][\"max_area_occlusion\"] if loc[\"stats\"][\"partial\"] else None,\n",
" \"Max Diameter Occlusion\": loc[\"stats\"][\"partial\"][\"max_diam_occlusion\"] if loc[\"stats\"][\"partial\"] else None,\n",
" }\n",
" for loc in locations\n",
" ]\n",
")\n",
"df.style.format({\"Total Eclipse Begin\": lambda x: x.strftime(\"%H:%M:%S\") if not pd.isnull(x) else 'N/A',\n",
" \"Total Eclipse End\": lambda x: x.strftime(\"%H:%M:%S\") if not pd.isnull(x) else 'N/A',\n",
" \"Partial Eclipse Begin\": lambda x: x.strftime(\"%H:%M:%S\") if not pd.isnull(x) else 'N/A',\n",
" \"Partial Eclipse End\": lambda x: x.strftime(\"%H:%M:%S\") if not pd.isnull(x) else 'N/A',\n",
" \"Partial Eclipse Peak\": lambda x: x.strftime(\"%H:%M:%S\") if not pd.isnull(x) else 'N/A',\n",
" \"Total Eclipse Duration (s)\": lambda x: f\"{x:.0f}\" if not pd.isnull(x) else 'N/A',\n",
" \"Partial Eclipse Duration (s)\": lambda x: f\"{x:.0f}\" if not pd.isnull(x) else 'N/A',\n",
" \"Min Separation (deg)\": lambda x: f\"{x:.2f}\" if not pd.isnull(x) else \"N/A\",\n",
" \"Max Area Occlusion\": lambda x: f\"{x:.2f}\" if not pd.isnull(x) else '1.00',\n",
" \"Max Diameter Occlusion\": lambda x: f\"{x:.2f}\" if not pd.isnull(x) else '1.00',\n",
" })"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.13.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}