<p>Here's a solution that uses only Astropy and a guess for the time of the eclipse.
It doesn't use any Besselian elements, instead it finds the latitude and time that minimizes the separation between the Sun and Moon.</p>
<pre><code>import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize
import astropy.units as u
import astropy.time
import astropy.coordinates
import astropy.constants
import cartopy

# Define the radius of the moon
radius_moon = 1737.4 * u.km

def get_sun_moon(
    location: astropy.coordinates.EarthLocation,
    time: astropy.time.Time,
) -&gt; tuple[astropy.coordinates.SkyCoord, astropy.coordinates.SkyCoord]:
  &quot;&quot;&quot;
  Given a location on the Earth and an observation time,
  get the coordinates of the Sun and Moon in the local
  coordinate system.

  Parameters
  ----------
  location
    A point on the surface of the earth
  time
    An observation time
  &quot;&quot;&quot;
  frame_local = astropy.coordinates.AltAz(
      location=location,
      obstime=time,
  )
  sun = astropy.coordinates.get_body(
      body=&quot;sun&quot;,
      location=location,
      time=time,
      ephemeris=&quot;de440&quot;,
  ).transform_to(frame_local)
  moon = astropy.coordinates.get_body(
      body=&quot;moon&quot;,
      location=location,
      time=time,
      ephemeris=&quot;de440&quot;,
  ).transform_to(frame_local)
  return sun, moon


def sun_moon_separation(
    location: astropy.coordinates.EarthLocation,
    time: astropy.time.Time,
) -&gt; u.Quantity:
  &quot;&quot;&quot;
  Find the angular separation of the Sun and Moon
  at the given Earth location and time.

  Parameters
  ----------
  location
    A point on the surface of the Earth
  time
    An observation time
  &quot;&quot;&quot;
  sun, moon = get_sun_moon(location, time)
  return sun.separation(moon)


def sun_moon_separation_argmin(
    location: astropy.coordinates.EarthLocation,
    time: astropy.time.Time,
) -&gt; tuple[
    astropy.coordinates.EarthLocation,
    astropy.time.Time,
]:
  &quot;&quot;&quot;
  Compute the time and latitude of the closest separation
  between the Sun and Moon using local minimization.

  Parameters
  ----------
  location
    A point on the Earth to use as an initial guess.
    The longitude is kept constant and the latitude is
    allowed to vary.
  time
    An initial guess for the time of closest approach.
  &quot;&quot;&quot;
  def factory(
      lat: float, 
      t: float
  ) -&gt; tuple[
      astropy.coordinates.EarthLocation,
      astropy.time.Time,
  ]:
    loc = astropy.coordinates.EarthLocation.from_geodetic(
        lat=lat * u.deg,
        lon=location.lon,
    )
    t = astropy.time.Time(t, format=&quot;jyear&quot;)
    return loc, t

  def fun(x: tuple[float, float]) -&gt; float:
    lat, t = x
    loc, t = factory(lat, t)
    res = sun_moon_separation(loc, t)
    return res

  lat, t = scipy.optimize.minimize(
      fun=fun,
      x0=[location.lat.to_value(u.deg), time.jyear],
      bounds=[
          (-90, 90),
          (time.jyear - 0.01, time.jyear + 0.01),
      ],
      method=&quot;nelder-mead&quot;,
  ).x

  return factory(lat, t)


def calc_path_center(
    longitude: u.Quantity,
    time_guess: astropy.time.Time,
) -&gt; tuple[
    astropy.coordinates.EarthLocation, 
    astropy.time.Time,
]:
  &quot;&quot;&quot;
  Compute the center of the eclipse path across the Earth

  Parameters
  ----------
  longitude
    A grid of longitudes at which to compute the center of the eclipse
  time_guess
    An initial guess of the time at which the eclipse occurs.
  &quot;&quot;&quot;
  time = astropy.time.Time(
      np.zeros(longitude.shape),
      format=&quot;jyear&quot;,
  )
  latitude = np.empty_like(longitude)

  lat_guess = 0 * u.deg

  for i in np.ndindex(longitude.shape):
    location_i, time_i = sun_moon_separation_argmin(
        location=astropy.coordinates.EarthLocation.from_geodetic(
            lat=lat_guess,
            lon=longitude[i],
        ),
        time=time_guess,
    )
    latitude[i] = lat_guess = location_i.lat
    time[i] = time_guess = time_i

  location = astropy.coordinates.EarthLocation.from_geodetic(
      lat=latitude,
      lon=longitude,
  )
  sun, moon = get_sun_moon(
      location=location,
      time=time,
  )
  separation = sun.separation(moon)
  angular_radius_sun = np.arctan2(
      astropy.constants.R_sun,
      sun.cartesian.norm(),
  )
  angular_radius_moon = np.arctan2(
      radius_moon,
      moon.cartesian.norm(),
  )

  where_eclipse = separation &lt; (angular_radius_moon - angular_radius_sun)
  where_eclipse = where_eclipse &amp; (sun.alt &gt; 0)

  location = location[where_eclipse]
  time = time[where_eclipse]

  return location, time


