# 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.

## 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

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
earth, body_ephemeride = ephemeris['earth'], ephemeris[body]
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

# 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)
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)]

# 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
skyfield_time = timescale.from_datetime(time)
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?)
greenwich_apparent_sidereal_time = skyfield_time.gast

terminator_longitudes = np.linspace(-np.pi, np.pi, samples)
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]:
skyfield_time = timescale.from_datetime(time)
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.

• The most comon definition of sunset is when the geometric center is .833deg below the horizon. This accounts for the size of the sun and refraction. I think the errors your seeing are too big for that to be the source of the problem, but keep that in mind for later. Commented Feb 22 at 17:02
• This may be of some use, it draws the terminator but using a different algorithm: celestialprogramming.com/snippets/terminator.html Commented Feb 27 at 2:54
• @GregMiller This gave me the results I'm looking for! Ok the difference between the website and my code is they define the hour angle as $h = \lambda_T - \lambda_S$ where $\lambda_S$ is the longitude of the sun's subpoint. Why was I wrong? My hour angle definition $h = \theta_G + \lambda_T - \alpha$ is directly from the wikipedia page: en.wikipedia.org/wiki/Hour_angle#Relation_with_right_ascension Commented Feb 27 at 5:23
• @Hunter have you got a solution? If so you are welcome to post it as an answer! The bounty ends in a few days, and it would be a shame for it to go to waste. It's always okay to answer your own question in Stack Exchange if you find the answer.
– uhoh
Commented Feb 29 at 3:17

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$$.