Author: Hal Snyder
Views : 29
Description: 2024 Solar Eclipse demo
%auto # This cell automatically evaluates on startup -- or run it manually if it didn't evaluate. # Here, it initializes the Jupyter kernel with the specified name and sets it as the default mode for this worksheet. jupyter_kernel = jupyter("python3") %default_mode jupyter_kernel
%sage sage_server.MAX_HTML_SIZE = 1000000 jupyter_kernel.smc_image_scaling = 0.5

Studying the 2024 April 8 Solar Eclipse with Astropy

Goals

• Learn about the 2024 Eclipse.

• Model it with Astropy.

• Check our work against the USNO.

Extra Sources Step 1. Get basic info from USNO.

(Ok, this is cheating.) Step 2. Follow the steps of Adrian Price-Whelan's 2017 notebook.

This will give us animations and the setup for basic stats.

Add geopy so it knows about the Palatine Library.

from datetime import datetime import astropy.coordinates as coord from astropy.time import Time import astropy.units as u import numpy as np import pytz from IPython.display import HTML import matplotlib.pyplot as plt from matplotlib.dates import HourLocator, MinuteLocator, DateFormatter import matplotlib.animation as animation import matplotlib as mpl import math
mpl.rcParams['figure.figsize'] = (8,6) mpl.rcParams['axes.titlesize'] = 26 mpl.rcParams['axes.labelsize'] = 22 mpl.rcParams['xtick.labelsize'] = 18 mpl.rcParams['ytick.labelsize'] = 18
# Set timezone here: mpl.rcParams['timezone'] = 'US/Central' # Enter address here: address = 'Palatine Public Library, Palatine, IL'

Choose a range of 1024 equally-spaced times that span the eclipse

17:00 UTC to 21:00 UTC

tz = pytz.timezone(mpl.rcParams['timezone']) times = (Time(datetime(2024, 4, 8, 17, 0, 0).astimezone(pytz.UTC)) + np.linspace(0, 4, 1024) * u.hour) dt = times.to_datetime(tz)
WARNING: ErfaWarning: ERFA function "dtf2d" yielded 1 of "dubious year (Note 6)" [astropy._erfa.core] WARNING: ErfaWarning: ERFA function "utctai" yielded 1 of "dubious year (Note 3)" [astropy._erfa.core] WARNING: ErfaWarning: ERFA function "taiutc" yielded 1024 of "dubious year (Note 4)" [astropy._erfa.core] WARNING: ErfaWarning: ERFA function "d2dtf" yielded 1024 of "dubious year (Note 5)" [astropy._erfa.core]
# geopy gets street address and county for the library as well as long/lat. Cool. import geopy from geopy.geocoders import Nominatim geolocator = Nominatim(user_agent="chi-eclipse") location = geolocator.geocode("Palatine Public Library, Palatine, IL, US") print(location.address) print(f'lon: {location.longitude}, lat: {location.latitude}')
Palatine Public Library, 700, North Court, Palatine, Cook County, Illinois, 60067, USA lon: -88.0366723932833, lat: 42.12225845
location.altitude
0.0
loc = coord.EarthLocation(lon=location.longitude*u.deg, lat=location.latitude*u.deg, height=0*u.m)
# use the astropy ephemeris to get coordinates of moon and sun # for the range of times moon = coord.get_moon(times) sun = coord.get_sun(times)
WARNING: ErfaWarning: ERFA function "utctai" yielded 1024 of "dubious year (Note 3)" [astropy._erfa.core] WARNING: ErfaWarning: ERFA function "taiutc" yielded 1024 of "dubious year (Note 4)" [astropy._erfa.core]
# adjust for our observation point alt_az = coord.AltAz(obstime=times, location=loc) moon_aa = moon.transform_to(alt_az) sun_aa = sun.transform_to(alt_az)
WARNING: Tried to get polar motions for times after IERS data is valid. Defaulting to polar motion from the 50-yr mean for those. This may affect precision at the 10s of arcsec level [astropy.coordinates.builtin_frames.utils] WARNING: ErfaWarning: ERFA function "utcut1" yielded 1024 of "dubious year (Note 3)" [astropy._erfa.core] WARNING: ErfaWarning: ERFA function "apio13" yielded 1024 of "dubious year (Note 2)" [astropy._erfa.core] WARNING: ErfaWarning: ERFA function "utctai" yielded 1024 of "dubious year (Note 3)" [astropy._erfa.core] WARNING: ErfaWarning: ERFA function "taiutc" yielded 1024 of "dubious year (Note 4)" [astropy._erfa.core]
# plot the course through the sky of moon and sun # they are very similar! plt.plot(moon_aa.az, moon_aa.alt, marker='None', linestyle='-', label='moon') plt.plot(sun_aa.az, sun_aa.alt, marker='None', linestyle='-', label='sun') plt.xlabel('Azimuth [deg]') plt.ylabel('Altitude [deg]') plt.legend(); # the key to eclipse calculation is the separation of the two bodies in the sky sun_moon_sep = moon_aa.separation(sun_aa)
WARNING: Tried to get polar motions for times after IERS data is valid. Defaulting to polar motion from the 50-yr mean for those. This may affect precision at the 10s of arcsec level [astropy.coordinates.builtin_frames.utils] WARNING: ErfaWarning: ERFA function "apio13" yielded 1024 of "dubious year (Note 2)" [astropy._erfa.core]
# plot the angular separation as a function of time fig, ax = plt.subplots(1, 1, figsize=(10, 6)) ax.plot(dt, sun_moon_sep.to(u.arcmin), marker='None') offset = 350 ax.xaxis_date(tz=tz) ax.set_xlim(dt[sun_moon_sep.argmin() - offset], dt[sun_moon_sep.argmin() + offset]) ax.xaxis.set_major_locator(HourLocator()) ax.xaxis.set_major_formatter(DateFormatter('%H:%M')) ax.xaxis.set_minor_locator(MinuteLocator(byminute=np.arange(15, 60, 15))) ax.xaxis.set_minor_formatter(DateFormatter(':%M')) ax.set_ylim(0, 30) ax.set_xlabel('Local time [{0:%Z}]'.format(dt)) ax.set_ylabel('Sun–Moon separation [arcmin]'); # Animate the eclipse for the time range # centered on the sun # approximate shapes of sun and moon by 1/2-degree ellipses i0 = sun_moon_sep.argmin() - 128 i1 = sun_moon_sep.argmin() + 128 fig, ax = plt.subplots(1, 1, figsize=(6, 6)) x = np.arange(0, 2*np.pi, 0.01) moon_pa = mpl.patches.Ellipse((0,0), width=0.5, height=0.5, color='#666666', zorder=10) sun_pa = mpl.patches.Ellipse((0,0), width=0.5, height=0.5, color='orange', zorder=1) ax.add_patch(moon_pa) ax.add_patch(sun_pa) ax.set_xlabel('Azimuth [deg]') ax.set_ylabel('Altitude [deg]') #ax.set_aspect(1.6) fig.tight_layout() def animate(i): moon_pa.center = [moon_aa.az[i].degree, moon_aa.alt[i].degree] sun_pa.center = [sun_aa.az[i].degree, sun_aa.alt[i].degree] moon_pa.height = 0.5 * np.cos(moon_aa.alt[i]) sun_pa.height = 0.5 * np.cos(sun_aa.alt[i]) az_lim = (sun_aa.az[i].to(u.degree).value - 2, sun_aa.az[i].to(u.degree).value + 2) alt_lim = (sun_aa.alt[i].to(u.degree).value - 2, sun_aa.alt[i].to(u.degree).value + 2) ax.set_xlim(az_lim) ax.set_ylim(alt_lim) return moon_pa, sun_pa def init(): return animate(i0) ani = animation.FuncAnimation(fig, animate, np.arange(i0, i1), init_func=init, interval=25, blit=True, repeat=False) HTML(ani.to_html5_video())
# Animation with frame fixed in the sky i0 = sun_moon_sep.argmin() - 128 i1 = sun_moon_sep.argmin() + 128 fig2, ax = plt.subplots(1, 1, figsize=(12, 4.8)) x = np.arange(0, 2*np.pi, 0.01) moon_pa = mpl.patches.Ellipse((0,0), width=0.5, height=0.5, color='#666666', zorder=10) sun_pa = mpl.patches.Ellipse((0,0), width=0.5, height=0.5, color='#fec44f', zorder=1) ax.add_patch(moon_pa) ax.add_patch(sun_pa) ax.set_xlabel('Azimuth [deg]') ax.set_ylabel('Altitude [deg]') ax.set_xlim(180, 255) ax.set_ylim(35, 65) fig2.tight_layout() def animate2(i): moon_pa.center = [moon_aa.az[i].degree, moon_aa.alt[i].degree] sun_pa.center = [sun_aa.az[i].degree, sun_aa.alt[i].degree] moon_pa.height = 0.5 * np.cos(moon_aa.alt[i]) sun_pa.height = 0.5 * np.cos(sun_aa.alt[i]) return moon_pa, sun_pa def init2(): return animate2(i0) ani2 = animation.FuncAnimation(fig2, animate2, np.arange(i0, i1), init_func=init2, interval=25, blit=True) ani2
<matplotlib.animation.FuncAnimation at 0x7f926888ec50> HTML(ani2.to_html5_video())

