<p>Your methodology seems to be based on an algorithm that uses successive approximations to obtain the desired values. I was not able to arrive quite at the same values as you mentioned, but I noticed an error in the posted script, <code>zeta = sqrt(ro^2 + x^2 + y^2)</code> clearly should be <code>zeta = sqrt(ro^2 - x^2 - y^2)</code>, so perhaps slightly different values or operations were used somewhere else too, but I couldn't find out where. I also noticed that you did not take into account the difference between Ephemeris Longitude and True Longitude in your calculation, True Longitude shifts towards the East by <code>1.002738 * At (Delta t) / 240</code> degrees. Moreover, executing more iterations could possibly converge to a better result.</p>

<p>However, to compute the intersection of the Moon's shadow axis from Besselian Elements, I suggest using the following closed-form solution, which is derived from solving the intersection of a parametric line (Moon's shadow axis) and the Earth ellipsoid.</p>

<pre><code>f = 1 / 298.257223563
b = 1 - f
sin_d = Sin(d * PI / 180)
cos_d = Cos(d * PI / 180)
y0_r = y * cos_d
z0_r = -y * sin_d
y1_r = y0_r + sin_d
z1_r = z0_r + cos_d
v = y1_r - y0_r
w = z1_r - z0_r
a_quadratic = v ^ 2 + b ^ 2 * w ^ 2
b_quadratic = 2 * z0_r * w * b ^ 2 + 2 * y0_r * v
c_quadratic = z0_r ^ 2 * b ^ 2 + y0_r ^ 2 + b ^ 2 * (x ^ 2 - 1)
p = (-b_quadratic + Sqrt(b_quadratic ^ 2 - 4 * a_quadratic * c_quadratic)) / (2 * a_quadratic)
y_i = y0_r + v * p
z_i = z0_r + w * p
dist_axis = Sqrt(x ^ 2 + z_i ^ 2)
lat = Atan(y_i / (dist_axis * (1 - f) ^ 2)) * 180 / PI
lon = Atan2(x, z_i) * 180 / PI - u + 1.002738 * At / 240
</code></pre>

<p>Here is the explanation of the formulas:</p>

<p>We use the flattening value to define the WGS84 ellipsoid, and calculate the semiminor-axis:</p>

<pre><code>f = 1 / 298.257223563
b = 1 - f
</code></pre>

<p>Convert declination angle d to radians and execute sine and cosine functions:</p>

<pre><code>sin_d = Sin(d * PI / 180)
cos_d = Cos(d * PI / 180)
</code></pre>

<p>Select 2 points (points 0 and 1 in the code) on the shadow axis, at z = 0 and z = 1, and apply the declination rotation (_r) on them:</p>

<pre><code>y0_r = y * cos_d
z0_r = -y * sin_d
y1_r = y0_r + sin_d
z1_r = z0_r + cos_d
</code></pre>

<p>The x value always remains the same, so we don't need to recalculate it. 
Next, evaluate the parametric line's vectors v,w (delta y, z) from point 0 as a base:</p>

<pre><code>v = y1_r - y0_r
w = z1_r - z0_r
</code></pre>

<p>By substituting the parametric line's equations in the ellipsoid equation, and rearranging to quadratic form and solve:</p>

<pre><code>a_quadratic = v ^ 2 + b ^ 2 * w ^ 2
b_quadratic = 2 * z0_r * w * b ^ 2 + 2 * y0_r * v
c_quadratic = z0_r ^ 2 * b ^ 2 + y0_r ^ 2 + b ^ 2 * (x ^ 2 - 1)
p = (-b_quadratic + Sqrt(b_quadratic ^ 2 - 4 * a_quadratic * c_quadratic)) / (2 * a_quadratic)
</code></pre>

<p>'p' here represents the distance from point 0 along the parametric line, expressed in Earth's equatorial radius.
The contents of the square root must ideally be tested beforehand for negative values to avoid exceptions if the line does not intersect the ellipsoid.
With this, we find the coordinates of our points of intersection:</p>

<pre><code>y_i = y0_r + v * p
z_i = z0_r + w * p
</code></pre>

<p>Then we convert this to geographic latitude and longitude:</p>

