<h1>Euclidean Approach</h1>
<p>I still am not quite sure why the euclidean approach is off.</p>
<h1>Spherical Approach</h1>
<p>The correction to this approach is thanks to the algorithm collection website posted by Greg in the comments.</p>
<p>The error in the spherical approach happens in my definition of the hour angle <span class="math-container">$h = \theta_G + \lambda_T - \alpha$</span>. The correct hour angle definition is given by <span class="math-container">$h^{\prime} = \lambda_T - \alpha$</span>.</p>

CC BY-SA 4.0
2024-02-29T20:58:30.690Z
2024-02-29T20:58:30.690Z
Hunter
<h1>Problem</h1>
<p>I would like to plot the small circle defining the visible portion of Earth for an observer at position <span class="math-container">$P$</span> 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.</p>
<p>I believe I have the math solved but I am getting errors in my plot and am not sure where the bug is at.</p>
<h2>Results</h2>
<p><a href="https://i.sstatic.net/oB21m.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/oB21m.png" alt="my result" /></a></p>
<h2>Expected</h2>
<p><a href="https://i.sstatic.net/eBaRW.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/eBaRW.png" alt="expected results" /></a></p>
<h2>Discrepancies</h2>
<h3>Euclidean case</h3>
<p>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?</p>
<p>Potential fixes:</p>
<ul>
<li>Maybe there is a difference with how I've defined my ECEF frame and how it is defined in skyfield</li>
<li>Maybe I'm not accounting for the codomains of numpy trig functions?</li>
<li>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 <a href="https://en.wikipedia.org/wiki/Antipodal_point" rel="nofollow noreferrer">antipodes</a> just a consequence of the algorithm assuming a spherical earth?</li>
<li>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?</li>
</ul>
<h3>Spherical case</h3>
<p>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.</p>
<p><strong>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.</strong></p>
<h1>The Code</h1>
<pre><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 = (&quot;Sun&quot;, &quot;Moon&quot;, &quot;Mercury Barycenter&quot;, &quot;Venus Barycenter&quot;, &quot;Mars Barycenter&quot;, &quot;Jupiter Barycenter&quot;, &quot;Saturn Barycenter&quot;, &quot;Uranus Barycenter&quot;, &quot;Neptune Barycenter&quot;, &quot;Pluto Barycenter&quot;)

def visible_earth_boundary_euclidean(body: str, time: dt.datetime, samples: int = 100) -&gt; 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) -&gt; 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) -&gt; 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) -&gt; 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 = { &quot;projection&quot;: transform_type })
   ax.stock_img()
   ax.grid()
   ax.set_xlabel('Longitude')
   ax.set_ylabel('Latitude')
   plt.title(f&quot;Solar terminator for {time.strftime('%d %b %Y, %I:%M%p %Z')}&quot;)

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

   plt.show()

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

   # Create an aware python datetime object for the observer's location
   date_string = f&quot;{date} {time}&quot;
   date_format = &quot;%Y-%m-%d %H:%M:%S&quot;
   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(&quot;Sun&quot;, aware_observer_time, samples = 1000)
   spherical = visible_earth_boundary_spherical(&quot;Sun&quot;, aware_observer_time, houston_longitude, houston_latitude, samples = 1000)
   antipode_longitudes = get_antipode_longitudes(&quot;Sun&quot;, aware_observer_time)

   # Plot the results
   visualize_results(aware_observer_time, houston_longitude, houston_latitude, antipode_longitudes, euclidean, spherical, &quot;Sun&quot;)
