<blockquote>
  <p>This got me interested to know if there is a gold standard for such calculations</p>
</blockquote>

<p>For the orbital mechanics and geometrical aspects, the <a href="https://en.wikipedia.org/wiki/Jet_Propulsion_Laboratory_Development_Ephemeris" rel="nofollow noreferrer">NASA JPL Development Ephemerides</a> are commonly used. These have errors for the Sun and inner planets of tens of kilometers or less. </p>

<ul>
<li><a href="https://english.stackexchange.com/q/370655/217285">How to pronounce “Ephemerides”?</a></li>
</ul>

<blockquote>
  <p>...and what is the range of error for it. For forecasts n years out is the error range n times more?</p>
</blockquote>

<p>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.</p>

<p>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 <a href="https://rhodesmill.org/skyfield/" rel="nofollow noreferrer">Skyfield</a> which in turn uses the NASA JPL Development Ephemerides. For this calculation I used DE421.</p>

<p>I used Caracas as the location since the whole event was visible from this location. The main calculation happens in these three lines. The <code>.observe()</code> method automatically does the light-time correction for the positions.</p>

<pre><code>sun_astrometric     = caracas.at(times).observe(sun) 
mercury_astrometric = caracas.at(times).observe(mercury)
angular_separation  = sun_astrometric.separation_from(mercury_astrometric).radians
</code></pre>

<hr>

<p>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.</p>

<p>70 km is a pretty small distance at the edge of the Sun. From <a href="https://astronomy.stackexchange.com/a/32737/7982">this answer</a> 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.</p>

<p>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!</p>

<pre><code>Equatorial radius   
  695,700 km  [7]
  696,342 km  [8]
</code></pre>

<p>From <a href="https://astronomy.stackexchange.com/a/32737/7982">this answer</a>:</p>

<blockquote>
  <p>From <a href="https://people.smp.uq.edu.au/TimoNieminen/papers/thesis/chapter2_photosphere.pdf" rel="nofollow noreferrer">Chapter 2: The Photosphere</a> of Timo Nieminen's thesis <a href="https://people.smp.uq.edu.au/TimoNieminen/thesis.html" rel="nofollow noreferrer">Solar Line Asymmetries: Modelling the Effect of Granulation on the Solar Spectrum</a></p>
  
  <p><a href="https://i.sstatic.net/JS8JR.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/JS8JR.png" alt="density versus hight above the Sun&#39;s surface"></a></p>
  
  <p><a href="https://i.sstatic.net/s6zLZ.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/s6zLZ.png" alt="density and other parameters versus hight above the Sun&#39;s surface"></a></p>
  
  <blockquote>
    <p>Figure 2-3: The Holweger-Müller Model Atmosphere</p>
    
    <p><sup>7</sup> 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.</p>
    
    <p><sup>8</sup><strong>The height scale is not arbitrary.</strong> The base of the photosphere (height = 0 km) is chosen to be at standard optical depth of one (i.e.  <sub>5000Å</sub> = 1 ).</p>
  </blockquote>
</blockquote>

<p><a href="https://i.sstatic.net/QAoRh.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/QAoRh.png" alt="approximate transit of Venus 11-11-2019 from Caracas"></a></p>

<p>Python script using the <a href="https://rhodesmill.org/skyfield/" rel="nofollow noreferrer">Skyfield</a> package to calculate some things about the 11-Nov-2019 transit of Mercury across the Sun as seen from the city of Caracas, Venezuela.</p>

<pre><code>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 &lt; 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()
</code></pre>