<pre><code>dist_axis = Sqrt(x ^ 2 + z_i ^ 2)
lat = Atan(y_i / (dist_axis * (1 - f) ^ 2)) * 180 / PI
lon = Atan2(x, z_i) * 180 / PI - u + 1.002738 * At / 240
</code></pre>

<p>Direct conversion for the Latitude is possible here because our point is on the surface of the ellipsoid.</p>

<p>At the end, you can use the following logic to capture out-of-bounds longitude values:</p>

<pre><code>If lon &lt;= -180 Then
    lon += 360
ElseIf lon &gt; 180 Then
    lon -= 360
End If
</code></pre>

<p>In terms of accuracy, you can expect pretty good results with this method, but it does not take into account some factors, like local topography, geoid height, refraction, etc. Topography can lead to differences of a few kilometers, the geoid, +-100 meters, and refraction can bend the shadow by a few meters to a few hundred meters depending on the incident angle. Also, there is no correction for the moon's center of figure, that can shift the umbral center by several hundred meters as well. </p>

<p>With the above formulas for UT1 = 18h and delta t = 72.3s, the computed intersection point is at 58.26828°N, 21.57331°W, which matches F. Espenak's values. </p>

CC BY-SA 4.0
2020-12-28T15:54:21.727Z
This is such an excellent piece of coding. Sorry but I don't fully understand all of the terms used, I think that d represents the declination of the moon but what does y refer to in the equations `y0_r = y * cos_d` and `z0_r = -y * sin_d` please. Thanks for your expertise and help Philip
Philip Sterling
CC BY-SA 4.0
2021-01-03T23:25:32.377Z
@PhilipSterling Thanks for your kind words. In Besselian elements and in my code,`d` represents the declination of the _direction of the shadow axis_ so it is not the same as the traditional geocentric declination of the moon. The Y in my code represents the y-coordinate of the moon's shadow axis on the fundamental plane. See [this image](https://en.wikipedia.org/wiki/Besselian_elements#/media/File:BesselianElementsForTotalEclipse3D-en.svg), it might help you visualize what the x,y variables represent. Feel free to ask if you have any other questions.
FSimardGIS
CC BY-SA 4.0
2021-01-03T23:43:25.483Z
To obtain the value of `y`, you can refer to the table of Besselian elements for your specific eclipse, and compute it just like in the original question: `y = y0 + y1 * t + y2 * t^2 + y3 * t^3`. Same thing for `x` which you'll need later on.
FSimardGIS
CC BY-SA 4.0
2018-10-03T03:53:45.057Z
2018-10-03T19:06:06.707Z
2018-10-03T19:06:06.707Z
FSimardGIS
<p>In your formula for lat:</p>
<p><span class="math-container">$$
lat = Atan\bigg(\frac{y_i}{(dist_{axis} \times (1 - f) ^ 2)} \bigg) \times \frac{180 }{ \pi}
$$</span></p>
<p>I think <span class="math-container">$(1-f)^2$</span> should be just <span class="math-container">$(1-f)$</span>, because <span class="math-container">$\sqrt{1-\epsilon^2}=1-f$</span>,
since <span class="math-container">$1-\epsilon^2=\frac{b^2}{a^2}=(1-f)^2$</span>.</p>
<p>See for instance Chauvenet, A Manual of Spherical and Practical Astronomy, formula 539, p. 491.</p>

