Skyfield: achieve sgp4 results with 1 second periodicity for given time interval

Question

I've implemented sgp4 algorithm using Skyfield in python.

I need the position vector for 1 day time interval with 1 second periodicity. For that, I have to calculate the sgp4 repeatedly, adding 1 second each time.

Is there a way to return a vector, which contains all position?

Answer 1

You can do this by creating a skyfield.EarthSatellite object form a TLE, and then generating astrometric objects for each time in a list of times. You need to specify your list of times as timezone-aware datetime objects. You can create your list of times with something like (probably lots of room to make this more efficient):

from datetime import timezone, timedelta, datetime
from skyfield.api import load

ts = load.timescale()  # create skyfield timescale object
tz = timezone(timedelta(hours=-4))  # whatever your timezone offset from UTC is
start = datetime(2018, 06, 21, 0, 0, 0, tzinfo=tz)  # timezone-aware start time
end = start + timedelta(hours=24)  # one day's worth of times
delta = timedelta(minutes=1)  # your interval over which you evaluate
times = [start]
now = start
while now <= end:
    now += delta
    times.append(now)

And evaluate it with:

astrometrics = my_earthsat.at(ts.utc(times))

This will give you a list of astrometrics at each time that you specified. You can use these objects to get positions (and I think velocities) of your object in just about any units or reference frame you like. You can get, for instance, the altitude and azimuth in radians for each element in the list with something like:

alt_az = []
for ast in astrometrics:
    alt_az.append(ast.altaz().radians)

Answer 2

Happily, Skyfield is designed for time objects that are NumPy arrays and will generate as many positions as there are objects in the array. Here's a section of the documentation describing the feature:

http://rhodesmill.org/skyfield/time.html#date-arrays

And here's an example that accomplishes the task you're describing:

from numpy import arange
from skyfield.api import load, EarthSatellite

iss_tle0 = """\
1 25544U 98067A   18184.80969102  .00001614  00000-0  31745-4 0  9993
2 25544  51.6414 295.8524 0003435 262.6267 204.2868 15.54005638121106
"""
my_earthsat = EarthSatellite(*iss_tle0.splitlines())

tz_offset = -4

ts = load.timescale()
t = ts.utc(2018, 6, 21, tz_offset, 0, arange(24 * 60 * 60))
astrometrics = my_earthsat.at(t)

This code isn't particularly fast at the moment because Earth satellite computations have yet to be optimized, but should be faster than constructing multiple time objects manually in a loop.