CC BY-SA 4.0
2019-11-16T00:32:21.453Z
The reference to the JPL model is useful.  Point taken that "edge" is not clearly defined for the sun.  Mm
Automaton
CC BY-SA 4.0
2019-11-16T02:19:58.783Z
@Automaton This isn't an exhaustive list of uncertainties at all, but I think it is the largest one. Uncertainties in the positions of the Sun, Mercury and Earth are much smaller in comparison. See for example ftp://ssd.jpl.nasa.gov/pub/eph/planets/ioms/de423.iom.pdf and also [This orientation is the limiting error source for the orbits of the terrestrial planets, and corresponds to orbit uncertainties of a few hundred meters.](https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de430_and_de431.pdf) for example for DE431.
uhoh
CC BY-SA 4.0
2019-11-16T02:23:26.767Z
and also [Uncertainties in the JPL Planetary Ephemeris](https://syrte.obspm.fr/jsr/journees2010/pdf/Folkner.pdf)
uhoh
CC BY-SA 4.0
2019-11-28T02:34:32.437Z
The graphs in Uncertainties in the JPL Planetary Ephemeris are very interesting.  I'm familiar with the N body physics problem which is characterised as chaotic, so I was expecting uncertainties rapidly go up with time e.g. exponentially but it isn't the case in these graphs.
Automaton
CC BY-SA 4.0
2019-11-28T03:27:10.100Z
@Automaton orbits in the solar system are definitely chaotic, but it becomes a problem over millions of years or more, way beyond where these ephemerides are calculated. There is tightly constrained and accurate enough observational data over enough time now to predict with fairly good accuracy over tens of thousands of years. As a random  example, see the [arXiv](https://arxiv.org/abs/1802.04718) paper linked in [this answer](https://space.stackexchange.com/a/25385/12102) and search for "chaotic".
uhoh
CC BY-SA 4.0
2019-11-14T05:06:55.827Z
2019-11-14T05:20:39.203Z
2019-11-14T05:20:39.203Z
uhoh
<p>What is the error range for calculations for the time of 1st-4th contacts for a given GPS coordinate on Earth?</p>

<p>I observed through a telescope the 3rd and 4th contacts and compared against SkySafari and <a href="https://stellarium-web.org/" rel="nofollow noreferrer">Stellarium</a>.  I found SkySafari showed contacts about 10 seconds earlier for the same location vs <a href="https://stellarium-web.org/" rel="nofollow noreferrer">Stellarium</a>.  Stellarium was very close to my measurements for 3rd contact (1-2 seconds).  To get the calculation I zoomed the GUI and tweaked the time until it was closest to the contact.  Estimates for 4th contact (transit end) are much harder to visually time.</p>

<p>This got me interested to know if there is a gold standard for such calculations, and what is the range of error for it.   For forecasts n years out is the error range n times more?</p>

<p><a href="https://i.sstatic.net/E5ePS.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/E5ePS.png" alt="Stellarium Web visual calculation of 3rd contact time"></a></p>

CC BY-SA 4.0
2019-11-13T21:56:57.150Z
@barrycarter I tried https://wgc.jpl.nasa.gov:8443/webgeocalc/#OccultationFinder but couldn't figure it out.  For kernel selection I picked "Solar System Kernels" for lack of a better guess, and picked ellipsoid shapes for mercury and sun.  I got CSPICE_N0066: CSPICE.gfoclt: SPICE(INVALIDFRAME): [gfoclt_c --> GFOCLT --> GFOCCE --> ZZGFOCIN] Supposed body-fixed frame ITRF93 for front target 199 is actually centered on body 399.  So it's clearly not meant for amateurs like myself.  However, I see there are some tutorials (https://naif.jpl.nasa.gov/naif/tutorials.html) which may explain things.
Automaton
<planetary-transits>
Accuracy of Mercury transit calculations
Automaton
# Accuracy of Mercury transit calculations

(Asked by: username_0 on 2019-11-13)

What is the error range for calculations for the time of 1st-4th contacts for a given GPS coordinate on Earth?

I observed through a telescope the 3rd and 4th contacts and compared against SkySafari and [Stellarium](https://stellarium-web.org/). I found SkySafari showed contacts about 10 seconds earlier for the same location vs [Stellarium](https://stellarium-web.org/). Stellarium was very close to my measurements for 3rd contact (1-2 seconds). To get the calculation I zoomed the GUI and tweaked the time until it was closest to the contact. Estimates for 4th contact (transit end) are much harder to visually time.

This got me interested to know if there is a gold standard for such calculations, and what is the range of error for it. For forecasts n years out is the error range n times more?

[![Stellarium Web visual calculation of 3rd contact time](https://i.sstatic.net/E5ePS.png)](https://i.sstatic.net/E5ePS.png)

**Tags:** planetary-transits
> **Comments:**
> **username_0**: @barrycarter I tried https://wgc.jpl.nasa.gov:8443/webgeocalc/#OccultationFinder but couldn't figure it out.  For kernel selection I picked "Solar System Kernels" for lack of a better guess, and picked ellipsoid shapes for mercury and sun.  I got CSPICE_N0066: CSPICE.gfoclt: SPICE(INVALIDFRAME): [gfoclt_c --> GFOCLT --> GFOCCE --> ZZGFOCIN] Supposed body-fixed frame ITRF93 for front target 199 is actually centered on body 399.  So it's clearly not meant for amateurs like myself.  However, I see there are some tutorials (https://naif.jpl.nasa.gov/naif/tutorials.html) which may explain things.  

---

## 1 answers

### Answer by username_1 on 2019-11-14. Score: 3 | Accepted answer

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. 

How to pronounce “Ephemerides”?


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

> **Comments:**
> **username_0**: The reference to the JPL model is useful.  Point taken that "edge" is not clearly defined for the sun.  Mm  
> **username_1**: @Automaton This isn't an exhaustive list of uncertainties at all, but I think it is the largest one. Uncertainties in the positions of the Sun, Mercury and Earth are much smaller in comparison. See for example ftp://ssd.jpl.nasa.gov/pub/eph/planets/ioms/de423.iom.pdf and also [This orientation is the limiting error source for the orbits of the terrestrial planets, and corresponds to orbit uncertainties of a few hundred meters.](https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/de430_and_de431.pdf) for example for DE431.  
> **username_1**: and also [Uncertainties in the JPL Planetary Ephemeris](https://syrte.obspm.fr/jsr/journees2010/pdf/Folkner.pdf)  
> **username_0**: The graphs in Uncertainties in the JPL Planetary Ephemeris are very interesting.  I'm familiar with the N body physics problem which is characterised as chaotic, so I was expecting uncertainties rapidly go up with time e.g. exponentially but it isn't the case in these graphs.  
> **username_1**: @Automaton orbits in the solar system are definitely chaotic, but it becomes a problem over millions of years or more, way beyond where these ephemerides are calculated. There is tightly constrained and accurate enough observational data over enough time now to predict with fairly good accuracy over tens of thousands of years. As a random  example, see the [arXiv](https://arxiv.org/abs/1802.04718) paper linked in [this answer](https://space.stackexchange.com/a/25385/12102) and search for "chaotic".  

---