CC BY-SA 4.0
2021-01-04T16:57:22.240Z
That is because Chauvenet uses the tangent of a reduced latitude `phi_1` in formula 539. In my formula, I use `y_i / dist_axis` which is equivalent to the tangent of the geocentric latitude, not reduced. So that's why I had to use `(1-f)^2`
FSimardGIS
CC BY-SA 4.0
2021-01-04T10:12:21.857Z
2021-01-04T12:31:16.390Z
2021-01-04T12:31:16.390Z
priis
<p>I'm trying to calculate <strong>the closest path</strong> of the center line of a given eclipse at a given time by using its besselian elements; However i'm getting an error range up to 30 km (18.6 mile) in the worst case and 5 km (3.1 mile) in the best case (both cases was compared with F. Espenak prediction)</p>
<p>Lets use eclipse of August 12, 2026 at 18:00 UTC.</p>
<p>F. Espenak <a href="http://www.eclipsewise.com/solar/SEpath/2001-2100/SE2026Aug12Tpath.html" rel="nofollow noreferrer">prediction</a>: 58°16.1'N   021°34.4'W (58.2683, 21.5733)</p>
<p>My prediction: 58°16.8'N 21°47.4'W (58.2795, 21.7901)</p>
<p>The difference is close to 13 km (8 mile) which is too big for my application so I need to refine/review the procedure/algorithm I'm using.</p>
<h2>Calculus</h2>
<p>Given the following besselian elements:</p>
<pre><code>e20260812 = {
    t0: 18.00, At: 72.3,
    x0: 0.475510, x1: 0.518923, x2: -7.727e-05, x3: -8.039e-06,
    y0: 0.771185, y1: -0.230167, y2: -1.246e-04, y3: 3.767e-06,
    d0: 14.796683, d1: -0.012065, d2: -3.098e-06, d3: 0.0,
    l10: 0.537974, l11: 0.000094, l12: -1.212e-05, l13: 0.0,
    l20: -0.008142, l21: 0.000093, l22: -1.206e-05, l23: 0.0,
    u0: 88.747810, u1: 15.003090, u2: 1.764e-06, u3: 0.0,
    tanf1: 0.004614,  tanf2: 0.004591
}
</code></pre>
<p>(<code>At</code> represents <span class="math-container">$\Delta T$</span>)</p>
<h3>0 - Reference timestamp.</h3>
<p>Where <code>t1 = 18.00</code></p>
<pre><code>t = t1 - t0 + At/3600.0
</code></pre>
<h3>1 - Time dependent calculus</h3>
<pre><code>x = x0 + x1 * t + x2 * t^2 + x3 * t^3
y = y0 + y1 * t + y2 * t^2 + y3 * t^3
d = d0 + d1 * t + d2 * t^2 + d3 * t^3
u = u0 + u1 * t + u2 * t^2 + u3 * t^3 
</code></pre>
<h3>2 - Location dependent calculus</h3>
<p>Initially I'm using <code>ro=1.0</code> and following the below sequence to obtain an approximation of <code>phi</code>, which is the Latitude.</p>
<pre><code>ro = 1.0
e = 0.0818192 // earth eccentricity 
sin_d = sin(d * pi / 180.0)
cos_d = cos(d * pi / 180.0)
zeta = sqrt(ro^2 + x^2 + y^2)
cos_theta = (-y * sin_d + zeta * cos_d) / ro
sin_phi_1 = (y * cos_d + zeta * sen_d) / ro
cos_phi_1 = sqrt( (x/ro)^2 + cos_theta^2 )
phi_1 = asin(sin_phi_1) / pi * 180.0
phi = atan(sin_phi_1 / cos_phi_1 / sqrt(1 - e^2)) / pi * 180.0
</code></pre>
<p>At this point <code>phi = 58.1224</code> for 18:00 hs UTC.</p>
<p>Once <code>phi</code> was calculated, we recalculate <code>ro</code> as follow to get a better approximation of <code>phi</code> and <code>varpi</code> (not yet calculated)</p>
<pre><code>sin_phi = sin(phi * pi / 180.0)
ro1 = sqrt( (1 - 2e^2 * sin_phi^2 + e^4 * sin_phi^2) / (1 - e^2 * sin_phi^2) )
</code></pre>
<p>Applying the updated <code>ro</code> (now <code>ro1</code>), we update <code>sin_d</code> and <code>cos_d</code> as follow:</p>
<pre><code>sin_d1 = sin_d / ro1;
cos_d1 = sqrt(1-e^2) * (cos_d / ro1)
</code></pre>
<p>And update <code>zeta</code> and <code>cos_theta</code></p>
<pre><code>zeta1 = sqrt(ro1^2 + x^2 + y^2)
cos_theta1 = (-y * sin_d1 + zeta1 * cos_d1) / ro1
</code></pre>
<p>Finally we calculate <code>vapri</code> (the Longitud) as follows:</p>
<pre><code>tan_theta = x / cos_theta1
theta = atan( tan_theta ) / Math.PI * 180.0
vapri = u - theta * ro1
</code></pre>
<p>and recalculate <code>phi</code> by doing:</p>
<pre><code>cos_phi2 =  cos_theta1 / cos(theta / 180.0 * pi)
phi1 = acos(cos_phi2) / pi * 180.0 
phi = atan(tan(phi1 / 180.0 * pi) / sqrt(1 - e^2)) / pi * 180.0
</code></pre>
<p>At this point <code>phi = 58.2795 , vapri = 21.7901</code> for 18:00 hs UTC (both values MUST be properly signed but I omitted that trivial part) a significant difference compared with the NASA prediction: 58.2683, 21.5733 (13 km or 8 mile)</p>
<p>My feeling is that <code>ro1</code> in not good enough precise, so the final results get shifted, but I have no idea at all about how to getting unblocked.</p>
<p>I'm using 64 bit floating point calculator; It is ok for this kind of calculus?</p>
<p>Is this method good enough to get an error range near to 500m? if not, could you provide a method/fix to achieve that aim.</p>
<p>Any help would be appreciated.</p>

