Calculation speed of skyfield orbit propagation for multiple EarthSatellite objects loaded from TLE

[dsholes]

Apologies if this is too similar to #30 and #183. I decided to investigate the bottleneck for orbit propagation (position and velocity) for multiple skyfield.sgp4lib.EarthSatellite objects, using EarthSatellite.at()

Summary:

• Bottleneck is related to skyfield.positionlib.ITRF_to_GCRS2(), which relies on "expensive" functions from skyfield.nutationlib and skyfield.precessionlib to calculate important attributes for skyfield.timelib.Time object.

  • If future work related to Numba-jit takes place (mentioned in Issue #30), it would be nice to address skyfield.nutationlib and skyfield.precessionlib

• In the meantime, for those who are interested in potential speed boosts, I have two ideas for reducing computation time for multiple satellite orbits

  • Option 1 - Ideally, we would use the same skyfield.timelib.Time object for every skyfield.sgp4lib.EarthSatellite object that we want to propagate. Then we could pre-compute skyfield.timelib.Time.MT and skyfield.timelib.Time.gast (the "expensive" attributes mentioned above), and the results would be cached for all of the proceeding satellite propagations.
    • Since accuracy of the sgp4 propagation rapidly degrades before/after TLE epoch, users must understand the risks associated with using Times "far away" from TLE epoch for each satellite, since TLE epoch is different for each satellite in the Celestrak NORAD elements (e.g. 'http://celestrak.com/NORAD/elements/planet.txt')
    • For 15 satellites and 5760 orbit positions, using a single time array (skyfield.timelib.Time object) for all satellites and pre-computing/caching t.MT and t.gast reduces computation time from 35 sec to 5 sec.
  • Option 2 - Another idea is to allow users (with a keyword argument) to return position and velocity in TEME or ITRF reference frames instead of GCRS, for reduced computation time. This allows for even greater speed (~3 sec for 15 satellites and 5760 orbit positions), and doesn't restrict user to single time array (skyfield.timelib.Time object) for all satellites. The obvious drawback is some users need pos and vel in GCRS

Orbit propagation, 1 satellite, 5760 orbit positions

```python from skyfield.api import Topos, load import numpy as np import time ``` ```python # Orbit propagation parameters, for 5760 orbit positions over time SIM_TIME = 1.6*60.*60. # in seconds TIME_STEP = 1. # in seconds ``` ### Total computation time for orbit propagation, 1 satellite, 5760 orbit positions ```python start_propagation = time.time() ##### Load TLE ts = load.timescale() stations_url = 'http://celestrak.com/NORAD/elements/planet.txt' satellites = load.tle(stations_url,reload=False,) cubesat = satellites['FLOCK 2K-08'] ##### ##### Build time array cubesat_tle_epoch = cubesat.epoch.utc_datetime() seconds = cubesat_tle_epoch.second + np.arange(0, SIM_TIME, TIME_STEP) # time_array = ts.utc(year = cubesat_tle_epoch.year, month = cubesat_tle_epoch.month, day = cubesat_tle_epoch.day, hour = cubesat_tle_epoch.hour, minute = cubesat_tle_epoch.minute, second = seconds) ##### ##### Perform EarthSatellite.at() orbit propagation cubesat_propagated = cubesat.at(time_array) ##### ##### Get GCRS Position in km position_gcrs = cubesat_propagated.position.km ##### ##### Get GCRS Velocity in km/s velocity_gcrs = cubesat_propagated.velocity.km_per_s ##### stop_propagation = time.time() print("--- Skyfield Orbit Propagation ---\n"\ "Total Computation Time (for 1 sat, {0} sequential orbit positions): {1} seconds" .format(int(SIM_TIME/TIME_STEP), stop_propagation-start_propagation)) ``` --- Skyfield Orbit Propagation --- Total Computation Time (for 1 sat, 5760 sequential orbit positions): 2.5103211402893066 seconds

Orbit propagation, 15 satellites, 5760 orbit positions

### Total computation time for orbit propagation, 15 satellites, 5760 orbit positions NOTE: Reinitializing `time_array` for each new `EarthSatellite` in `for` loop ```python start_propagation = time.time() ##### Load TLE ts = load.timescale() stations_url = 'http://celestrak.com/NORAD/elements/planet.txt' satellites = load.tle(stations_url,reload=False,) flock_2k_keys = [k for k in satellites.keys() if "FLOCK 2K" in str(k)] flock_2k_keys.sort() for satkey in flock_2k_keys[0:15]: cubesat = satellites[satkey] ##### Build time array cubesat_tle_epoch = cubesat.epoch.utc_datetime() seconds = cubesat_tle_epoch.second + np.arange(0, SIM_TIME, TIME_STEP) # time_array = ts.utc(year = cubesat_tle_epoch.year, month = cubesat_tle_epoch.month, day = cubesat_tle_epoch.day, hour = cubesat_tle_epoch.hour, minute = cubesat_tle_epoch.minute, second = seconds) ##### ##### Perform EarthSatellite.at() orbit propagation cubesat_propagated = cubesat.at(time_array) ##### ##### Get GCRS Position in km position_gcrs = cubesat_propagated.position.km ##### ##### Get GCRS Velocity in km/s velocity_gcrs = cubesat_propagated.velocity.km_per_s ##### stop_propagation = time.time() print("--- Skyfield Orbit Propagation ---\n"\ "Total Computation Time (for 15 sats, {0} sequential orbit positions): {1} seconds" .format(int(SIM_TIME/TIME_STEP), stop_propagation-start_propagation)) ``` --- Skyfield Orbit Propagation --- Total Computation Time (for 15 sats, 5760 sequential orbit positions): 35.4270076751709 seconds

Speed Up: OPTION 1

# Option 1 (scroll up to see computation time without caching) ### Total computation time for orbit propagation, 15 satellites, 5760 orbit positions, pre-computing attributes for `time_array` and caching for successive runs ```python start_propagation = time.time() ##### Load TLE ts = load.timescale() stations_url = 'http://celestrak.com/NORAD/elements/planet.txt' satellites = load.tle(stations_url,reload=False,) flock_2k_keys = [k for k in satellites.keys() if "FLOCK 2K" in str(k)] flock_2k_keys.sort() cubesat = satellites[flock_2k_keys[0]] ##### ##### Build time array cubesat_tle_epoch = cubesat.epoch.utc_datetime() seconds = cubesat_tle_epoch.second + np.arange(0, SIM_TIME, TIME_STEP) time_array = ts.utc(year = cubesat_tle_epoch.year, month = cubesat_tle_epoch.month, day = cubesat_tle_epoch.day, hour = cubesat_tle_epoch.hour, minute = cubesat_tle_epoch.minute, second = seconds) ##### ##### Pre-compute expensive attributes start_mt_gast = time.time() time_array.MT time_array.gast stop_mt_gast = time.time() ##### for satkey in flock_2k_keys[0:15]: cubesat = satellites[satkey] ##### Perform EarthSatellite.at() orbit propagation cubesat_propagated = cubesat.at(time_array) ##### ##### Get GCRS Position in km position_gcrs = cubesat_propagated.position.km ##### ##### Get GCRS Velocity in km/s velocity_gcrs = cubesat_propagated.velocity.km_per_s ##### stop_propagation = time.time() print("--- Skyfield Orbit Propagation ---\n"\ "Total Computation Time (for 15 sats, {0} sequential orbit positions): {1} seconds" .format(int(SIM_TIME/TIME_STEP), stop_propagation-start_propagation)) print("(NOTE: computation time to pre-compute t.MT and t.gast = {0} seconds)".format(stop_mt_gast-start_mt_gast)) ``` --- Skyfield Orbit Propagation --- Total Computation Time (for 15 sats, 5760 sequential orbit positions): 4.998683214187622 seconds (NOTE: computation time to pre-compute t.MT and t.gast = 2.1935174465179443 seconds)

Speed Up: OPTION 2

# Option 2 ### Returning pos and vel in TEME and ITRF only NOTE: leads to further speed up, and allows for unique `time_array` for each satellite (no caching needed). Obviously not an option if user needs pos and vel in GCRS ```python from sgp4.propagation import sgp4 from skyfield.sgp4lib import TEME_to_ITRF,theta_GMST1982 from skyfield.positionlib import ITRF_to_GCRS2, ITRF_to_GCRS from skyfield.constants import AU_KM, DAY_S, T0, tau ``` ```python def at_teme_and_itrf_only(t,sat_model): ##### sgp4 loop min_past_epoch = (t._utc_float() - sat_model.jdsatepoch)*1440. position_teme = [] velocity_teme = [] for i,m in enumerate(min_past_epoch): p, v = sgp4(sat_model,m) position_teme.append(p) velocity_teme.append(v) position_teme_km = np.array(position_teme).T velocity_teme_km_per_sec = np.array(velocity_teme).T ##### position_teme_au = position_teme_km/AU_KM velocity_teme_au_per_day = velocity_teme_km_per_sec/AU_KM velocity_teme_au_per_day *= DAY_S ##### TEME_to_ITRF Transformation position_itrf_au,velocity_itrf_au_per_day = TEME_to_ITRF(t.ut1,position_teme_au,velocity_teme_au_per_day) ##### position_itrf_km = position_itrf_au*AU_KM velocity_itrf_km_per_sec = velocity_itrf_au_per_day*AU_KM velocity_itrf_km_per_sec /= DAY_S return (position_teme_km,velocity_teme_km_per_sec,position_itrf_km,velocity_itrf_km_per_sec) ``` ```python start_propagation = time.time() ##### Load TLE ts = load.timescale() stations_url = 'http://celestrak.com/NORAD/elements/planet.txt' satellites = load.tle(stations_url,reload=False,) flock_2k_keys = [k for k in satellites.keys() if "FLOCK 2K" in str(k)] flock_2k_keys.sort() for satkey in flock_2k_keys[0:15]: cubesat = satellites[satkey] ##### Build time array cubesat_tle_epoch = cubesat.epoch.utc_datetime() seconds = cubesat_tle_epoch.second + np.arange(0, SIM_TIME, TIME_STEP) # time_array = ts.utc(year = cubesat_tle_epoch.year, month = cubesat_tle_epoch.month, day = cubesat_tle_epoch.day, hour = cubesat_tle_epoch.hour, minute = cubesat_tle_epoch.minute, second = seconds) ##### ##### Perform EarthSatellite.at() orbit propagation position_teme,velocity_teme,position_itrf,velocity_itrf = at_teme_and_itrf_only(time_array,cubesat.model) ##### stop_propagation = time.time() print("--- Skyfield Orbit Propagation ---\n"\ "Total Computation Time (for 15 sats, {0} sequential orbit positions): {1} seconds" .format(int(SIM_TIME/TIME_STEP), stop_propagation-start_propagation)) ``` --- Skyfield Orbit Propagation --- Total Computation Time (for 15 sats, 5760 sequential orbit positions): 2.73648738861084 seconds

Detailed Profiling

# --- DETAILED ANALYSIS --- ### Computation time for different components of orbit propagation ```python ##### Load TLE start_loadtle = time.time() ts = load.timescale() stations_url = 'http://celestrak.com/NORAD/elements/planet.txt' satellites = load.tle(stations_url,reload=False,) cubesat = satellites['FLOCK 2K-08'] stop_loadtle = time.time() ##### ##### Build time array start_buildtimearray = time.time() cubesat_tle_epoch = cubesat.epoch.utc_datetime() seconds = cubesat_tle_epoch.second + np.arange(0, SIM_TIME, TIME_STEP) time_array = ts.utc(year = cubesat_tle_epoch.year, month = cubesat_tle_epoch.month, day = cubesat_tle_epoch.day, hour = cubesat_tle_epoch.hour, minute = cubesat_tle_epoch.minute, second = seconds) stop_buildtimearray = time.time() ##### ##### Perform EarthSatellite.at() orbit propagation start_earthsatobj_at = time.time() cubesat_propagated = cubesat.at(time_array) stop_earthsatobj_at = time.time() ##### ##### Get GCRS Position in km start_geocentricobj_pos = time.time() position_gcrs = cubesat_propagated.position.km stop_geocentricobj_pos = time.time() ##### ##### Get GCRS Velocity in km/s start_geocentricobj_vel = time.time() velocity_gcrs = cubesat_propagated.velocity.km_per_s stop_geocentricobj_vel = time.time() ##### print("Load TLE: {0} seconds".format(stop_loadtle-start_loadtle)) print("Build time array: {0} seconds".format(stop_buildtimearray-start_buildtimearray)) print("Perform EarthSatellite.at() orbit propagation: {0} seconds".format(stop_earthsatobj_at-start_earthsatobj_at)) print("Get GCRS Position in km: {0} seconds".format(stop_geocentricobj_pos-start_geocentricobj_pos)) print("Get GCRS Velocity in km/s: {0} seconds".format(stop_geocentricobj_vel-start_geocentricobj_vel)) print("") ``` Load TLE: 0.028937339782714844 seconds Build time array: 0.0010006427764892578 seconds Perform EarthSatellite.at() orbit propagation: 2.3472838401794434 seconds Get GCRS Position in km: 0.0 seconds Get GCRS Velocity in km/s: 0.0 seconds ### Separate *Build time array* from `%%timeit` cell NOTE: The big difference between the first run and the rest of the runs indicates some attribute related to the `skyfield.timelib.Time` is being cached. ```python ##### Load TLE ts = load.timescale() stations_url = 'http://celestrak.com/NORAD/elements/planet.txt' satellites = load.tle(stations_url,reload=False,) cubesat = satellites['FLOCK 2K-08'] ##### ##### Build time array cubesat_tle_epoch = cubesat.epoch.utc_datetime() seconds = cubesat_tle_epoch.second + np.arange(0, SIM_TIME, TIME_STEP) time_array = ts.utc(year = cubesat_tle_epoch.year, month = cubesat_tle_epoch.month, day = cubesat_tle_epoch.day, hour = cubesat_tle_epoch.hour, minute = cubesat_tle_epoch.minute, second = seconds) ##### ``` ```python %%timeit -n 1 ##### Perform EarthSatellite.at() orbit propagation start_earthsatobj_at = time.time() cubesat_propagated = cubesat.at(time_array) stop_earthsatobj_at = time.time() ##### ##### Get GCRS Position in km start_geocentricobj_pos = time.time() position_gcrs = cubesat_propagated.position.km stop_geocentricobj_pos = time.time() ##### ##### Get GCRS Velocity in km/s start_geocentricobj_vel = time.time() velocity_gcrs = cubesat_propagated.velocity.km_per_s stop_geocentricobj_vel = time.time() ##### #print("Load TLE: {0} seconds".format(stop_loadtle-start_loadtle)) #print("Build time array: {0} seconds".format(stop_buildtimearray-start_buildtimearray)) print("Perform EarthSatellite.at() orbit propagation: {0} seconds".format(stop_earthsatobj_at-start_earthsatobj_at)) print("Get GCRS Position in km: {0} seconds".format(stop_geocentricobj_pos-start_geocentricobj_pos)) print("Get GCRS Velocity in km/s: {0} seconds".format(stop_geocentricobj_vel-start_geocentricobj_vel)) print("") ``` Perform EarthSatellite.at() orbit propagation: 2.3373708724975586 seconds Get GCRS Position in km: 0.0 seconds Get GCRS Velocity in km/s: 0.0 seconds Perform EarthSatellite.at() orbit propagation: 0.18749761581420898 seconds Get GCRS Position in km: 0.0 seconds Get GCRS Velocity in km/s: 0.0 seconds Perform EarthSatellite.at() orbit propagation: 0.1791830062866211 seconds Get GCRS Position in km: 0.0 seconds Get GCRS Velocity in km/s: 0.0 seconds Perform EarthSatellite.at() orbit propagation: 0.19774484634399414 seconds Get GCRS Position in km: 0.0 seconds Get GCRS Velocity in km/s: 0.0 seconds Perform EarthSatellite.at() orbit propagation: 0.18579626083374023 seconds Get GCRS Position in km: 0.0 seconds Get GCRS Velocity in km/s: 0.0 seconds Perform EarthSatellite.at() orbit propagation: 0.1874992847442627 seconds Get GCRS Position in km: 0.0 seconds Get GCRS Velocity in km/s: 0.0 seconds Perform EarthSatellite.at() orbit propagation: 0.19822025299072266 seconds Get GCRS Position in km: 0.0 seconds Get GCRS Velocity in km/s: 0.0 seconds The slowest run took 12.52 times longer than the fastest. This could mean that an intermediate result is being cached. 496 ms ± 752 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) ### Test components of *EarthSatellite.at()* NOTE: The bulk of the computation time occurs in the `skyfield.positionlib.ITRF_to_GCRS2` coordinate transformation function. ```python %%timeit -n 1 ##### Load TLE and build time array ts = load.timescale() stations_url = 'http://celestrak.com/NORAD/elements/planet.txt' satellites = load.tle(stations_url,reload=False,) cubesat = satellites['FLOCK 2K-08'] cubesat_tle_epoch = cubesat.epoch.utc_datetime() seconds = cubesat_tle_epoch.second + np.arange(0, SIM_TIME, TIME_STEP) time_array = ts.utc(year = cubesat_tle_epoch.year, month = cubesat_tle_epoch.month, day = cubesat_tle_epoch.day, hour = cubesat_tle_epoch.hour, minute = cubesat_tle_epoch.minute, second = seconds) ##### ##### sgp4 loop start_sgp4loop = time.time() min_past_epoch = time_array._utc_float() - cubesat.model.jdsatepoch position_teme = [] velocity_teme = [] for i,m in enumerate(min_past_epoch): p, v = sgp4(cubesat.model,m) position_teme.append(p) velocity_teme.append(v) position_teme = np.array(position_teme).T velocity_teme = np.array(velocity_teme).T stop_sgp4loop = time.time() ##### ##### TEME_to_ITRF Transformation start_itrf_calcs = time.time() position_itrf,velocity_itrf = TEME_to_ITRF(time_array.ut1,position_teme,velocity_teme) stop_itrf_calcs = time.time() ##### ##### ITRF_to_GCRS2 calcs start_gcrs_calcs = time.time() position_gcrs,velocity_gcrs = ITRF_to_GCRS2(time_array,position_itrf,velocity_itrf) stop_gcrs_calcs = time.time() ##### print("sgp4 loop: {0} seconds".format(stop_sgp4loop-start_sgp4loop)) print("TEME_to_ITRF Transformation: {0} seconds".format(stop_itrf_calcs-start_itrf_calcs)) print("ITRF_to_GCRS2 Transformation: {0} seconds".format(stop_gcrs_calcs-start_gcrs_calcs)) print("") ``` sgp4 loop: 0.17188477516174316 seconds TEME_to_ITRF Transformation: 0.0 seconds ITRF_to_GCRS2 Transformation: 2.143871545791626 seconds sgp4 loop: 0.1718745231628418 seconds TEME_to_ITRF Transformation: 0.0 seconds ITRF_to_GCRS2 Transformation: 2.1249587535858154 seconds sgp4 loop: 0.17186927795410156 seconds TEME_to_ITRF Transformation: 0.0 seconds ITRF_to_GCRS2 Transformation: 2.188218832015991 seconds sgp4 loop: 0.18648433685302734 seconds TEME_to_ITRF Transformation: 0.0 seconds ITRF_to_GCRS2 Transformation: 2.140382766723633 seconds sgp4 loop: 0.17915558815002441 seconds TEME_to_ITRF Transformation: 0.0 seconds ITRF_to_GCRS2 Transformation: 2.2860541343688965 seconds sgp4 loop: 0.20894432067871094 seconds TEME_to_ITRF Transformation: 0.0020017623901367188 seconds ITRF_to_GCRS2 Transformation: 2.264150857925415 seconds sgp4 loop: 0.1843414306640625 seconds TEME_to_ITRF Transformation: 0.002002239227294922 seconds ITRF_to_GCRS2 Transformation: 2.21419095993042 seconds 2.41 s ± 66 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) ### Separate *Build Time Array* from %%timeit cell NOTE: Caching of `skyfield.timelib.Time` is related to `skyfield.positionlib.ITRF_to_GCRS2` ```python ##### Load TLE ts = load.timescale() stations_url = 'http://celestrak.com/NORAD/elements/planet.txt' satellites = load.tle(stations_url,reload=False,) cubesat = satellites['FLOCK 2K-08'] ##### ##### Build time array cubesat_tle_epoch = cubesat.epoch.utc_datetime() seconds = cubesat_tle_epoch.second + np.arange(0, SIM_TIME, TIME_STEP) time_array = ts.utc(year = cubesat_tle_epoch.year, month = cubesat_tle_epoch.month, day = cubesat_tle_epoch.day, hour = cubesat_tle_epoch.hour, minute = cubesat_tle_epoch.minute, second = seconds) ##### ``` ```python %%timeit -n 1 ##### sgp4 loop start_sgp4loop = time.time() min_past_epoch = time_array._utc_float() - cubesat.model.jdsatepoch position_teme = [] velocity_teme = [] for i,m in enumerate(min_past_epoch): p, v = sgp4(cubesat.model,m) position_teme.append(p) velocity_teme.append(v) position_teme = np.array(position_teme).T velocity_teme = np.array(velocity_teme).T stop_sgp4loop = time.time() ##### ##### TEME_to_ITRF Transformation start_itrf_calcs = time.time() position_itrf,velocity_itrf = TEME_to_ITRF(time_array.ut1,position_teme,velocity_teme) stop_itrf_calcs = time.time() ##### ##### ITRF_to_GCRS2 calcs start_gcrs_calcs = time.time() position_gcrs,velocity_gcrs = ITRF_to_GCRS2(time_array,position_itrf,velocity_itrf) stop_gcrs_calcs = time.time() ##### print("sgp4 loop: {0} seconds".format(stop_sgp4loop-start_sgp4loop)) print("TEME_to_ITRF Transformation: {0} seconds".format(stop_itrf_calcs-start_itrf_calcs)) print("ITRF_to_GCRS2 Transformation: {0} seconds".format(stop_gcrs_calcs-start_gcrs_calcs)) print("") ``` sgp4 loop: 0.17531275749206543 seconds TEME_to_ITRF Transformation: 0.01563429832458496 seconds ITRF_to_GCRS2 Transformation: 2.188530683517456 seconds sgp4 loop: 0.18713903427124023 seconds TEME_to_ITRF Transformation: 0.0010008811950683594 seconds ITRF_to_GCRS2 Transformation: 0.0 seconds sgp4 loop: 0.1857895851135254 seconds TEME_to_ITRF Transformation: 0.002001047134399414 seconds ITRF_to_GCRS2 Transformation: 0.0 seconds sgp4 loop: 0.18268346786499023 seconds TEME_to_ITRF Transformation: 0.0 seconds ITRF_to_GCRS2 Transformation: 0.0 seconds sgp4 loop: 0.19740915298461914 seconds TEME_to_ITRF Transformation: 0.0010008811950683594 seconds ITRF_to_GCRS2 Transformation: 0.0010006427764892578 seconds sgp4 loop: 0.18395042419433594 seconds TEME_to_ITRF Transformation: 0.0010008811950683594 seconds ITRF_to_GCRS2 Transformation: 0.0010006427764892578 seconds sgp4 loop: 0.1882617473602295 seconds TEME_to_ITRF Transformation: 0.0010001659393310547 seconds ITRF_to_GCRS2 Transformation: 0.0 seconds The slowest run took 12.81 times longer than the fastest. This could mean that an intermediate result is being cached. 502 ms ± 767 ms per loop (mean ± std. dev. of 7 runs, 1 loop each) ### Testing `skyfield.timelib.Time.MT` and `skyfield.timelib.Time.gast` NOTE: Method `skyfield.positionlib.ITRF_to_GCRS2()` calls `t.gast` and `t.MT` in the spin and position calculations. A closer look at `skyfield.timelib` shows `@reify` decorators used for caching. Methods `t.gast` and `t.MT` both call "expensive" functions from `skyfield.nutationlib` and `skyfield.precessionlib`, such as `skyfield.nutationlib.compute_nutation()`, `skyfield.precessionlib.compute_precession()` and `skyfield.nutationlib.earth_tilt()`. ```python ##### Load TLE ts = load.timescale() stations_url = 'http://celestrak.com/NORAD/elements/planet.txt' satellites = load.tle(stations_url,reload=False,) cubesat = satellites['FLOCK 2K-08'] ##### ##### Build time array cubesat_tle_epoch = cubesat.epoch.utc_datetime() seconds = cubesat_tle_epoch.second + np.arange(0, SIM_TIME, TIME_STEP) time_array = ts.utc(year = cubesat_tle_epoch.year, month = cubesat_tle_epoch.month, day = cubesat_tle_epoch.day, hour = cubesat_tle_epoch.hour, minute = cubesat_tle_epoch.minute, second = seconds) ##### ``` ```python %%timeit -n 1 ##### Run skyfield.timelib.Time.MT start_time_MT = time.time() time_array.MT stop_time_MT = time.time() ##### ##### Run skyfield.timelib.Time.gast start_time_gast = time.time() time_array.gast stop_time_gast = time.time() ##### print("skyfield.timelib.Time.MT Computation: {0} seconds".format(stop_time_MT-start_time_MT)) print("skyfield.timelib.Time.gast Computation: {0} seconds".format(stop_time_gast-start_time_gast)) print("") ``` skyfield.timelib.Time.MT Computation: 1.102233648300171 seconds skyfield.timelib.Time.gast Computation: 1.0522608757019043 seconds skyfield.timelib.Time.MT Computation: 0.0 seconds skyfield.timelib.Time.gast Computation: 0.0 seconds skyfield.timelib.Time.MT Computation: 0.0 seconds skyfield.timelib.Time.gast Computation: 0.0 seconds skyfield.timelib.Time.MT Computation: 0.0 seconds skyfield.timelib.Time.gast Computation: 0.0 seconds skyfield.timelib.Time.MT Computation: 0.0 seconds skyfield.timelib.Time.gast Computation: 0.0 seconds skyfield.timelib.Time.MT Computation: 0.0 seconds skyfield.timelib.Time.gast Computation: 0.0 seconds skyfield.timelib.Time.MT Computation: 0.0 seconds skyfield.timelib.Time.gast Computation: 0.0 seconds The slowest run took 79284.54 times longer than the fastest. This could mean that an intermediate result is being cached. 308 ms ± 754 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

[brandon-rhodes]

Here's a further idea: satellite positions are inherently uncertain enough that it might not make any practical difference if you had Skyfield skip the expensive — but very slight — corrections for nutation and precession.

What if instead of waiting until the time object computes the expensive matrix t.M you instead set it to the identity matrix? Like:

import numpy as np
t.M = np.identity(3)

See if that speeds things up and, if it does, see how many meters or kilometers difference it makes in the satellite positions that come out.

[dsholes]

Great idea, I'll take a look and let you know. By the way, I imagine that will only address the position part of skyfield.positionlib.ITRF_to_GCRS2(), but the spin calculation takes around the same amount of time (it relies on t.gast). Any similar ideas for t.gast?

Thank you for your quick responses (and of course all of your wonderful work thus far on skyfield!)

[brandon-rhodes]

Good question! I don't have any immediate ideas about gast, since knowing which way the Earth is pointing is an unavoidable part of the calculation. I wonder if a lower-precision guess for gast could be provided instead of the full precision result?

[dsholes]

I wanted to provide a more detailed response, but haven't had the time the last couple weeks. Just wanted to share the basic idea in case someone else wants to investigate further. From some initial tests, it looks like changing t.MT leads to bigger differences in position that's returned than changing t.gast to t.gmst. Both changes lead to around the same speed up in computation time. Using t.gmst instead of t.gast in the spin calculation of positionlib.ITRF_to_GCRS2() (i.e. spin = rot_z(t.gast * tau / 24.0)) leads to speed up, and very small change in the final calculated position.

In summary, don't use t.M = np.identity(3), but instead use t.gmst in place of t.gast. I need to run some more tests to confirm, but perhaps this will get the ball rolling for someone else.

[brandon-rhodes]

Thanks for the report! When you say "leads to bigger differences", do you mean differences in computation time, or differences in the position that's returned?

[dsholes]

Sorry for the confusion, I've updated my original comment.

[dsholes]

I noticed that t.M changes very little, at least over the time scales that I'm interested in (and for which the TLEs are valid). Perhaps instead of t.M = np.identity(3), a better guess is t.M = t[0].M. As long as t spans less than a few days, then I think the error will be very low. Needs more thorough testing, but for now I'll pre-compute t.gast = t.gmst and t.M = t[0].M whenever I want a speed boost. Thanks again for the ideas.

[dsholes]

A somewhat different approach:

I noticed that sgp4.propagation.sgp4() only accepted scalars. So anytime .at(TimeArray) is called, skyfield must loop over all times in the TimeArray. For-loops in Python are obviously much slower than vectorized calculations when math heavy computations are being performed (as in sgp4). I've started to "vectorize" the sgp4 code, which is mostly just switching the math module functions for numpy functions, but there are a number of if-then statements that will need to be updated as well.

I've changed enough to run simple tests on LEO satellites (Flock). I've submitted a pull request here that outlines the changes, and any insight is appreciated.

[brandon-rhodes]

I have just released a new version of Skyfield that (a) uses a new version of sgp4 that compiles the official C++ code for speedy computation of positions, and (b) that can search for satellite passes fairly quickly:

https://rhodesmill.org/skyfield/earth-satellites.html#finding-when-a-satellite-rises-and-sets

I found (at least if I measured correctly!) that if I set both t.M to the identity matrix and t.gast to 0.0, the resulting azimuth and altitude were unchanged — because those two values cancel out completely (I think?) when Skyfield goes from an Earth location, out to ICRF coordinates, and then back in to the position of the satellite above the Earth.

Oh, and, the new sgp4 library if used raw has the ability to take whole arrays of satellites, so n dates × m satellites can be computed in a single call.

If future work related to Numba-jit takes place (mentioned in Issue #30), it would be nice to address skyfield.nutationlib and skyfield.precessionlib

I should answer this comment from the original issue, which I'm not sure I ever addressed: I suspect the expense of those routines is already almost entirely in NumPy's C code. I don't think there's much that Numba could do because so little of the work is in Python?

We should put together a list of requirements for getting this issue closed. I'm imagining a new section of the Earth Satellites documentation page, that at least shares the hints and tips we have for speeding up computation when accuracy can be sacrificed.