# Define the start and stop times
# of the search
time_start = astropy.time.Time(&quot;2024-04-08T18:00&quot;)

# Define a grid of longitudes to search
longitude = np.linspace(-180, 180, num=101) * u.deg

# Compute the eclipse path
location_center, time_center = calc_path_center(
    longitude=longitude,
    time_guess=time_start,
)

# Plot the eclipse path
fig, ax = plt.subplots(
    subplot_kw=dict(projection=cartopy.crs.Robinson()),
)
ax.stock_img();
ax.coastlines();
ax.plot(
    location_center.lon,
    location_center.lat,
    transform=cartopy.crs.PlateCarree(),
);
</code></pre>
<p><a href="https://i.sstatic.net/XIJnUAzc.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/XIJnUAzc.png" alt="April 2024 Eclipse" /></a></p>

CC BY-SA 4.0
2024-05-03T02:41:52.117Z
I started with a solution like this, but the accuracy using Astropy was off by about a degree or so, and I ended up giving up on it.
PearsonArtPhoto
CC BY-SA 4.0
2024-05-03T03:09:08.113Z
Astropy should really be just as close as anything else, if it's not that would be a big problem. Were you using the default ephemeris? I would expect the default ephemeris (as opposed to the JPL ephemeris) to be off by a little bit, but not a whole degree.
Roy Smart
CC BY-SA 4.0
2024-05-03T19:05:38.930Z
I tried both JPL and the default, they were both fairly similar, but off. Maybe my code was off, however...
PearsonArtPhoto
CC BY-SA 4.0
2024-05-04T00:46:35.247Z
Not sure what I am doing, but I get `25.613N, 103.832 W time 2024-04-08T18:18:29` as the peak with my effectively identical method, while the peak location at that time according to the catalog was 25.3 N, 104.1 W. This with the de440.
PearsonArtPhoto
CC BY-SA 4.0
2024-05-04T21:54:14.900Z
It seems like the numbers `25.613N, 103.832 W time 2024-04-08T18:18:29` match the values on @GregMiller's site, can you show me what tabulated values you're using?
Roy Smart
CC BY-SA 4.0
2024-05-04T21:56:08.777Z
Also, confusingly, my code right now gives `25.802 N, 103.832 W time 2024-04-08T18:20:02.846` for that longitude, so I'm not sure why our codes are different.
Roy Smart
CC BY-SA 4.0
2024-05-05T22:07:04.223Z
I don't have exactly the same code, just functionally the same (Same key functions, different working goods). I'm comparing it to the values reported in the NASA catalog for the event.
PearsonArtPhoto
CC BY-SA 4.0
2024-05-02T22:50:19.013Z
2024-05-02T22:50:19.013Z
Roy Smart
<p>As you've already found out there aren't many good resources for using Besselian Elements to compute solar eclipses.  While there are some thorough resources, they are a bit terse, and careful study is required in order to figure out exactly which equations to use.  Also Besselian Elements is a bit of a rabbit hole, as there are quite a few different algorithms for computing various phenomenon.</p>
<p>I have had it on my todo list to write a simplified article that can be followed by a typical programmer without much astronomy knowledge.  So far, <a href="https://www.celestialprogramming.com/articles/eclipsePredictions/predictions.html" rel="nofollow noreferrer">I just have a stub</a> which has implementations of many examples from various sources.  An article on <a href="https://www.celestialprogramming.com/eclipsesGeneratingBesselianElements/generatingBesselianPolynomials.html" rel="nofollow noreferrer">computing Besselian Elements</a>. And an <a href="https://www.celestialprogramming.com/apps/SolarEclipseViewer/viewer.html" rel="nofollow noreferrer">example implementation of an eclipse map</a> which implements most of the algorithms.</p>
<p>Below I'll just answer the question you asked, which is computing the central line of an eclipse which is simple enough to reproduce here, and enough to get your feet wet implementing other algorithms.</p>
<p>In the interest of brevity, I'm not going to try to explain what each of the variables represent, as most of them are intermediate values without a simple explanation.  For a more thorough coverage, consult one (or many) of the below resources which are what I used.  The below algorithm is from Meeus.</p>
<ul>
<li>A Manual of Spherical and Practical Astronomy - Chauvenet 1863</li>
<li>The Mathematical Theory of Eclipses According to Chauvenet's Transformation of Bessel's Method - Buchanan 1904</li>
<li>Explanatory Supplement to the Astronomical Ephemeris 1961</li>
<li>Explanatory Supplement to the Astronomical Almanac 1992</li>
<li>Elements of Solar Eclipses - Meeus 1989</li>
<li>Prediction and Analysis of Solar Eclipse Circumstances - Williams 1971</li>
</ul>
<p>You will first need a set of Besselian Elements.  These can be obtained from the <a href="https://eclipse.gsfc.nasa.gov/SEcat5/SEcatalog.html" rel="nofollow noreferrer">Five Millennium Catalog of Solar Eclipses</a>.  I have extracted <a href="https://github.com/gmiller123456/FiveMillenniumCanonOfSolarEclipses-Besselian-Elements" rel="nofollow noreferrer">all of the elements from there</a> and put them in a Github repository to make that process easier.  The elements below are from the <a href="https://eclipse.gsfc.nasa.gov/SEsearch/SEdata.php?Ecl=%2020240408" rel="nofollow noreferrer">2024 Apr 8 eclipse</a>.</p>
<div class="s-table-container"><table class="s-table">
<thead>
<tr>
<th>n</th>
<th>x</th>
<th>y</th>
<th>d</th>
<th>L1</th>
<th>L2</th>
<th>μ</th>
</tr>
</thead>
<tbody>
<tr>
<td>0</td>
<td>-0.3182440</td>
<td>0.2197640</td>
<td>7.5862002</td>
<td>0.5358140</td>
<td>-0.0102720</td>
<td>89.591217</td>
</tr>
<tr>
<td>1</td>
<td>0.5117116</td>
<td>0.2709589</td>
<td>0.0148440</td>
<td>0.0000618</td>
<td>0.0000615</td>
<td>15.004080</td>
</tr>
<tr>
<td>2</td>
<td>0.0000326</td>
<td>-0.0000595</td>
<td>-0.0000020</td>
<td>-0.0000128</td>
<td>-0.0000127</td>
<td>0.000000</td>
</tr>
<tr>
<td>3</td>
<td>-0.0000084</td>
<td>-0.0000047</td>
<td>0.0000000</td>
<td>0.0000000</td>
<td>0.0000000</td>
<td>0.000000</td>
</tr>
</tbody>
</table></div>
<pre><code>tan f1 = 0.0046683        tan f2 = 0.0046450 
&quot;Eclipse Predictions by Fred Espenak, NASA's GSFC&quot;
</code></pre>
<p>In addition to the elements, some auxiliary information is also provided, namely <span class="math-container">$ t0 $</span>, which is the basis time the elements were computed around, any time computed with the elements will be relative to <span class="math-container">$ t0 $</span>, and is in hours. E.g. for 2024 Apr 8, <span class="math-container">$ t0 $</span> is 18.000, or 6pm Terrestrial Dynamic Time (TDT).  There is also <span class="math-container">$ \Delta t $</span>, used to convert TDT to UTC, and is in seconds.</p>
<p>The algorithm below computes the central line for a given time.  So you specify a time <span class="math-container">$ t $</span>, and it produces the latitude <span class="math-container">$ \phi $</span> and longitude <span class="math-container">$ \lambda $</span>, and a few other values.  Generally <span class="math-container">$t$</span> takes the values of -2 to +2, that's two hours before and two hours after <span class="math-container">$t0$</span>.</p>
<p>The first step in most algorithms is to compute the elements for the specified time:
<span class="math-container">$$
\begin{align*}
X &amp;= X_0 +X_1t + X_2t^2 + X_3t^3 \\
Y &amp;= Y_0 + Y_1t + Y_2t^2 + Y_3t^3 \\
d &amp;= d_0 + d_1t + d_2^2 \\
L1 &amp;= L1_0 + L1_1t + L1_2t^2 \\
L2 &amp;= L2_0 + L2_1t + L2_2t^2 \\
\mu &amp;= \mu_0 + \mu_1t \\
\\
X' &amp;= X_1 + 2X_2t + 3X_3t^2 \\
Y' &amp;= Y_1 + 2Y_2t + 3Y_3t^2
\end{align*}
$$</span></p>
<p>So, using the table above, <span class="math-container">$X_0$</span> = -0.3182440, <span class="math-container">$X_1$</span> = 0.5117116, <span class="math-container">$X_2$</span> = 0.0000326, etc. Note that the trigonometric functions are expressed in degrees.</p>
<p><span class="math-container">$$
\begin{align*}
\omega &amp;=\frac{1}{\sqrt{1 - 0.006694385 \cos^2 d}} \\
p &amp;=\mu_1/57.2957795 \\
b &amp;=Y' - pX \sin d \\
c &amp;=X' + pY \sin d \\
y1 &amp;=\omega Y \\
b1 &amp;=\omega \sin d \\
b2 &amp;=0.99664719 \omega \cos d \\
B &amp;= \sqrt{1 -  X^2 - y1^2} \\
\end{align*}
$$</span>
If <span class="math-container">$B$</span> does not exist, there is no central eclipse at given time.
<span class="math-container">$$
\begin{align*}
\sin \phi_1 &amp;=B*b1 + y1*b2 \\
\sin H &amp;= \frac{X}{\cos \phi_1} \\
\cos H &amp;= \frac{B*b2 - y1*b1}{\cos \phi_1} \\
\tan \phi &amp;= 1.00336409 \tan \phi_1 \\
\lambda &amp;= \mu - H - 0.00417807 Δt \\
\end{align*}
$$</span>
<span class="math-container">$\lambda$</span> is the longitude (Positive = West, opposite normal convention), and <span class="math-container">$\phi$</span> is the latitude of the central point for the given time. H should be the number that is in the correct quadrant for the given <span class="math-container">$sin H$</span> and <span class="math-container">$cos H$</span>. Additionally, the width, duration, and magnitude (<span class="math-container">$A$</span>) can be computed without much more work. Note if duration is less than 0 then it is the duration of total eclipse, if greater than 0 it is the duration of annular eclipse.
<span class="math-container">$$
\begin{align*}
L1'&amp;=L1 - B \tan f1 \\
L2'&amp;=L2 - B \tan f2 \\
a &amp;= c - p*B* \cos d \\
n &amp;= \sqrt{a^2 + b^2} \\
duration &amp;= \frac{7200*L2'}{n} \\
\sin h &amp;=\sin d \sin \phi + \cos d \cos \phi \cos H \\
K^2 &amp;= B^2 + \frac{(Xa + Yb)^2}{n^2} \\
width &amp;= 12756 \frac{|L2'|}{K} \\
A &amp;= \frac{L1' - L2'}{L1' + L2'} \\
\end{align*}
$$</span></p>

CC BY-SA 4.0
2024-05-01T02:33:22.863Z
I've actually found your entry, and it is something I've been anxiously awaiting to help me understand it as well! My ultimate goal is to have this be something more accessible, not just to programmers, but even to historians and others who might be interested in such. I got as far as the first set of equations, but the rest will no doubt be quite helpful!
PearsonArtPhoto
CC BY-SA 4.0
2024-05-01T02:34:28.813Z
If you are interested, we could work together to put a library to do this together, perhaps in Python. Thanks again!
PearsonArtPhoto
CC BY-SA 4.0
2024-05-01T11:56:16.887Z
I have already written a library which implements most of the algorithms in the eclipse map linked above.  You can view the source from that page, or the Github repo: https://github.com/gmiller123456/solareclipseviewer .  The article would be explaining how to use it.
Greg Miller
CC BY-SA 4.0
2024-05-02T11:51:41.593Z
That library helped me find and clarify a few issues for your equations above, so hopefully they are right now!
PearsonArtPhoto
CC BY-SA 4.0
2024-05-02T12:36:34.597Z
@PearsonArtPhoto Thanks.  That's the good thing about code, you can actually test it to make sure it's right.  There's still a few issues like $ b1 $ should be $b_1$, but it's consistent, it just doesn't match the variables used in most texts.
Greg Miller
CC BY-SA 4.0
2024-04-30T21:38:15.547Z
2024-05-02T12:01:11.183Z
2024-05-02T12:01:11.183Z
Greg Miller
<p>NASA has published 5000 years of future solar eclipses at <a href="https://eclipse.gsfc.nasa.gov/SEpubs/5MCSE.html" rel="noreferrer">https://eclipse.gsfc.nasa.gov/SEpubs/5MCSE.html</a> . They have included in them a number of low fidelity maps, and the parameters to calculated them directly, specifically the Besselian elements. I have found from this site that there are a few key things, interpolating the values in time over a short period of time.</p>
<p>How can I use these elements to determine the path of an eclipse over the globe? Everything I have found that is online has some pieces of this, with reference to an out-of-print book that seems to have the other formulas, namely &quot;Elements of Solar Eclipses 1951-2200&quot;. Can anyone point me to such a resource, be it online or even via a recently in-print book? A software library, preferably in Python, would also be great!</p>

CC BY-SA 4.0
2024-04-29T02:49:25.200Z
This you might find interesting: https://www.mathworks.com/matlabcentral/fileexchange/71132-a-matlab-implementation-of-elements-of-solar-eclipses?s_tid=FX_rc2_behav
Pierre Paquette
CC BY-SA 4.0
2024-04-29T10:52:17.473Z
That might work. Annoying you have to have an account to download it, but...
PearsonArtPhoto
CC BY-SA 4.0
2024-04-29T13:38:48.473Z
Fred Espenak has retired from NASA, but he has his own site. Some of the data on Fred's site uses more precise calculations, and may cover a greater time span, eg https://eclipsewise.com/solar/SEcatalog/SEcatalog.html spans 6000 years.
PM 2Ring
CC BY-SA 4.0
2024-04-29T13:55:25.793Z
Fred gives 2 other references for Besselian elements, but they're even older. https://eclipsewise.com/solar/SEhelp/SEbeselm.html
PM 2Ring
CC BY-SA 4.0
2024-04-29T14:03:56.317Z
Will see what I can find. Good to see there is more information out there. I swear I'm going to make an open-source Python library that covers this when I am done...
PearsonArtPhoto
CC BY-SA 4.0
2024-04-30T00:47:16.050Z
The algorithms are published in a few places.  This is a public domain implementation of many of the relevant algorithms, sources are cited in the code.  Start with the "eclipse.js" file.  https://celestialprogramming.com/apps/SolarEclipseViewer/viewer.html
Greg Miller
CC BY-SA 4.0
2024-04-30T15:53:49.457Z
"Elements of Solar Eclipses" is not out of print.  The publisher closed a few years ago and Sky Publishing took over.  The book is availabe from them, though occasionally goes out of stock. https://shopatsky.com/products/elements-of-solar-eclipses-1951-2200
Greg Miller
CC BY-SA 4.0
2024-05-01T02:35:44.650Z
@uhoh I've looked, they are ~5 hours away at the nearest, but I did request it from interlibrary loan, so I'll eventually get it.
PearsonArtPhoto
<solar-eclipse><mathematics>
Plotting the path of an eclipse
PearsonArtPhoto
# Plotting the path of an eclipse

(Asked by: username_0 on 2024-04-29)

NASA has published 5000 years of future solar eclipses at <https://eclipse.gsfc.nasa.gov/SEpubs/5MCSE.html> . They have included in them a number of low fidelity maps, and the parameters to calculated them directly, specifically the Besselian elements. I have found from this site that there are a few key things, interpolating the values in time over a short period of time.

How can I use these elements to determine the path of an eclipse over the globe? Everything I have found that is online has some pieces of this, with reference to an out-of-print book that seems to have the other formulas, namely "Elements of Solar Eclipses 1951-2200". Can anyone point me to such a resource, be it online or even via a recently in-print book? A software library, preferably in Python, would also be great!

**Tags:** solar-eclipse, mathematics

> **Comments:**
> **username_1**: This you might find interesting: https://www.mathworks.com/matlabcentral/fileexchange/71132-a-matlab-implementation-of-elements-of-solar-eclipses?s_tid=FX_rc2_behav  
> **username_0**: That might work. Annoying you have to have an account to download it, but...  
> **username_2**: Fred Espenak has retired from NASA, but he has his own site. Some of the data on Fred's site uses more precise calculations, and may cover a greater time span, eg https://eclipsewise.com/solar/SEcatalog/SEcatalog.html spans 6000 years.  
> **username_2**: Fred gives 2 other references for Besselian elements, but they're even older. https://eclipsewise.com/solar/SEhelp/SEbeselm.html  
> **username_0**: Will see what I can find. Good to see there is more information out there. I swear I'm going to make an open-source Python library that covers this when I am done...  
> **username_3**: The algorithms are published in a few places.  This is a public domain implementation of many of the relevant algorithms, sources are cited in the code.  Start with the "eclipse.js" file.  https://celestialprogramming.com/apps/SolarEclipseViewer/viewer.html  
> **username_3**: "Elements of Solar Eclipses" is not out of print.  The publisher closed a few years ago and Sky Publishing took over.  The book is availabe from them, though occasionally goes out of stock. https://shopatsky.com/products/elements-of-solar-eclipses-1951-2200  
> **username_0**: @uhoh I've looked, they are ~5 hours away at the nearest, but I did request it from interlibrary loan, so I'll eventually get it.  

---

## 2 answers

### Answer by username_3 on 2024-04-30. Score: 4 | Accepted answer

As you've already found out there aren't many good resources for using Besselian Elements to compute solar eclipses. While there are some thorough resources, they are a bit terse, and careful study is required in order to figure out exactly which equations to use. Also Besselian Elements is a bit of a rabbit hole, as there are quite a few different algorithms for computing various phenomenon.

I have had it on my todo list to write a simplified article that can be followed by a typical programmer without much astronomy knowledge. So far, [I just have a stub](https://www.celestialprogramming.com/articles/eclipsePredictions/predictions.html) which has implementations of many examples from various sources. An article on [computing Besselian Elements](https://www.celestialprogramming.com/eclipsesGeneratingBesselianElements/generatingBesselianPolynomials.html). And an [example implementation of an eclipse map](https://www.celestialprogramming.com/apps/SolarEclipseViewer/viewer.html) which implements most of the algorithms.

Below I'll just answer the question you asked, which is computing the central line of an eclipse which is simple enough to reproduce here, and enough to get your feet wet implementing other algorithms.

In the interest of brevity, I'm not going to try to explain what each of the variables represent, as most of them are intermediate values without a simple explanation. For a more thorough coverage, consult one (or many) of the below resources which are what I used. The below algorithm is from Meeus.

* A Manual of Spherical and Practical Astronomy - Chauvenet 1863
* The Mathematical Theory of Eclipses According to Chauvenet's Transformation of Bessel's Method - Buchanan 1904
* Explanatory Supplement to the Astronomical Ephemeris 1961
* Explanatory Supplement to the Astronomical Almanac 1992
* Elements of Solar Eclipses - Meeus 1989
* Prediction and Analysis of Solar Eclipse Circumstances - Williams 1971

You will first need a set of Besselian Elements. These can be obtained from the [Five Millennium Catalog of Solar Eclipses](https://eclipse.gsfc.nasa.gov/SEcat5/SEcatalog.html). I have extracted [all of the elements from there](https://github.com/gmiller123456/FiveMillenniumCanonOfSolarEclipses-Besselian-Elements) and put them in a Github repository to make that process easier. The elements below are from the [2024 Apr 8 eclipse](https://eclipse.gsfc.nasa.gov/SEsearch/SEdata.php?Ecl=%2020240408).

| n | x | y | d | L1 | L2 | μ |
| --- | --- | --- | --- | --- | --- | --- |
| 0 | -0.3182440 | 0.2197640 | 7.5862002 | 0.5358140 | -0.0102720 | 89.591217 |
| 1 | 0.5117116 | 0.2709589 | 0.0148440 | 0.0000618 | 0.0000615 | 15.004080 |
| 2 | 0.0000326 | -0.0000595 | -0.0000020 | -0.0000128 | -0.0000127 | 0.000000 |
| 3 | -0.0000084 | -0.0000047 | 0.0000000 | 0.0000000 | 0.0000000 | 0.000000 |

```
tan f1 = 0.0046683        tan f2 = 0.0046450 
"Eclipse Predictions by Fred Espenak, NASA's GSFC"

```

In addition to the elements, some auxiliary information is also provided, namely $ t0 $, which is the basis time the elements were computed around, any time computed with the elements will be relative to $ t0 $, and is in hours. E.g. for 2024 Apr 8, $ t0 $ is 18.000, or 6pm Terrestrial Dynamic Time (TDT). There is also $ \Delta t $, used to convert TDT to UTC, and is in seconds.

The algorithm below computes the central line for a given time. So you specify a time $ t $, and it produces the latitude $ \phi $ and longitude $ \lambda $, and a few other values. Generally $t$ takes the values of -2 to +2, that's two hours before and two hours after $t0$.

The first step in most algorithms is to compute the elements for the specified time:
$$
\begin{align\*}
X &= X\_0 +X\_1t + X\_2t^2 + X\_3t^3 \\
Y &= Y\_0 + Y\_1t + Y\_2t^2 + Y\_3t^3 \\
d &= d\_0 + d\_1t + d\_2^2 \\
L1 &= L1\_0 + L1\_1t + L1\_2t^2 \\
L2 &= L2\_0 + L2\_1t + L2\_2t^2 \\
\mu &= \mu\_0 + \mu\_1t \\
\\
X' &= X\_1 + 2X\_2t + 3X\_3t^2 \\
Y' &= Y\_1 + 2Y\_2t + 3Y\_3t^2
\end{align\*}
$$

So, using the table above, $X\_0$ = -0.3182440, $X\_1$ = 0.5117116, $X\_2$ = 0.0000326, etc. Note that the trigonometric functions are expressed in degrees.

$$
\begin{align\*}
\omega &=\frac{1}{\sqrt{1 - 0.006694385 \cos^2 d}} \\
p &=\mu\_1/57.2957795 \\
b &=Y' - pX \sin d \\
c &=X' + pY \sin d \\
y1 &=\omega Y \\
b1 &=\omega \sin d \\
b2 &=0.99664719 \omega \cos d \\
B &= \sqrt{1 - X^2 - y1^2} \\
\end{align\*}
$$
If $B$ does not exist, there is no central eclipse at given time.
$$
\begin{align\*}
\sin \phi\_1 &=B\*b1 + y1\*b2 \\
\sin H &= \frac{X}{\cos \phi\_1} \\
\cos H &= \frac{B\*b2 - y1\*b1}{\cos \phi\_1} \\
\tan \phi &= 1.00336409 \tan \phi\_1 \\
\lambda &= \mu - H - 0.00417807 Δt \\
\end{align\*}
$$
$\lambda$ is the longitude (Positive = West, opposite normal convention), and $\phi$ is the latitude of the central point for the given time. H should be the number that is in the correct quadrant for the given $sin H$ and $cos H$. Additionally, the width, duration, and magnitude ($A$) can be computed without much more work. Note if duration is less than 0 then it is the duration of total eclipse, if greater than 0 it is the duration of annular eclipse.
$$
\begin{align\*}
L1'&=L1 - B \tan f1 \\
L2'&=L2 - B \tan f2 \\
a &= c - p\*B\* \cos d \\
n &= \sqrt{a^2 + b^2} \\
duration &= \frac{7200\*L2'}{n} \\
\sin h &=\sin d \sin \phi + \cos d \cos \phi \cos H \\
K^2 &= B^2 + \frac{(Xa + Yb)^2}{n^2} \\
width &= 12756 \frac{|L2'|}{K} \\
A &= \frac{L1' - L2'}{L1' + L2'} \\
\end{align\*}
$$

> **Comments:**
> **username_0**: I've actually found your entry, and it is something I've been anxiously awaiting to help me understand it as well! My ultimate goal is to have this be something more accessible, not just to programmers, but even to historians and others who might be interested in such. I got as far as the first set of equations, but the rest will no doubt be quite helpful!  
> **username_0**: If you are interested, we could work together to put a library to do this together, perhaps in Python. Thanks again!  
> **username_3**: I have already written a library which implements most of the algorithms in the eclipse map linked above.  You can view the source from that page, or the Github repo: https://github.com/gmiller123456/solareclipseviewer .  The article would be explaining how to use it.  
> **username_0**: That library helped me find and clarify a few issues for your equations above, so hopefully they are right now!  
> **username_3**: @username_0 Thanks.  That's the good thing about code, you can actually test it to make sure it's right.  There's still a few issues like $ b1 $ should be $b_1$, but it's consistent, it just doesn't match the variables used in most texts.  

---

### Answer by username_4 on 2024-05-02. Score: 1

Here's a solution that uses only Astropy and a guess for the time of the eclipse.
It doesn't use any Besselian elements, instead it finds the latitude and time that minimizes the separation between the Sun and Moon.

```
import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize
import astropy.units as u
import astropy.time
import astropy.coordinates
import astropy.constants
import cartopy

# Define the radius of the moon
radius_moon = 1737.4 * u.km

def get_sun_moon(
    location: astropy.coordinates.EarthLocation,
    time: astropy.time.Time,
) -> tuple[astropy.coordinates.SkyCoord, astropy.coordinates.SkyCoord]:
  """
  Given a location on the Earth and an observation time,
  get the coordinates of the Sun and Moon in the local
  coordinate system.

  Parameters
  ----------
  location
    A point on the surface of the earth
  time
    An observation time
  """
  frame_local = astropy.coordinates.AltAz(
      location=location,
      obstime=time,
  )
  sun = astropy.coordinates.get_body(
      body="sun",
      location=location,
      time=time,
      ephemeris="de440",
  ).transform_to(frame_local)
  moon = astropy.coordinates.get_body(
      body="moon",
      location=location,
      time=time,
      ephemeris="de440",
  ).transform_to(frame_local)
  return sun, moon

def sun_moon_separation(
    location: astropy.coordinates.EarthLocation,
    time: astropy.time.Time,
) -> u.Quantity:
  """
  Find the angular separation of the Sun and Moon
  at the given Earth location and time.

  Parameters
  ----------
  location
    A point on the surface of the Earth
  time
    An observation time
  """
  sun, moon = get_sun_moon(location, time)
  return sun.separation(moon)

def sun_moon_separation_argmin(
    location: astropy.coordinates.EarthLocation,
    time: astropy.time.Time,
) -> tuple[
    astropy.coordinates.EarthLocation,
    astropy.time.Time,
]:
  """
  Compute the time and latitude of the closest separation
  between the Sun and Moon using local minimization.

  Parameters
  ----------
  location
    A point on the Earth to use as an initial guess.
    The longitude is kept constant and the latitude is
    allowed to vary.
  time
    An initial guess for the time of closest approach.
  """
  def factory(
      lat: float, 
      t: float
  ) -> tuple[
      astropy.coordinates.EarthLocation,
      astropy.time.Time,
  ]:
    loc = astropy.coordinates.EarthLocation.from_geodetic(
        lat=lat * u.deg,
        lon=location.lon,
    )
    t = astropy.time.Time(t, format="jyear")
    return loc, t

  def fun(x: tuple[float, float]) -> float:
    lat, t = x
    loc, t = factory(lat, t)
    res = sun_moon_separation(loc, t)
    return res

  lat, t = scipy.optimize.minimize(
      fun=fun,
      x0=[location.lat.to_value(u.deg), time.jyear],
      bounds=[
          (-90, 90),
          (time.jyear - 0.01, time.jyear + 0.01),
      ],
      method="nelder-mead",
  ).x

  return factory(lat, t)

def calc_path_center(
    longitude: u.Quantity,
    time_guess: astropy.time.Time,
) -> tuple[
    astropy.coordinates.EarthLocation, 
    astropy.time.Time,
]:
  """
  Compute the center of the eclipse path across the Earth

  Parameters
  ----------
  longitude
    A grid of longitudes at which to compute the center of the eclipse
  time_guess
    An initial guess of the time at which the eclipse occurs.
  """
  time = astropy.time.Time(
      np.zeros(longitude.shape),
      format="jyear",
  )
  latitude = np.empty_like(longitude)

  lat_guess = 0 * u.deg

  for i in np.ndindex(longitude.shape):
    location_i, time_i = sun_moon_separation_argmin(
        location=astropy.coordinates.EarthLocation.from_geodetic(
            lat=lat_guess,
            lon=longitude[i],
        ),
        time=time_guess,
    )
    latitude[i] = lat_guess = location_i.lat
    time[i] = time_guess = time_i

  location = astropy.coordinates.EarthLocation.from_geodetic(
      lat=latitude,
      lon=longitude,
  )
  sun, moon = get_sun_moon(
      location=location,
      time=time,
  )
  separation = sun.separation(moon)
  angular_radius_sun = np.arctan2(
      astropy.constants.R_sun,
      sun.cartesian.norm(),
  )
  angular_radius_moon = np.arctan2(
      radius_moon,
      moon.cartesian.norm(),
  )

  where_eclipse = separation < (angular_radius_moon - angular_radius_sun)
  where_eclipse = where_eclipse & (sun.alt > 0)

  location = location[where_eclipse]
  time = time[where_eclipse]

  return location, time

# Define the start and stop times
# of the search
time_start = astropy.time.Time("2024-04-08T18:00")

# Define a grid of longitudes to search
longitude = np.linspace(-180, 180, num=101) * u.deg

# Compute the eclipse path
location_center, time_center = calc_path_center(
    longitude=longitude,
    time_guess=time_start,
)

# Plot the eclipse path
fig, ax = plt.subplots(
    subplot_kw=dict(projection=cartopy.crs.Robinson()),
)
ax.stock_img();
ax.coastlines();
ax.plot(
    location_center.lon,
    location_center.lat,
    transform=cartopy.crs.PlateCarree(),
);

```

[![April 2024 Eclipse](https://i.sstatic.net/XIJnUAzc.png)](https://i.sstatic.net/XIJnUAzc.png)

> **Comments:**
> **username_0**: I started with a solution like this, but the accuracy using Astropy was off by about a degree or so, and I ended up giving up on it.  
> **username_4**: Astropy should really be just as close as anything else, if it's not that would be a big problem. Were you using the default ephemeris? I would expect the default ephemeris (as opposed to the JPL ephemeris) to be off by a little bit, but not a whole degree.  
> **username_0**: I tried both JPL and the default, they were both fairly similar, but off. Maybe my code was off, however...  
> **username_0**: Not sure what I am doing, but I get `25.613N, 103.832 W time 2024-04-08T18:18:29` as the peak with my effectively identical method, while the peak location at that time according to the catalog was 25.3 N, 104.1 W. This with the de440.  
> **username_4**: It seems like the numbers `25.613N, 103.832 W time 2024-04-08T18:18:29` match the values on @username_3's site, can you show me what tabulated values you're using?  
> **username_4**: Also, confusingly, my code right now gives `25.802 N, 103.832 W time 2024-04-08T18:20:02.846` for that longitude, so I'm not sure why our codes are different.  
> **username_0**: I don't have exactly the same code, just functionally the same (Same key functions, different working goods). I'm comparing it to the values reported in the NASA catalog for the event.  

---