CC BY-SA 3.0
2016-09-28T17:27:18.720Z
Hi @barrycarter, thanks for commenting. As I mentioned, I'm open to other methods.
rnrneverdies
<solar-eclipse>
Closest path of solar eclipse central line
rnrneverdies
# Closest path of solar eclipse central line

(Asked by: username_0 on 2016-09-24)

I'm trying to calculate the closest path of the center line of a given eclipse at a given time by using its besselian elements; However i'm getting an error range up to 30 km (18.6 mile) in the worst case and 5 km (3.1 mile) in the best case (both cases was compared with F. Espenak prediction)
Lets use eclipse of August 12, 2026 at 18:00 UTC.
F. Espenak prediction: 58°16.1'N   021°34.4'W (58.2683, 21.5733)
My prediction: 58°16.8'N 21°47.4'W (58.2795, 21.7901)
The difference is close to 13 km (8 mile) which is too big for my application so I need to refine/review the procedure/algorithm I'm using.
Calculus
Given the following besselian elements:
e20260812 = {
    t0: 18.00, At: 72.3,
    x0: 0.475510, x1: 0.518923, x2: -7.727e-05, x3: -8.039e-06,
    y0: 0.771185, y1: -0.230167, y2: -1.246e-04, y3: 3.767e-06,
    d0: 14.796683, d1: -0.012065, d2: -3.098e-06, d3: 0.0,
    l10: 0.537974, l11: 0.000094, l12: -1.212e-05, l13: 0.0,
    l20: -0.008142, l21: 0.000093, l22: -1.206e-05, l23: 0.0,
    u0: 88.747810, u1: 15.003090, u2: 1.764e-06, u3: 0.0,
    tanf1: 0.004614,  tanf2: 0.004591
}

(At represents $\Delta T$)
0 - Reference timestamp.
Where t1 = 18.00
t = t1 - t0 + At/3600.0

1 - Time dependent calculus
x = x0 + x1 * t + x2 * t^2 + x3 * t^3
y = y0 + y1 * t + y2 * t^2 + y3 * t^3
d = d0 + d1 * t + d2 * t^2 + d3 * t^3
u = u0 + u1 * t + u2 * t^2 + u3 * t^3 

2 - Location dependent calculus
Initially I'm using ro=1.0 and following the below sequence to obtain an approximation of phi, which is the Latitude.
ro = 1.0
e = 0.0818192 // earth eccentricity 
sin_d = sin(d * pi / 180.0)
cos_d = cos(d * pi / 180.0)
zeta = sqrt(ro^2 + x^2 + y^2)
cos_theta = (-y * sin_d + zeta * cos_d) / ro
sin_phi_1 = (y * cos_d + zeta * sen_d) / ro
cos_phi_1 = sqrt( (x/ro)^2 + cos_theta^2 )
phi_1 = asin(sin_phi_1) / pi * 180.0
phi = atan(sin_phi_1 / cos_phi_1 / sqrt(1 - e^2)) / pi * 180.0

At this point phi = 58.1224 for 18:00 hs UTC.
Once phi was calculated, we recalculate ro as follow to get a better approximation of phi and varpi (not yet calculated)
sin_phi = sin(phi * pi / 180.0)
ro1 = sqrt( (1 - 2e^2 * sin_phi^2 + e^4 * sin_phi^2) / (1 - e^2 * sin_phi^2) )

