# How to calculate the true distance of the sun and the earth using skyfield?

I'm trying to calculate the true distance between the sun and the Earth using skyfield.

To estimate the distance at a time instance I use the barycentric coordinates of the sun and earth using the following Python code:

eph = load('de440.bsp');
earth = eph['earth'];
sun = eph['sun'];
pos_earth = earth.at(t).position.au;
pos_sun = sun.at(t).position.au


I then use the (xyz)-coordinates of the sun and earth to calculate the distance with the formula

sqrt((xE - xS)**2+(yE - yS)**2 + (zE - zS)**2)


where xE, xS, ... zE, zS are the coordinates of the sun and the earth.

Uisng this procedure I get for the maximum rep. minimum distance of the sun from the earth in 2030 the following values:

minimum 0.98331 AU at 31.12.2030
maximum 1.01672 AU at 04.07.2030

The date for the minimum seems a little bit strange to me.

Is the above-described procedure to calculate the sun-earth distance correct or do I miss something?

• The exact time (and distance) of Earth's perihelion is a little wobbly, mostly because of the influence of the Moon. See my answers astronomy.stackexchange.com/a/49546/16685 & astronomy.stackexchange.com/a/49605/16685 Commented Mar 10 at 22:47
• FWIW, the 2031 perihelion is ~2031-Jan-04 20:50 UTC, at a distance of ~0.98326643 au. ssd.jpl.nasa.gov/api/… Commented Mar 10 at 23:27
• A script I threw together as an answer to astronomy.stackexchange.com/a/55631/23813 uses skyfield's find_minima and find_maxima functions to hunt down apoapses and periapses. It pulls 2030-01-03 10:12:26 +0000 and 2030-07-04 12:57:38 +0000, when set to search 2030. but it's also using de421.bsp Commented Mar 11 at 14:21

Your procedure is sound! Here is the same calculation, using Skyfield's ability to build a vector of all the hours of the year 2030, and NumPy's ability to find where the maximum and minimum occur in a sequence — and it gives the same odd answer as your script:

import numpy as np
from skyfield.api import load

t = ts.utc(2030, 1, 1, range(365 * 24))
sun = eph['sun']
earth = eph['earth']
difference = (earth - sun).at(t).distance().km

i = np.argmax(difference)
j = np.argmin(difference)

print('Aphelion:')
print(t[i].utc_strftime())
print(difference[i], 'km')

print()

print('Perihelion:')
print(t[j].utc_strftime())
print(difference[j], 'km')


The output:

Aphelion:
2030-07-04 13:00:00 UTC
152099543.53658378 km

Perihelion:
2030-12-31 23:00:00 UTC
147101197.3727886 km


So what's going on?

The problem is that the perihelion of 2031 happens to be a much "deeper" perihelion than that of 2030. So, even though 2030-12-31 23:00:00 UTC is still several days away from the 2031 perihelion, the Earth at that point is already so close to the sun that it "wins" the contest against the 2030 perihelion and gets returned instead.

What you really want for perihelion is "the moment that is closer than the two adjacent moments", so that we rule out returning a point that's at the very end of our range. We want only a point where the Sun distance is less than both of the points around it. Here's one way NumPy can do that:

import numpy as np
from skyfield.api import load

t = ts.utc(2030, 1, 1, range(365 * 24))
sun = eph['sun']
earth = eph['earth']
difference = (earth - sun).at(t).distance().km

dprev = difference[:-2]
d = difference[1:-1]
dnext = difference[2:]
smaller_than_previous = (d < dprev)
smaller_than_next = (d < dnext)
both = (smaller_than_previous & smaller_than_next)
i = np.argmax(both) + 1  # the '1' corrects for our trimming the array

print('Perihelion:')
print(t[i].utc_strftime())
print(difference[i], 'km')


The result is the moment of perihelion:

Perihelion:
2030-01-03 10:00:00 UTC
147105837.56789857 km


There might be even more elegant ways to do this in NumPy; a web search for 'numpy local minima' should suggest several.

Skyfield also has maxima and minima search functions, if you want to try them out. Maybe I should add built-in perihelion and aphelion routines someday?

• Hi Brandon, thank you very much for the help. Because the Earth orbit is near circluar it is a challenge to get the Aphelion and Perihelion. May be it also helps to evalutate the velocity. Kind Regards Klaus Commented Mar 10 at 22:10
• Note that Brandon is the author of Skyfield, so he knows what he's talking about. ;) Commented Mar 10 at 22:40
• There is a nice article about finding local extrema in numpy arrays: blog.finxter.com/… Commented Mar 11 at 12:33

There is a nice article about finding local extrema in numpy arrays: https://blog.finxter.com/how-to-find-local-minima-in-1d-and-2d-numpy-arrays/