Attempting to plot solar terminator using Python and Skyfield library

Question

Problem

I would like to plot the small circle defining the visible portion of Earth for an observer at position $P$ away from the Earth, on a 2D map. In the case of the Sun, this would be the solar terminator. Without loss of generality, we can focus strictly on the Sun case.

I believe I have the math solved but I am getting errors in my plot and am not sure where the bug is at.

Results

Expected

Discrepancies

Euclidean case

The euclidean approach seems to be pretty fine, but it is consistently not centered around the poles. How can I get it centered around the poles?

Potential fixes:

• Maybe there is a difference with how I've defined my ECEF frame and how it is defined in skyfield
• Maybe I'm not accounting for the codomains of numpy trig functions?
• I'm using the geodetic latitude and longitude of the subsolar point as inputs to an algorithm that assumes a spherical Earth. Is the plot being off-center of the antipodes just a consequence of the algorithm assuming a spherical earth?
• Mathematically speaking, I have calculated a small circle with poles that are supposed to align with the subsolar point, but they don't. If this difference in poles is simply due to my algorithm assuming a spherical earth, then I suppose I need to check the error in the poles and rotate the small circle to align with the correct poles?

Spherical case

The spherical approach seems to be consistently flipped and out of phase by about 180 degrees (or maybe it's just out of phase by 180 degrees?). The spherical approach also keeps giving latitudes outside of the [-90, 90] range and therefore the spherical approach often can not be plotted.

Where am I going wrong with either one of these approaches? I only need a correction to just ONE of the approaches, I do not need both approaches.

The Code

import skyfield
from skyfield import positionlib, api
import numpy as np
from numpy import cos, sin, tan
import datetime as dt
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from pytz import timezone

MEAN_EARTH_RADIUS_KM = 6371
SOUTHERN_MOST_LATITUDE = -90
NORTHERN_MOST_LATITUDE = 90
CELESTIAL_INTERESTS = ("Sun", "Moon", "Mercury Barycenter", "Venus Barycenter", "Mars Barycenter", "Jupiter Barycenter", "Saturn Barycenter", "Uranus Barycenter", "Neptune Barycenter", "Pluto Barycenter")

def visible_earth_boundary_euclidean(body: str, time: dt.datetime, samples: int = 100) -> tuple[tuple[float, float]]:
   import numpy as np
   from numpy import cos, sin
   from scipy.spatial.transform import Rotation as R

   # Get the body's geocentric, astrometric position
   ephemeris = skyfield.api.load('de421.bsp')
   earth, body_ephemeride = ephemeris['earth'], ephemeris[body]
   timescale = skyfield.api.load.timescale()
   skyfield_time = timescale.from_datetime(time)
   astrometric_position = earth.at(skyfield_time).observe(body_ephemeride)

   # The cone aperature = 2 * cone_half_angle; https://en.wikipedia.org/wiki/Cone#Further_terminology
   cone_half_angle = np.arccos(MEAN_EARTH_RADIUS_KM / astrometric_position.distance().km)

   # The rotation matrix from cone frame to earth frame
   # The code below is the transpose found here: https://shorturl.at/iCWZ7
   position_subpoint = skyfield.api.wgs84.subpoint_of(astrometric_position)
   lon = position_subpoint.longitude.radians
   lat = position_subpoint.latitude.radians
   row1 = [ cos(lon)*sin(lat), sin(lon), -cos(lon)*cos(lat)]
   row2 = [-sin(lon)*sin(lat), cos(lon),  sin(lon)*cos(lat)]
   row3 = [ cos(lat),          0       ,  sin(lat)         ]
   cone_frame_to_earth_frame_rotation =  R.from_matrix(np.array([row1, row2, row3]))

   # Psi is the longitude of the position along the disk
   # That is, it is the angle between the cone's positive x-axis and the projeciton of the disk point position vector on the xy-plane
   psi = np.linspace(0, 2 * np.pi, samples)

   # Make the z value match the length of the x and y values since np.array(list) will expect list to be a multidimensional array
   # In other words, np.array([ [x1, x2], [y1, y2], z1 ]) will not work, but np.array([ [x1, x2], [y1, y2], [z1, z2] ]) will.
   cartesian_coordinates_in_cone_frame = [sin(cone_half_angle)*cos(psi), sin(cone_half_angle)*sin(psi), [cos(cone_half_angle)] * len(psi)]

   
   position_along_small_circle_in_cone_frame = MEAN_EARTH_RADIUS_KM * np.array(cartesian_coordinates_in_cone_frame)

   # Since np.array(cartesian_coordinates_in_cone_frame) has shape [ [x1, x2, ..., xn], [y1, y2, ..., yn], [z1, z2, ..., zn] ]
   # and Rotation.apply() expects an input with the form [ [x1, y1, z1], [x2, y2, z2], ... [xn, yn, zn] ]
   # We take the transpose to achieve this
   position_along_small_circle_in_earth_frame = cone_frame_to_earth_frame_rotation.apply(position_along_small_circle_in_cone_frame.T)

   # In order to index the position's x, y, and z values in an intutive way (e.g. r[0] for x) we take the transpose
   r = position_along_small_circle_in_earth_frame.T

   # Compute the geographic coordiantes from the cartesian coordinates
   small_circle_latlons_in_earth_frame = np.degrees(np.array((np.arctan2(r[1], r[0]), (np.pi/2.0) - np.arccos(r[2]/MEAN_EARTH_RADIUS_KM))))

   return small_circle_latlons_in_earth_frame

def visible_earth_boundary_spherical(body: str, time: dt.datetime, longitude: float, latitude: float, samples: int = 100) -> np.array:
   # Get astrometric position of body relative to earth observer
   timescale = skyfield.api.load.timescale()
   skyfield_time = timescale.from_datetime(time)
   ephemeris = skyfield.api.load('de421.bsp')
   earth, body_ephemeride = ephemeris['earth'], ephemeris[body]
   observer_location = earth + skyfield.api.wgs84.latlon(latitude, longitude)
   astrometric_position = observer_location.at(skyfield_time).observe(body_ephemeride)

   # Relative to ICRS (i think?)
   right_ascension, declination, distance = astrometric_position.radec()
   hour_angle, _, _ = astrometric_position.hadec()
   greenwich_apparent_sidereal_time = skyfield_time.gast

   terminator_longitudes = np.linspace(-np.pi, np.pi, samples)
   terminator_latitudes = np.arctan2(-cos(greenwich_apparent_sidereal_time + terminator_longitudes - right_ascension.radians), tan(declination.radians))
   terminator_coordinates = np.degrees(np.array([terminator_longitudes, terminator_latitudes]))
   return terminator_coordinates

def get_antipode_longitudes(body: str, time: dt.datetime) -> tuple[float, float]:
   timescale = skyfield.api.load.timescale()
   skyfield_time = timescale.from_datetime(time)
   ephemeris = skyfield.api.load('de421.bsp')
   earth, body_ephemeride = ephemeris['earth'], ephemeris[body]
   astrometric_position = earth.at(skyfield_time).observe(body_ephemeride)
   subpoint = skyfield.api.wgs84.subpoint_of(astrometric_position)
   subpoint_longitude = subpoint.longitude.degrees
   antipode_longitude = (subpoint_longitude % 360) - 180
   return (subpoint_longitude, antipode_longitude)

def visualize_results(time: dt.datetime, longitude: float, latitude: float, antipode_longitudes: tuple[float, float], euclidean: np.array, spherical: np.array, body: str) -> None:
   subpoint_coordinates = ((antipode_longitudes[0], antipode_longitudes[0]), (SOUTHERN_MOST_LATITUDE, NORTHERN_MOST_LATITUDE))
   antipode_coordinates = ((antipode_longitudes[1], antipode_longitudes[1]), (SOUTHERN_MOST_LATITUDE, NORTHERN_MOST_LATITUDE))

   transform_type = ccrs.PlateCarree()
   fig, ax = plt.subplots(subplot_kw = { "projection": transform_type })
   ax.stock_img()
   ax.grid()
   ax.set_xlabel('Longitude')
   ax.set_ylabel('Latitude')
   plt.title(f"Solar terminator for {time.strftime('%d %b %Y, %I:%M%p %Z')}")

   ax.plot(longitude, latitude, marker = ".")
   ax.plot(subpoint_coordinates[0], subpoint_coordinates[1], transform = transform_type, color = "green", label = "Solar noon meridian")
   ax.plot(antipode_coordinates[0], antipode_coordinates[1], transform = transform_type, color = "yellow", label = "Solar midnight meridian")
   ax.scatter(euclidean[0], euclidean[1], transform = transform_type, color = "red", marker = ".", s = 1.5, label = "Euclidian approach")
   ax.scatter(spherical[0], spherical[1], transform = transform_type, color = "blue", marker = ".", s = 1.5, label = "Spherical approach")
   plt.legend()

   plt.show()

if __name__ == "__main__":
   # Define the input parameters
   houston_latitude = 29.7604
   houston_longitude = -95.3698
   time = "12:00:00"
   houston_timezone = "America/Chicago"
   date = "2021-07-22"

   # Create an aware python datetime object for the observer's location
   date_string = f"{date} {time}"
   date_format = "%Y-%m-%d %H:%M:%S"
   naive_observer_time = dt.datetime.strptime(date_string, date_format)
   time_zone = timezone(houston_timezone)
   aware_observer_time = time_zone.localize(naive_observer_time)

   # Calculate the sun terminator
   euclidean = visible_earth_boundary_euclidean("Sun", aware_observer_time, samples = 1000)
   spherical = visible_earth_boundary_spherical("Sun", aware_observer_time, houston_longitude, houston_latitude, samples = 1000)
   antipode_longitudes = get_antipode_longitudes("Sun", aware_observer_time)

   # Plot the results
   visualize_results(aware_observer_time, houston_longitude, houston_latitude, antipode_longitudes, euclidean, spherical, "Sun")

The Math

Spherical Geometry Approach

For an observer on Earth at $(\text{longitude } \lambda_O, \text{latitude } \phi_O)$, the altitude/elevation $a$ above the horizon of a celestial object with hour angle $h$, declination $\delta$, and right ascension $\alpha$ is given by

$$ \sin(a) = \sin(\phi_O)\sin(\delta) + \cos(\phi_O)\cos(\delta)\cos(h) $$

We wish to find all latitude $\phi_O$ such that $a = 0$ (the geometric center of the celestial object is at the mathematical horizon). Note that $h = \theta_G + \lambda_O - \alpha$, where $\theta_G$ is the greenwich sidereal time. Therefore, the latitude of the terminator $\phi_T$ as a function of longitude $\lambda_T$ is given by

$$ \phi_T = \arctan\left[\frac{-\cos(\theta_G + \lambda_T - \alpha)}{\tan(\delta)}\right], \lambda_T \in [-180^{\circ}, 180^{\circ}] $$

Euclidean Geometry Approach

Definitions

$O$ : the center of the Earth; used as the origin for both reference frames

$P$ : the celestial interest

$\vec{r}_{P/O}$ : the position of $P$ relative to $O$

$S$ : the subpoint of $P$

$(\lambda_S, \phi_S)$ : the longitude and latitude of $S$, respectively

$A$ : the center of the disk which forms the base of spherical cap

$B, C$ : points on the sphere where rays coming from $P$ are tangent to the sphere

$D$ : an arbitrary point on the disk

$R_E$ : radius of the earth

$\theta_C$ : the cone half-angle; the polar angle of $D$; the angle between positive $\hat{z}_C$ and $\vec{OD}$

$(\psi_D, \gamma_D)$ : the longitude and latitude of $D$, respectively

Reference Frames

1. ECEF: $F_E$
2. Cone frame: $F_C$.
   • $\hat{z}_C$ points to the center of the small circle

$F_C$ is achieved by the following rotations:

1. A rotation about $\hat{z}_E$ by the longitude $\lambda$ of the subsolar point $S$
2. A rotation about $\hat{y}_E$ by the colatitude $\phi^{'}$ of the subsolar point $S$

That is,

$$ \begin{align} T_{E->C} &= R_y(\phi^{'}_S)R_z(\lambda_S)\\ &= R_y(90^{\circ} - \phi_s)R_z(\lambda_S)\\ &= \begin{bmatrix} \sin(\phi_S) & 0 & \cos(\phi_S) \\ 0 & 1 & 0 \\ -\cos(\phi_S)& 0 & \sin(\phi_S) \end{bmatrix} \begin{bmatrix} \cos(\lambda_S) & -\sin(\lambda_S) & 0 \\ \sin(\lambda_S) & \cos(\lambda_S) & 0 \\ 0 & 0 & 1 \end{bmatrix} \end{align} $$

and therefore,

$$ T_{C->E} = T_{E->C}^{T} $$

Solution

The coordinates $^C\vec{x}_D$ of an arbitrary point $D$ on the disk, in the cone frame $F_C$ are

$$ \begin{align} \text{Spherical} &= (R_E, \theta_C, \psi_D)\\ \text{Cartesian} &= R_E \begin{bmatrix} \sin(\theta_C)\cos(\psi_D) \\ \sin(\theta_C)\sin(\psi_D) \\ \cos(\theta_C) \end{bmatrix}\\ \text{Geographic in C} &= (\psi, \gamma)_D = (\psi, 90^{\circ} - \theta_C) \end{align} $$

The coordinates $^E\vec{x}_D$ of an arbitrary point $D$ on the disk, in the ECEF frame $F_E$ are

$$ \begin{align} \text{Spherical} &= (R_E, \phi^{'}_D, \lambda_D)\\ \text{Cartesian} &= R_E \begin{bmatrix} \sin(\phi^{'}_D)\cos(\lambda_D) \\ \sin(\phi^{'}_D)\sin(\lambda_D) \\ \cos(\phi^{'}_D) \end{bmatrix}\\ \text{Geographic in E} &= (\lambda, \phi^{'})_D\\ &= (\lambda, 90^{\circ} - \phi)_D\\ &= (\text{arctan2}(y. x), 90^{\circ} - \text{arccos}(z/R_E)) \end{align} $$

I am plotting the geographic coordinates in E.

Edit: Ruling out possibilities

• Reference frames differences: I have been working on correcting the euclidean case. I believe I have ruled out the possibility that my reference frames are different than the skyfield reference frames. I don't believe I have an error in this area.

For as of right now, I am marking up my difference (in the euclidean case) as a result of assuming a spherical earth, and as a work around I have calculated the angle difference between the disk normal and the subpoint vector, and then rotate the disk such that it's normal vector points in the subpoint vector direction. This seems to give the expected results but I am not sure if this is just a bandaid on a more egregious error somewhere.

Answer 1

Euclidean Approach

I still am not quite sure why the euclidean approach is off.

Spherical Approach

The correction to this approach is thanks to the algorithm collection website posted by Greg in the comments.

The error in the spherical approach happens in my definition of the hour angle $h = \theta_G + \lambda_T - \alpha$. The correct hour angle definition is given by $h^{\prime} = \lambda_T - \alpha$.