I'm trying to use astropy to plot a simulation of the solar eclipse that we had in Exmouth 20 April 2023 11:31 Local time. I'm a programmer not a physist! So after much messing around I tried this:
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
from matplotlib.animation import FuncAnimation
from astropy.coordinates import EarthLocation, AltAz, get_body,solar_system_ephemeris
from astropy.time import Time
from astropy import units as u
import os
from matplotlib.transforms import Bbox
#SERRURIER ISLAND 21.6087° S, 114.6802° E
# Define the location of Exmouth
with solar_system_ephemeris.set('jpl'):
location = EarthLocation.from_geodetic(114.6802*u.deg,-21.6087*u.deg,0*u.m)
# Define the time range for the eclipse
start_time = Time('2023-04-20 02:00:00') # Start time (UTC)
end_time = Time('2023-04-20 06:00:00') # End time (UTC)
time_resolution = 60 * u.second # Time resolution (1 second)
times = start_time + np.arange(0, (end_time - start_time).to(u.second).value + 1, time_resolution.value) * u.second
# Convert times to local time in Exmouth (UTC+8)
local_times = times + 8*u.hour
# Get the positions of the sun and moon
alt_az = AltAz(obstime=times, location=location)
sun_altaz = get_body('sun',times,ephemeris='de440s').transform_to(alt_az)
moon_altaz = get_body('moon',times,ephemeris='de440s').transform_to(alt_az)
sun_moon_sep = moon_altaz.separation(sun_altaz)
i0 = sun_moon_sep.argmin() - 128
i1 = sun_moon_sep.argmin() + 128
# Calculate the angular diameters
sun_diameter = 1_391_000 * u.km
moon_diameter = 3_474 * u.km
distance_to_sun = 149_600_000 * u.km
distance_to_moon = 384_400 * u.km
sun_angular_radius = ((sun_diameter / sun_altaz.distance)*u.rad).to(u.deg).value
moon_angular_radius = ((moon_diameter / moon_altaz.distance)*u.rad).to(u.deg).value
# Calculate the angular sizes of the sun and moon
#sun_angular_radius = np.arctan2(696340*u.km, get_body('sun',times).distance).to(u.deg)
#moon_angular_radius = np.arctan2(1737.4*u.km, get_body('moon',times).distance).to(u.deg)
# Create a directory to save the frames
output_dir = '/home/mor582/eclipse_frames'
os.makedirs(output_dir, exist_ok=True)
# Set up the plot
fig, ax = plt.subplots(figsize=(19.2,10.8), dpi=100)
plt.gca().set_position([0, 0, 1, 1])
ax.set_xlim(-8, 8)
#ax.set_ylim(-8*0.5625, 8*0.5625)
ax.set_aspect('equal')
ax.set_xlabel('Degrees')
ax.set_ylabel('Degrees')
ax.set_axis_off()
xlim = ax.get_xlim()
ylim = ax.get_ylim()
# Plot the sun and moon
moon_circle = mpl.patches.Ellipse((0,0), width=moon_angular_radius[0], height=moon_angular_radius[0],
color='#666666', zorder=10)
sun_circle = mpl.patches.Ellipse((0,0), width=sun_angular_radius[0], height=sun_angular_radius[0],
color='orange', zorder=1)
# sun_circle = plt.Circle((0.0, 0.0), sun_angular_radius.value, color='yellow', ec='black', lw=2)
# # Set the new position for the time text
#moon_circle = plt.Circle((0.0, 0.0), moon_angular_radius.value, color='gray', ec='black', lw=2)
ax.add_artist(sun_circle)
ax.add_artist(moon_circle)
# Add text annotation for the local time with a semi-transparent background
time_text = ax.text(0, 0, '2023-04-20 11:00:00', fontsize=12, ha='center', color='black',
bbox=dict(facecolor='white', alpha=0.6, edgecolor='none'),fontweight='bold')
time_text_bbox = time_text.get_window_extent(renderer=fig.canvas.get_renderer(),dpi=100)
bbox_text = Bbox(ax.transData.inverted().transform(time_text_bbox))
time_text_x = xlim[1] - 0.1 -bbox_text.width/2
time_text_y = ylim[1] -0.01 -bbox_text.height
# Set the new position for the time text
time_text.set_position((time_text_x, time_text_y))
sun_bbox = sun_circle.get_window_extent(renderer=fig.canvas.get_renderer())
bbox_sun = Bbox(ax.transData.inverted().transform(sun_bbox))
sun_x = xlim[1] - bbox_sun.width/2
sun_y = time_text_y -ax.margins()[1] -bbox_sun.height
#sun_circle.center =(sun_x,sun_y)
# Set the background color to be transparent
fig.patch.set_alpha(0)
ax.patch.set_alpha(0)
ax.axis('off')
def wrap_north(bearing):
if bearing>180:
bearing = bearing -360
return bearing
# Function to update the position of the moon and the sun's transparency in the plot
def update(frame):
# Update the moon's position
moon_circle.center = (wrap_north(moon_altaz.az[frame].deg),moon_altaz.alt[frame].deg)
sun_circle.center = (wrap_north(sun_altaz.az[frame].deg),sun_altaz.alt[frame].deg)
az_lim = (wrap_north(sun_altaz.az[frame].deg) - 2,
wrap_north(sun_altaz.az[frame].deg) + 2)
#alt_lim = (sun_aa.alt[i].to(u.degree).value - 2,
# sun_aa.alt[i].to(u.degree).value + 2)
alt_lim = (sun_altaz.alt[frame].deg - 2,
sun_altaz.alt[frame].deg + 2)
ax.set_xlim(az_lim)
ax.set_ylim(alt_lim)
time_text_bbox = time_text.get_window_extent(renderer=fig.canvas.get_renderer(),dpi=100)
bbox_text = Bbox(ax.transData.inverted().transform(time_text_bbox))
time_text_x = az_lim[1]-bbox_text.width/2
time_text_y = alt_lim[1] -bbox_text.height *2
#moon_circle.height = moon_angular_radius * np.cos(moon_altaz.alt[frame])
#sun_circle.height = sun_angular_radius * np.cos(sun_altaz.alt[frame])
# Calculate the alpha value for the sun's transparency based on its altitude
#alpha = max(0, min(1, (sun_altitude - sun_altaz.alt.min().deg) / (sun_altaz.alt.max().deg - sun_altaz.alt.min().deg)))
#sun_circle.set_alpha(alpha)
# Update the local time text
local_time_str = local_times[frame].iso[:-4] # Format time string without seconds
time_text.set_text(local_time_str)
time_text.set_position([time_text_x,time_text_y])
# Position the time text in the top-right corner
# Save the frame as a PNG file with a transparent background
file_time=local_times[frame].iso.replace('-', '').replace(':', '').replace(' ', 'T').split('.')[0]
hour_dir =f"{output_dir}/{file_time[:-4]}"
# if not os.path.exists(hour_dir):
# os.mkdir(hour_dir)
#/{file_time[:-4]}
plt.savefig(f"{output_dir}/frame_{file_time}_{frame:04d}.png", transparent=True)
return sun_circle, moon_circle, time_text
# Create animation
ani = FuncAnimation(fig, update, frames=len(times), blit=True, interval=1000)
# Save each frame as an image
for i in range(len(times)):
update(i)
plt.close()
Here is a link to what I'm trying to recreate
https://www.timeanddate.com/eclipse/in/australia/exmouth?iso=20230420
There seems to be a difference in the Sun's position between Astropy and this rendering.
If I use sunpy and the sun.eclipse_amount(observer, moon_radius='minimum') it gives the correct times for the start and end of the eclipse.
