# Closest path of solar eclipse central line

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.

• I could show you how to calculate this, but it would be an entirely different non-Besselian method so wouldn't really answer your question. – user21 Sep 28 '16 at 16:32
• Hi @barrycarter, thanks for commenting. As I mentioned, I'm open to other methods. – rnrneverdies Sep 28 '16 at 17:27
• Are you familiar with pyephem, skyfield, SPICE libraries and other tools used for precision astronomy? I was going to suggest one of those. – user21 Sep 28 '16 at 20:54

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)
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' 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.

• 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 Dec 28 '20 at 15:54
• @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, it might help you visualize what the x,y variables represent. Feel free to ask if you have any other questions. – FSimardGIS Jan 3 at 23:25
• 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 Jan 3 at 23:43

$$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$$.
• 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 Jan 4 at 16:57