Inter-satellites visibility is not commutative

[aballet]

Hello,

I am working with inter-satellites visibilities, i.e. when satellite1 could have a line of sight towards satellite2 without being blocked by the Earth. I noticed that sometimes depending on the order of the vector subtraction (satellite1 - satellite2, or satellite2 - satellite1), we get different results.

I think this could be a bug. I cannot find a reasonable explanation on why the order of operations matters: if satellite1 sees satellite2, then the opposite should be true as well.

Here's a working piece of code, with the results below. This should be run in Python 3.10. For the figure to be displayed, the package kaleido must be installed as well.

# Imports
from skyfield.api import EarthSatellite
from datetime import datetime, timedelta, timezone
import plotly.express as px
import pandas as pd
import numpy as np

# Utility method to generate datetimes evenly-spaced
def steps(start: datetime, end: datetime, step: timedelta):
    while start < end:
        yield start
        start += step

# Returns a numpy array of booleans representing if the satellite1 has a line of sight towards satellite2 at each datetime in the list
def get_visibilities(satellite1: EarthSatellite, satellite2: EarthSatellite, datetimes=list[datetime]):
    times = satellite1.ts.from_datetimes(datetimes)
    vector = satellite1 - satellite2
    return ~vector.at(times).is_behind_earth()

# Returns a list of time intervals when satellite1 sees satellite2
def get_visibility_intervals(satellite1: EarthSatellite, satellite2: EarthSatellite, datetimes=list[datetime]):
    visibilities = get_visibilities(satellite1, satellite2, datetimes=datetimes)
    indices = np.where(visibilities[:-1] != visibilities[1:])[0]

    if visibilities[0]:
        indices = np.insert(indices, 0, 0)

    if visibilities[-1]:
        indices = np.insert(indices, len(indices), len(visibilities)-1)

    for start, end in zip(indices[::2], indices[1::2]):
        yield datetimes[start], datetimes[end]

# Load two satellites
satellite1 = EarthSatellite(
    "1 39418U 13066C   23052.81509870  .00005076  00000+0  39672-3 0  9998",
    "2 39418  97.5384 122.0838 0033543 136.4642 223.9239 15.01076184506049",
    name="SKYSAT-A",
)
satellite2 = EarthSatellite(
    "1 45788U 20038BL  23052.87050788  .00050160  00000+0  70337-3 0  9994",
    "2 45788  52.9867 191.0024 0003500 313.9003  46.1717 15.56275066153575",
    name="SKYSAT-C14"
)

# Propagate for a few hours
start = datetime(2023, 1, 1, 10, tzinfo=timezone.utc)
end = start + timedelta(hours=4)
step = timedelta(seconds=10)

# Compute the visibilities in both directions
# We swap the satellites here
datetimes = list(steps(start, end, step))
intervals1 = list(get_visibility_intervals(satellite1, satellite2, datetimes))
intervals2 = list(get_visibility_intervals(satellite2, satellite1, datetimes))
print(f"{satellite1.name} --> {satellite2.name}: {len(intervals1)} visibility intervals.")
print(f"{satellite2.name} --> {satellite1.name}: {len(intervals2)} visibility intervals.")

data = []
for s, e in intervals1:
    data.append(dict(start=s, end=e, name=f"{satellite1.name} --> {satellite2.name}"))
for s, e in intervals2:
    data.append(dict(start=s, end=e, name=f"{satellite2.name} --> {satellite1.name}"))
df = pd.DataFrame(data)

# Plot the data
fig = px.timeline(df, x_start="start", x_end="end", y="name")
fig.write_image("visibilities.png")
fig.show()

Below, we notice a gap of visibility of two minutes in one direction. Depending on the pair of satellites and the time horizon, the duration of the gap can vary.

[brandon-rhodes]

Drat—I definitely agree that visibility should be commutative, so something must be wrong here.

Thanks so much for providing a script to get an investigation started.

I have a busy next few days, but hopefully I'll have some time to sit down, run the script and confirm that I get the same results on my own computer, and then take a look to see what the math is doing behind the scenes. I'll report back once I've found something! (Though any other investigators are welcome to weigh in in the meantime.)

[Bernmeister]

@aballet I added

#!/usr/bin/env python3

to the top and installed via pip3 the module plotly and then ran with the result:

$ python3 test.py 
Traceback (most recent call last):
  File "test.py", line 18, in <module>
    def get_visibilities(satellite1: EarthSatellite, satellite2: EarthSatellite, datetimes=list[datetime]):
TypeError: 'type' object is not subscriptable

[aballet]

Oh right, I should have mentioned that this was written in Python 3.10. You can remove all the type annotations, they have no impact on the results.

Sorry about that.

[aballet]

To see the plot in an easy way, you can also run this in a Jupyter notebook. If you use python3 ..., the figure will not be saved. I will provide the edited piece of code in an hour to fix this.

[aballet]

Hi, I updated the code in my first message. The problem is now fully reproducible with hard-coded TLEs. To display the figure, the package kaleido must be installed as well. If you don't want to install anything else, just delete the last 3 lines of code, the problem will show up in the standard output as well.

[brandon-rhodes]

@aballet — I had a bit of time this evening, so I was able to sit down and run your script, and it turns out that the is_behind_earth() method wasn't properly accounting for the situation where the Earth is on the line joining the two satellite positions, but where the Earth is actually out past the other satellite rather than sitting in between them. So I've edited the routine to properly recognize when the Earth's disc is in the background behind the other satellite, rather than in front of it.

The improvement should come out in the next version of Skyfield. But if you want to try it now, install the development version of this project:

pip install -U https://github.com/skyfielders/python-skyfield/archive/master.zip

Thanks again for letting me know about the problem!

[aballet]

And thank you for taking the time to investigate and fix it!