In :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 :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: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 :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)
The computation below should give something close to zero, as we are approximating the position at
t1 by the position plus velocity at
In :rv1.position.km - rv0.position.km - rv0.velocity.km_per_s
Out: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 :rv1.position.km - rv0.position.km
If we ask the satellite to compute its position in ITRF instead of GCRS, then the problem doesn't happen.
In :rv0_ITRF = sat.ITRF_position_velocity_error(t0)[:2] r0_ITRF = rv0_ITRF * AU_KM v0_ITRF = rv0_ITRF * AU_KM / DAY_S rv1_ITRF = sat.ITRF_position_velocity_error(t1)[:2] r1_ITRF = rv1_ITRF * AU_KM v1_ITRF = rv1_ITRF * AU_KM / DAY_S
The same calculation for the Z component now gives a much smaller error.
In :r1_ITRF - r0_ITRF - v0_ITRF
Out: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 :r1_ITRF - r0_ITRF
So it seems there is a problem with how the velocity is handled by
Update 2019-07-24: This shows how this bug affects Doppler computations.
In :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 :np.sum(los0.velocity.km_per_s * los0.position.km) / length_of(los0.position.km)
In :np.sum(los1.velocity.km_per_s * los1.position.km) / length_of(los1.position.km)
If instead we compute the velocity vector subtracting the positions at t1 and t0:
In :v = los1.position.km - los0.position.km np.sum(v * los0.position.km) / length_of(los0.position.km)
There is a huge difference (a factor of 2) between both ways of computing the Doppler. The second way gives the correct result.