<p>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:</p>
<pre><code>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')
</code></pre>
<p>The output:</p>
<pre><code>Aphelion:
2030-07-04 13:00:00 UTC
152099543.53658378 km

Perihelion:
2030-12-31 23:00:00 UTC
147101197.3727886 km
</code></pre>
<p>So what's going on?</p>
<p>The problem is that the perihelion of 2031 happens to be a much &quot;deeper&quot; perihelion than that of 2030. So, even though <code>2030-12-31 23:00:00 UTC</code> is still several days away from the 2031 perihelion, the Earth at that point is already so close to the sun that it &quot;wins&quot; the contest against the 2030 perihelion and gets returned instead.</p>
<p>What you really want for perihelion is &quot;the moment that is closer than <em>the two adjacent moments</em>&quot;, 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:</p>
<pre><code>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 &lt; dprev)
smaller_than_next = (d &lt; dnext)
both = (smaller_than_previous &amp; 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')
</code></pre>
<p>The result is the moment of perihelion:</p>
<pre><code>Perihelion:
2030-01-03 10:00:00 UTC
147105837.56789857 km
</code></pre>
<p>There might be even more elegant ways to do this in NumPy; a web search for 'numpy local minima' should suggest several.</p>
<p>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?</p>

CC BY-SA 4.0
2024-03-10T22:10:43.790Z
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
Klaus Rohe
CC BY-SA 4.0
2024-03-10T22:40:05.003Z
Note that Brandon is the author of Skyfield, so he knows what he's talking about. ;)
PM 2Ring
CC BY-SA 4.0
2024-03-11T12:33:24.543Z
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/
Klaus Rohe
CC BY-SA 4.0
2024-03-10T16:49:01.403Z
2024-03-10T16:49:01.403Z
Brandon Rhodes
<p>There is a nice article about finding local extrema in numpy arrays:
<a href="https://blog.finxter.com/how-to-find-local-minima-in-1d-and-2d-numpy-arrays/" rel="nofollow noreferrer">https://blog.finxter.com/how-to-find-local-minima-in-1d-and-2d-numpy-arrays/</a></p>

CC BY-SA 4.0
2024-03-11T12:36:42.467Z
2024-03-11T12:36:42.467Z
Klaus Rohe
<p>I'm trying to calculate the true distance between the sun and the Earth using skyfield.</p>
<p>To estimate the distance at a time instance I use the barycentric coordinates of the sun and earth using the following Python code:</p>
<pre><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
</code></pre>
<p>I then use the (xyz)-coordinates of the sun and earth to calculate the distance with the
formula</p>
<pre><code>sqrt((xE - xS)**2+(yE - yS)**2 + (zE - zS)**2) 
</code></pre>
<p>where <code>xE, xS, ... zE, zS</code> are the
coordinates of the sun and the earth.</p>
<p>Uisng this procedure I get for the maximum rep. minimum distance of the sun from the earth in 2030 the following values:</p>
<p>minimum 0.98331 AU at 31.12.2030<br />
maximum 1.01672 AU at 04.07.2030</p>
<p>The date for the minimum seems a little bit strange to me.</p>
<p>Is the above-described procedure to calculate the sun-earth distance correct or do I miss something?</p>

CC BY-SA 4.0
2024-03-10T22:47:05.017Z
The exact time (and distance) of Earth's perihelion is a little wobbly, mostly because of the influence of the Moon. See my answers https://astronomy.stackexchange.com/a/49546/16685 & https://astronomy.stackexchange.com/a/49605/16685
PM 2Ring
CC BY-SA 4.0
2024-03-10T23:27:28.733Z
FWIW, the 2031 perihelion is ~2031-Jan-04 20:50 UTC, at a distance of ~0.98326643 au. https://ssd.jpl.nasa.gov/api/horizons.api?format=text&COMMAND=10&CENTER=%40399&QUANTITIES=20&START_TIME=%272031-Jan-4+20%3A00+UT%27&STOP_TIME=%272031-Jan-4+22%3A00%27&STEP_SIZE=10m
PM 2Ring
CC BY-SA 4.0
2024-03-11T14:21:09.727Z
A script I threw together as an answer to https://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`
notovny
<python><ephemeris><skyfield><jpl>
How to calculate the true distance of the sun and the earth using skyfield?
Klaus Rohe
# How to calculate the true distance of the sun and the earth using skyfield?

(Asked by: username_0 on 2024-03-10)

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?

**Tags:** python, ephemeris, skyfield, jpl

> **Comments:**
> **username_1**: The exact time (and distance) of Earth's perihelion is a little wobbly, mostly because of the influence of the Moon. See my answers https://astronomy.stackexchange.com/a/49546/16685 & https://astronomy.stackexchange.com/a/49605/16685  
> **username_1**: FWIW, the 2031 perihelion is ~2031-Jan-04 20:50 UTC, at a distance of ~0.98326643 au. https://ssd.jpl.nasa.gov/api/horizons.api?format=text&COMMAND=10&CENTER=%40399&QUANTITIES=20&START_TIME=%272031-Jan-4+20%3A00+UT%27&STOP_TIME=%272031-Jan-4+22%3A00%27&STEP_SIZE=10m  
> **username_2**: A script I threw together as an answer to https://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`  

---

## 2 answers

### Answer by username_3 on 2024-03-10. Score: 10

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?

> **Comments:**
> **username_0**: Hi username_3, 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 username_0  
> **username_1**: Note that username_3 is the author of Skyfield, so he knows what he's talking about. ;)  
> **username_0**: 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/  

---

### Answer by username_0 on 2024-03-11. Score: 1

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


---