</code></pre>
<h1>The Math</h1>
<h2>Spherical Geometry Approach</h2>
<p>For an observer on Earth at <span class="math-container">$(\text{longitude } \lambda_O, \text{latitude } \phi_O)$</span>, the altitude/elevation <span class="math-container">$a$</span> above the horizon of a celestial object with hour angle <span class="math-container">$h$</span>, declination <span class="math-container">$\delta$</span>, and right ascension <span class="math-container">$\alpha$</span> is given by</p>
<p><span class="math-container">$$
\sin(a) = \sin(\phi_O)\sin(\delta) + \cos(\phi_O)\cos(\delta)\cos(h)
$$</span></p>
<p>We wish to find all latitude <span class="math-container">$\phi_O$</span> such that <span class="math-container">$a = 0$</span> (the geometric center of the celestial object is at the mathematical horizon). Note that <span class="math-container">$h = \theta_G + \lambda_O - \alpha$</span>, where <span class="math-container">$\theta_G$</span> is the greenwich sidereal time. Therefore, the latitude of the terminator <span class="math-container">$\phi_T$</span> as a function of longitude <span class="math-container">$\lambda_T$</span> is given by</p>
<p><span class="math-container">$$
\phi_T = \arctan\left[\frac{-\cos(\theta_G + \lambda_T - \alpha)}{\tan(\delta)}\right], \lambda_T \in [-180^{\circ}, 180^{\circ}]
$$</span></p>
<h2>Euclidean Geometry Approach</h2>
<p><a href="https://i.sstatic.net/2Z8kL.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/2Z8kL.jpg" alt="enter image description here" /></a></p>
<h3>Definitions</h3>
<p><span class="math-container">$O$</span> : the center of the Earth; used as the origin for both reference frames</p>
<p><span class="math-container">$P$</span> : the celestial interest</p>
<p><span class="math-container">$\vec{r}_{P/O}$</span> : the position of <span class="math-container">$P$</span> relative to <span class="math-container">$O$</span></p>
<p><span class="math-container">$S$</span> : the <a href="https://en.wikipedia.org/wiki/Subsolar_point" rel="nofollow noreferrer">subpoint</a> of <span class="math-container">$P$</span></p>
<p><span class="math-container">$(\lambda_S, \phi_S)$</span> : the longitude and latitude of <span class="math-container">$S$</span>, respectively</p>
<p><span class="math-container">$A$</span> : the center of the disk which forms the base of spherical cap</p>
<p><span class="math-container">$B, C$</span> : points on the sphere where rays coming from <span class="math-container">$P$</span> are tangent to the sphere</p>
<p><span class="math-container">$D$</span> : an arbitrary point on the disk</p>
<p><span class="math-container">$R_E$</span> : radius of the earth</p>
<p><span class="math-container">$\theta_C$</span> : the cone half-angle; the polar angle of <span class="math-container">$D$</span>; the angle between positive <span class="math-container">$\hat{z}_C$</span> and <span class="math-container">$\vec{OD}$</span></p>
<p><span class="math-container">$(\psi_D, \gamma_D)$</span> : the longitude and latitude of <span class="math-container">$D$</span>, respectively</p>
<h3>Reference Frames</h3>
<ol>
<li>ECEF: <span class="math-container">$F_E$</span></li>
<li>Cone frame: <span class="math-container">$F_C$</span>.
<ul>
<li><span class="math-container">$\hat{z}_C$</span> points to the center of the small circle</li>
</ul>
</li>
</ol>
<p><span class="math-container">$F_C$</span> is achieved by the following rotations:</p>
<ol>
<li>A rotation about <span class="math-container">$\hat{z}_E$</span> by the longitude <span class="math-container">$\lambda$</span> of the subsolar point <span class="math-container">$S$</span></li>
<li>A rotation about <span class="math-container">$\hat{y}_E$</span> by the colatitude <span class="math-container">$\phi^{'}$</span> of the subsolar point <span class="math-container">$S$</span></li>
</ol>
<p>That is,</p>
<p><span class="math-container">$$
\begin{align}
T_{E-&gt;C} &amp;= R_y(\phi^{'}_S)R_z(\lambda_S)\\
         &amp;= R_y(90^{\circ} - \phi_s)R_z(\lambda_S)\\
         &amp;=
\begin{bmatrix}
    \sin(\phi_S) &amp; 0 &amp; \cos(\phi_S) \\
    0            &amp; 1 &amp; 0            \\
    -\cos(\phi_S)&amp; 0 &amp; \sin(\phi_S)
\end{bmatrix}
\begin{bmatrix}
    \cos(\lambda_S) &amp; -\sin(\lambda_S) &amp; 0 \\
    \sin(\lambda_S) &amp;  \cos(\lambda_S) &amp; 0 \\
    0               &amp; 0                &amp; 1
\end{bmatrix}
\end{align}
$$</span></p>
<p>and therefore,</p>
<p><span class="math-container">$$
T_{C-&gt;E} = T_{E-&gt;C}^{T}
$$</span></p>
<h3>Solution</h3>
<p>The coordinates <span class="math-container">$^C\vec{x}_D$</span> of an arbitrary point <span class="math-container">$D$</span> on the disk, in the cone frame <span class="math-container">$F_C$</span> are</p>
<p><span class="math-container">$$
\begin{align}
\text{Spherical} &amp;= (R_E, \theta_C, \psi_D)\\
\text{Cartesian} &amp;= 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} &amp;= (\psi, \gamma)_D = (\psi, 90^{\circ} - \theta_C)
\end{align}
$$</span></p>
<p>The coordinates <span class="math-container">$^E\vec{x}_D$</span> of an arbitrary point <span class="math-container">$D$</span> on the disk, in the ECEF frame <span class="math-container">$F_E$</span> are</p>
<p><span class="math-container">$$
\begin{align}
\text{Spherical} &amp;= (R_E, \phi^{'}_D, \lambda_D)\\
\text{Cartesian} &amp;= 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} &amp;= (\lambda, \phi^{'})_D\\
&amp;= (\lambda, 90^{\circ} - \phi)_D\\
&amp;= (\text{arctan2}(y. x), 90^{\circ} - \text{arccos}(z/R_E))
\end{align}
$$</span></p>
<p><strong>I am plotting the geographic coordinates in E.</strong></p>
<h1>Edit: Ruling out possibilities</h1>
<ul>
<li>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.</li>
</ul>
<p>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.</p>

CC BY-SA 4.0
2024-02-22T17:02:36.473Z
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.
Greg Miller
CC BY-SA 4.0
2024-02-27T02:54:42.603Z
This may be of some use, it draws the terminator but using a different algorithm: https://www.celestialprogramming.com/snippets/terminator.html
Greg Miller
CC BY-SA 4.0
2024-02-27T05:23:20.577Z
@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: https://en.wikipedia.org/wiki/Hour_angle#Relation_with_right_ascension
Hunter
CC BY-SA 4.0
2024-02-29T03:17:04.797Z
@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
<observational-astronomy><the-sun><earth><python><skyfield>
Attempting to plot solar terminator using Python and Skyfield library
Hunter
# Attempting to plot solar terminator using Python and Skyfield library

(Asked by: username_0 on 2024-02-20)

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

ECEF: $F_E$
Cone frame: $F_C$.

$\hat{z}_C$ points to the center of the small circle



$F_C$ is achieved by the following rotations:

A rotation about $\hat{z}_E$ by the longitude $\lambda$ of the subsolar point $S$
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.

**Tags:** observational-astronomy, the-sun, earth, python, skyfield

> **Comments:**
> **username_1**: 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.  
> **username_1**: This may be of some use, it draws the terminator but using a different algorithm: https://www.celestialprogramming.com/snippets/terminator.html  
> **username_0**: @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: https://en.wikipedia.org/wiki/Hour_angle#Relation_with_right_ascension  
> **username_2**: @username_0 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.  

---

## 1 answers

### Answer by username_0 on 2024-02-29. Score: 2 | Accepted answer

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


---

