```
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]:
```

```
In [3]:
```ts = skyfield.api.load.timescale()
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)

`t1`

by the position plus velocity at `t0`

.

```
In [4]:
```rv1.position.km - rv0.position.km - rv0.velocity.km_per_s

```
Out[4]:
```

`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]:
```

```
In [6]:
```rv0.velocity.km_per_s[2]

```
Out[6]:
```

```
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]:
```

`t0`

, we see that they are very close.

```
In [9]:
```r1_ITRF[2] - r0_ITRF[2]

```
Out[9]:
```

```
In [10]:
```v0_ITRF[2]

```
Out[10]:
```

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]:
```

```
In [13]:
```np.sum(los1.velocity.km_per_s * los1.position.km) / length_of(los1.position.km)

```
Out[13]:
```

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]:
```