Step 3. Summary stats

# maximum obscuration time max_obsc = dt[sun_moon_sep.argmin()] print(max_obsc.strftime("%Y-%m-%d %H:%M:%S"))
2024-04-08 14:07:37

look up subtended arcs of sun and moon

• sun: The sun subtends an angle of approximately 0.52° (31 arc-minutes)
• moon: The moon subtends an angle of approximately 0.54° (32 arc-minutes)

https://www.mathopenref.com/subtend.html

# maximum obscuration percent # area of the lune # https://rechneronline.de/pi/lune.php r_sun = 31 * 60 / 2 r_moon = 32 * 60 / 2 s = sun_moon_sep.arcsec ld = lunedelta=0.25*math.sqrt((r_sun+r_moon+s)*(r_moon+s-r_sun)*(s+r_sun-r_moon)*(r_sun+r_moon-s)) lune_area=2*ld + r_sun*r_sun*(math.acos(((r_moon*r_moon)-(r_sun*r_sun)-(s*s))/(2*r_sun*s))) - r_moon*r_moon*(math.acos(((r_moon*r_moon)+(s*s)-(r_sun*r_sun))/(2*r_moon*s))) # Calculate percentage of sun's disc eclipsed using lune area and sun size percent_eclipse=(1-(lune_area/(math.pi*r_sun*r_sun)))*100 print(percent_eclipse)
90.39218524334261
# first and last contact times # look for angular separation around 31.5 arcmin for ix in range(len(dt)): aax = sun_moon_sep[ix].arcmin if abs(aax - 31.5) < 0.07: dtx = dt[ix] print(aax,dtx.strftime("%Y-%m-%d %H:%M:%S"))
31.497463254119648 2024-04-08 12:54:39
31.474618929868615 2024-04-08 15:18:56

Results

US Naval Observatory posted times have been converted from UTC to US/Central.

value USNO calculations in this worksheet
First contact time 12:51:13.9 12:54:39
Maximum obscuration time 14:07:14.4 14:07:37
Maximum obscuration 92.4% 90.4%
Last contact time 15:21:29.0 15:18:56