{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Eclipse 2024\n", "\n", "For the complete solar eclipse in North America on April 8, 2024:\n", "\n", "* Generate eclipse statistic for any given location\n", "* Compute and plot the centerline of totalitys" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Compute & plot center line of eclipse totality" ] }, { "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 plotly.graph_objects as go\n", "\n", "lat, lon = zip(*[(c.latitude_deg, c.longitude_deg) for c in coords])\n", "\n", "fig = go.Figure()\n", "fig.add_trace(\n", " go.Scattergeo(\n", " locationmode = 'USA-states',\n", " lon = lon,\n", " lat = lat,\n", " mode = 'lines',\n", " line = dict(width = 2,color = 'red'),\n", " )\n", ")\n", "fig.update_layout(\n", " title_text = 'Centerline of 2024 Solar Eclipse',\n", " showlegend = False,\n", " geo = dict(\n", " projection_type = 'mercator',\n", " resolution=50,\n", " lataxis_range=[20,55],\n", " lonaxis_range=[-140,-40],\n", " showcountries=True,\n", " countrycolor='black',\n", " ),\n", " width=650,\n", " height=500,\n", ")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Eclipse statistics for arbitrary locations on Earth" ] }, { "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]].datetime(), # type: ignore\n", " \"end\": timearr[eidx[-1][0]].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]].datetime(),\n", " \"end\": timearr[pidx[-1][0]].datetime(),\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]].datetime(), # type: ignore\n", " \"end\": timearr[pidx[-1][0]].datetime(), # type: ignore\n", " \"peak\": timearr[idx].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 }