Contact
CoCalc Logo Icon
StoreFeaturesDocsShareSupport News AboutSign UpSign In
| Download

2024 Solar Eclipse demo

Views: 306
%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.

References

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[0]))
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[544].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.

Source: Solar Eclipse Computer - Solar Eclipse of 2024 April 08 - Palatine, IL

valueUSNOcalculations in this worksheet
First contact time12:51:13.912:54:39
Maximum obscuration time14:07:14.414:07:37
Maximum obscuration92.4%90.4%
Last contact time15:21:29.015:18:56