``````

In [1]:

import skyfield.api
from skyfield.sgp4lib import EarthSatellite
from skyfield.constants import AU_KM, DAY_S
from skyfield.functions import length_of
import numpy as np

``````

Update 2020-07-17: this has been re-run with the last version of Skyfield, where the problem has been fixed (see #224). I have mantained all the comments about the problems that Skyfield had, so now they don't match the results.

As a test TLE we use a TLE for Es'hail, which is in a nearly GEO orbit at 24ºE.

``````

In [2]:

sat = EarthSatellite('1 43700U 18090A   18335.89431171 +.00000133 +00000-0 +00000-0 0  9993',\
'2 43700 000.0858 245.4352 0001094 006.6237 164.6135 01.00274015000309')

sat.epoch.utc

``````
``````

Out[2]:

CalendarTuple(year=2018, month=12, day=1, hour=21, minute=27, second=48.5317263007164)

``````

The TLE epoch is near 2018-12-01 21 UTC. For simplicity, we compute the position and velocity of the satellite at 21:00:00 UTC and 21:00:01 UTC.

``````

In [3]:

t0 = ts.utc(2018, 12, 1, 21, 0, 0)
t1 = ts.utc(2018, 12, 1, 21, 0, 1)
rv0 = sat.at(t0)
rv1 = sat.at(t1)

``````

The computation below should give something close to zero, as we are approximating the position at `t1` by the position plus velocity at `t0`.

``````

In [4]:

rv1.position.km - rv0.position.km - rv0.velocity.km_per_s

``````
``````

Out[4]:

array([-1.08401336e-05, -1.40169930e-04,  2.14398141e-06])

``````

The Z component should be much smaller. Indeed, let us compare the Z component of the velocity estimated by taking the position at `t1` minus the position at `t0`, with the velocity computed at `t0`. We see that they differ by an order of magnitude.

``````

In [5]:

rv1.position.km[2] - rv0.position.km[2]

``````
``````

Out[5]:

-0.00026495564615558465

``````
``````

In [6]:

rv0.velocity.km_per_s[2]

``````
``````

Out[6]:

-0.0002670996275668849

``````

If we ask the satellite to compute its position in ITRF instead of GCRS, then the problem doesn't happen.

``````

In [7]:

rv0_ITRF = sat.ITRF_position_velocity_error(t0)[:2]
r0_ITRF = rv0_ITRF[0] * AU_KM
v0_ITRF = rv0_ITRF[1] * AU_KM / DAY_S
rv1_ITRF = sat.ITRF_position_velocity_error(t1)[:2]
r1_ITRF = rv1_ITRF[0] * AU_KM
v1_ITRF = rv1_ITRF[1] * AU_KM / DAY_S

``````

The same calculation for the Z component now gives a much smaller error.

``````

In [8]:

r1_ITRF - r0_ITRF - v0_ITRF

``````
``````

Out[8]:

array([ 2.56276996e-05, -6.11424304e-05,  2.16614035e-06])

``````

If we compare the Z components of the velocity estimated by taking differences or by getting the velocity at `t0`, we see that they are very close.

``````

In [9]:

r1_ITRF[2] - r0_ITRF[2]

``````
``````

Out[9]:

-0.00453944263882633

``````
``````

In [10]:

v0_ITRF[2]

``````
``````

Out[10]:

-0.00454160877917917

``````

So it seems there is a problem with how the velocity is handled by `ITRF_to_GCRS2()`.

Update 2019-07-24: This shows how this bug affects Doppler computations.

``````

In [11]:

observer  = skyfield.api.Topos(latitude = 40.595865, longitude = -3.699069, elevation_m = 800)
los0 = (sat - observer).at(t0)
los1 = (sat - observer).at(t1)

``````

Doppler computed as projection of velocity vector onto line of sight vector:

``````

In [12]:

np.sum(los0.velocity.km_per_s * los0.position.km) / length_of(los0.position.km)

``````
``````

Out[12]:

0.0005891289567179068

``````
``````

In [13]:

np.sum(los1.velocity.km_per_s * los1.position.km) / length_of(los1.position.km)

``````
``````

Out[13]:

0.0005891225228701842

``````

If instead we compute the velocity vector subtracting the positions at t1 and t0:

``````

In [14]:

v = los1.position.km - los0.position.km
np.sum(v * los0.position.km) / length_of(los0.position.km)

``````
``````

Out[14]:

0.00048325541071942595

``````

There is a huge difference (a factor of 2) between both ways of computing the Doppler. The second way gives the correct result.