Difference of ICRF position flips velocity sign

[shalen-original]

Hi :)

I think I found a bug. Take this example:

from datetime import datetime
from skyfield.api import Topos, load
from skyfield.sgp4lib import EarthSatellite

loc = Topos(latitude_degrees=45, longitude_degrees=45, elevation_m=10)
sat = EarthSatellite(name="NOAA 9", line1="1 15427U 84123A   20101.49917043  .00000013  00000-0  29655-4 0  9993", line2="2 15427  98.9347  55.5631 0013793 271.1939 225.9918 14.16094114823900")
t = load.timescale().utc(datetime.fromisoformat("2020-04-10T15:00:00.000000+00:00"))

v = (sat.at(t) - loc.at(t)).velocity.km_per_s
v_diff = sat.at(t).velocity.km_per_s - loc.at(t).velocity.km_per_s

print("diff.velocity.km_per_s = {}".format(v))
print("v_diff                 = {}".format(v_diff))

It outputs:

diff.velocity.km_per_s = [ 2.76918265  5.43919105 -3.85544706]
v_diff                 = [-2.76918265 -5.43919105  3.85544706]

After a bit of investigation I reached this line in the ICRF class:

def __sub__(self, body):
        p = self.position.au - body.position.au
        if self.velocity is None or body.velocity is None:
            v = None
        else:
            v = body.velocity.au_per_d - self.velocity.au_per_d  # <------ HERE
        return ICRF(p, v, self.t)

Shouldn't the two terms be flipped? That is, self.velocity - body.velocity, as done for the position component above?

If it is actually a bug and you want, I can put together a PR with the fix and a test to check it.

[brandon-rhodes]

Thanks very much for letting me know! Yes, it's a bug, and I've just added a fix that should appear in the next version of Skyfield.

[shalen-original]

Great, thank you :)