Applying the updated ro (now ro1), we update sin_d and cos_d as follow:
sin_d1 = sin_d / ro1;
cos_d1 = sqrt(1-e^2) * (cos_d / ro1)

And update zeta and cos_theta
zeta1 = sqrt(ro1^2 + x^2 + y^2)
cos_theta1 = (-y * sin_d1 + zeta1 * cos_d1) / ro1

Finally we calculate vapri (the Longitud) as follows:
tan_theta = x / cos_theta1
theta = atan( tan_theta ) / Math.PI * 180.0
vapri = u - theta * ro1

and recalculate phi by doing:
cos_phi2 =  cos_theta1 / cos(theta / 180.0 * pi)
phi1 = acos(cos_phi2) / pi * 180.0 
phi = atan(tan(phi1 / 180.0 * pi) / sqrt(1 - e^2)) / pi * 180.0

At this point phi = 58.2795 , vapri = 21.7901 for 18:00 hs UTC (both values MUST be properly signed but I omitted that trivial part) a significant difference compared with the NASA prediction: 58.2683, 21.5733 (13 km or 8 mile)
My feeling is that ro1 in not good enough precise, so the final results get shifted, but I have no idea at all about how to getting unblocked.
I'm using 64 bit floating point calculator; It is ok for this kind of calculus?
Is this method good enough to get an error range near to 500m? if not, could you provide a method/fix to achieve that aim.
Any help would be appreciated.

**Tags:** solar-eclipse
> **Comments:**
> **username_0**: Hi @barrycarter, thanks for commenting. As I mentioned, I'm open to other methods.  

---

## 2 answers

### Answer by username_1 on 2018-10-03. Score: 2

Your methodology seems to be based on an algorithm that uses successive approximations to obtain the desired values. I was not able to arrive quite at the same values as you mentioned, but I noticed an error in the posted script, zeta = sqrt(ro^2 + x^2 + y^2) clearly should be zeta = sqrt(ro^2 - x^2 - y^2), so perhaps slightly different values or operations were used somewhere else too, but I couldn't find out where. I also noticed that you did not take into account the difference between Ephemeris Longitude and True Longitude in your calculation, True Longitude shifts towards the East by 1.002738 * At (Delta t) / 240 degrees. Moreover, executing more iterations could possibly converge to a better result.
However, to compute the intersection of the Moon's shadow axis from Besselian Elements, I suggest using the following closed-form solution, which is derived from solving the intersection of a parametric line (Moon's shadow axis) and the Earth ellipsoid.
f = 1 / 298.257223563
b = 1 - f
sin_d = Sin(d * PI / 180)
cos_d = Cos(d * PI / 180)
y0_r = y * cos_d
z0_r = -y * sin_d
y1_r = y0_r + sin_d
z1_r = z0_r + cos_d
v = y1_r - y0_r
w = z1_r - z0_r
a_quadratic = v ^ 2 + b ^ 2 * w ^ 2
b_quadratic = 2 * z0_r * w * b ^ 2 + 2 * y0_r * v
c_quadratic = z0_r ^ 2 * b ^ 2 + y0_r ^ 2 + b ^ 2 * (x ^ 2 - 1)
p = (-b_quadratic + Sqrt(b_quadratic ^ 2 - 4 * a_quadratic * c_quadratic)) / (2 * a_quadratic)
y_i = y0_r + v * p
z_i = z0_r + w * p
dist_axis = Sqrt(x ^ 2 + z_i ^ 2)
lat = Atan(y_i / (dist_axis * (1 - f) ^ 2)) * 180 / PI
lon = Atan2(x, z_i) * 180 / PI - u + 1.002738 * At / 240

Here is the explanation of the formulas:
We use the flattening value to define the WGS84 ellipsoid, and calculate the semiminor-axis:
f = 1 / 298.257223563
b = 1 - f

Convert declination angle d to radians and execute sine and cosine functions:
sin_d = Sin(d * PI / 180)
cos_d = Cos(d * PI / 180)

Select 2 points (points 0 and 1 in the code) on the shadow axis, at z = 0 and z = 1, and apply the declination rotation (_r) on them:
y0_r = y * cos_d
z0_r = -y * sin_d
y1_r = y0_r + sin_d
z1_r = z0_r + cos_d

