Closest path of solar eclipse central line

Question

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.

Answer 1

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.

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