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

Question

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?

Answer 1

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
ts = load.timescale()
t = ts.utc(2030, 1, 1, range(365 * 24))
eph = load('de440.bsp')
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
ts = load.timescale()
t = ts.utc(2030, 1, 1, range(365 * 24))
eph = load('de440.bsp')
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?

Answer 2

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/