The x value always remains the same, so we don't need to recalculate it. 
Next, evaluate the parametric line's vectors v,w (delta y, z) from point 0 as a base:
v = y1_r - y0_r
w = z1_r - z0_r

By substituting the parametric line's equations in the ellipsoid equation, and rearranging to quadratic form and solve:
a_quadratic = v ^ 2 + b ^ 2 * w ^ 2
b_quadratic = 2 * z0_r * w * b ^ 2 + 2 * y0_r * v
c_quadratic = z0_r ^ 2 * b ^ 2 + y0_r ^ 2 + b ^ 2 * (x ^ 2 - 1)
p = (-b_quadratic + Sqrt(b_quadratic ^ 2 - 4 * a_quadratic * c_quadratic)) / (2 * a_quadratic)

'p' here represents the distance from point 0 along the parametric line, expressed in Earth's equatorial radius.
The contents of the square root must ideally be tested beforehand for negative values to avoid exceptions if the line does not intersect the ellipsoid.
With this, we find the coordinates of our points of intersection:
y_i = y0_r + v * p
z_i = z0_r + w * p

Then we convert this to geographic latitude and longitude:
dist_axis = Sqrt(x ^ 2 + z_i ^ 2)
lat = Atan(y_i / (dist_axis * (1 - f) ^ 2)) * 180 / PI
lon = Atan2(x, z_i) * 180 / PI - u + 1.002738 * At / 240

Direct conversion for the Latitude is possible here because our point is on the surface of the ellipsoid.
At the end, you can use the following logic to capture out-of-bounds longitude values:
If lon <= -180 Then
    lon += 360
ElseIf lon > 180 Then
    lon -= 360
End If

In terms of accuracy, you can expect pretty good results with this method, but it does not take into account some factors, like local topography, geoid height, refraction, etc. Topography can lead to differences of a few kilometers, the geoid, +-100 meters, and refraction can bend the shadow by a few meters to a few hundred meters depending on the incident angle. Also, there is no correction for the moon's center of figure, that can shift the umbral center by several hundred meters as well. 
With the above formulas for UT1 = 18h and delta t = 72.3s, the computed intersection point is at 58.26828°N, 21.57331°W, which matches F. Espenak's values.

> **Comments:**
> **username_2**: This is such an excellent piece of coding. Sorry but I don't fully understand all of the terms used, I think that d represents the declination of the moon but what does y refer to in the equations `y0_r = y * cos_d` and `z0_r = -y * sin_d` please. Thanks for your expertise and help Philip  
> **username_1**: @PhilipSterling Thanks for your kind words. In Besselian elements and in my code,`d` represents the declination of the _direction of the shadow axis_ so it is not the same as the traditional geocentric declination of the moon. The Y in my code represents the y-coordinate of the moon's shadow axis on the fundamental plane. See [this image](https://en.wikipedia.org/wiki/Besselian_elements#/media/File:BesselianElementsForTotalEclipse3D-en.svg), it might help you visualize what the x,y variables represent. Feel free to ask if you have any other questions.  
> **username_1**: To obtain the value of `y`, you can refer to the table of Besselian elements for your specific eclipse, and compute it just like in the original question: `y = y0 + y1 * t + y2 * t^2 + y3 * t^3`. Same thing for `x` which you'll need later on.  

---

### Answer by username_3 on 2021-01-04. Score: 2

In your formula for lat:

$$
lat = Atan\bigg(\frac{y\_i}{(dist\_{axis} \times (1 - f) ^ 2)} \bigg) \times \frac{180 }{ \pi}
$$

I think $(1-f)^2$ should be just $(1-f)$, because $\sqrt{1-\epsilon^2}=1-f$,
since $1-\epsilon^2=\frac{b^2}{a^2}=(1-f)^2$.

See for instance Chauvenet, A Manual of Spherical and Practical Astronomy, formula 539, p. 491.

> **Comments:**
> **username_1**: That is because Chauvenet uses the tangent of a reduced latitude `phi_1` in formula 539. In my formula, I use `y_i / dist_axis` which is equivalent to the tangent of the geocentric latitude, not reduced. So that's why I had to use `(1-f)^2`  

---

