This got me interested to know if there is a gold standard for such calculations
For the orbital mechanics and geometrical aspects, the NASA JPL Development Ephemerides are commonly used. These have errors for the Sun and inner planets of tens of kilometers or less.
...and what is the range of error for it. For forecasts n years out is the error range n times more?
That depends on the specific ephemeris you use and a lot of other details. The ephemerides are build from a huge body of observations and calculations and it's not a simple thing to say. I think the errors will be dominated by the definition of the edge of the Sun for several decades or even centuries. Farther out in time, the ephemeris itself will have some error that may increase with time in either direction.
Below I calculate that the apparent radial speed at which mercury crosses the "edge" of the Sun's disk is about -70.6 and +70.6 km/sec, assuming a "radius" of 695,700 km, using the Python package Skyfield which in turn uses the NASA JPL Development Ephemerides. For this calculation I used DE421.
I used Caracas as the location since the whole event was visible from this location. The main calculation happens in these three lines. The .observe() method automatically does the light-time correction for the positions.
sun_astrometric = caracas.at(times).observe(sun)
mercury_astrometric = caracas.at(times).observe(mercury)
angular_separation = sun_astrometric.separation_from(mercury_astrometric).radians
I put the words "edge" and "radius" in quotes earlier because these are not well defined for things like stars and gas giant planets. There are definitions that can be used to define nominals for these, in the case of the Sun it's when the opacity reaches 2/3 or something like that.
70 km is a pretty small distance at the edge of the Sun. From this answer you can see that there is a region of a few hundred kilometers where things are changing. It's hard to say exactly what defines the "edge" of the Sun in order to say when the hard edge of rocky Mercury will be seen to touch it. I can't be sure, but I have a hunch that you can get several seconds of ambiguity based on how the edge of the Sun is defined.
Wikipedia even lists **two different values for the radius that differ by 642 km*!* At 70 km/sec that's almost ten seconds right there!
Equatorial radius
695,700 km [7]
696,342 km [8]
From this answer:
From Chapter 2: The Photosphere of Timo Nieminen's thesis Solar Line Asymmetries: Modelling the Effect of Granulation on the Solar Spectrum
Figure 2-3: The Holweger-Müller Model Atmosphere
7 Holweger, H. and Müller, E. A. “The Photospheric Barium Spectrum: Solar Abundance and Collision Broadening of Ba II Lines by Hydrogen”, Solar Physics 39, pg 19-30 (1974). Extra points have been cubic spline interpolated by J. E. Ross. The optical properties (such as the optical depth and the opacity) of a model atmosphere are, obviously, very important, and will be considered later. See table C-4 for complete details of the Holweger-Müller model atmosphere including all depth points used.
8The height scale is not arbitrary. The base of the photosphere (height = 0 km) is chosen to be at standard optical depth of one (i.e. 𝜏 5000Å = 1 ).
Python script using the Skyfield package to calculate some things about the 11-Nov-2019 transit of Mercury across the Sun as seen from the city of Caracas, Venezuela.
import numpy as np
import matplotlib.pyplot as plt
from skyfield.api import Topos, Loader
# https://astronomy.stackexchange.com/questions/33918/accuracy-of-mercury-transit-calculations
load = Loader('~/Documents/fishing/SkyData') # avoids multiple copies of large files
data = load('de421.bsp')
ts = load.timescale()
bodies = ('sun', 'mercury', 'earth')
seconds = 12*3600. + np.arange(0, 390*60+1, 10)
hours = seconds/3600.
times = ts.utc(2019, 11, 11, hours)
sun, mercury, earth = [data[x] for x in bodies]
caracas = earth + Topos(latitude_degrees = 10.48,
longitude_degrees = -66.90,
elevation_m = 900.)
sun_astrometric = caracas.at(times).observe(sun)
mercury_astrometric = caracas.at(times).observe(mercury)
angular_separation = sun_astrometric.separation_from(mercury_astrometric).radians
distance_to_sun = sun_astrometric.distance().km
distance = distance_to_sun * np.tan(angular_separation)
r_sun = 695700*np.ones_like(hours)
thing = distance < r_sun
i_in = np.argmax(thing)
thing[:i_in+1] = True
i_out = np.argmax(~thing)
hours_in, hours_out = hours[[i_in, i_out]]
dt = seconds[1] - seconds[0]
ddistance_dt = (distance[1:] - distance[:-1])/dt
if True:
plt.figure()
fs = 12
plt.subplot(2, 1, 1)
plt.plot(hours, distance)
plt.plot(hours[[i_in, i_out]], distance[[i_in, i_out]], 'ok')
plt.plot(hours, r_sun, '-k')
plt.ylim(0, None)
plt.title('projected distance from Sun center (km,)', fontsize=fs)
plt.subplot(2, 1, 2)
plt.plot(hours[:-1], ddistance_dt)
plt.plot(hours[[i_in, i_out]], ddistance_dt[[i_in, i_out]], 'ok')
plt.xlabel('hours UTC 2019-11-11', fontsize=fs)
plt.title('d projected distance/dt (km/s)', fontsize=fs)
print(ddistance_dt[[i_in, i_out]])
plt.show()
