Hello!
I have been fiddling with the Python wrapper of Orekit for a little while now, specifically with propagating GPS satellites from SEM almanac files and integrating the code with FastAPI and Cesium JS. Recently, I got my hands on an STK license and decided to use it to validate my Orekit implementation, for some reassurance that I’m invoking Orekit methods correctly. Initially, I used both STK and Orekit to propagate the 31 GPS sats from this almanac file:
current_sem.al3 (6.7 KB) for a month and record the ephemeris at every time step (arbitrarily chose 10 minutes/600 seconds). The ephemeris from each are almost always within 1 meter (usually better) of each other, even out to a month in the future, as I was hoping/expecting.
Now, I tried to take it a step further and have both programs yield the line-of-sight rise/set times of each of the 31 GPS sats, from the perspective of a point I chose on Earth; [lat: 27.9517, lon: -82.4588, alt: 0.] which is roughly located in Tampa, Florida. Knowing that the propagation is at least in agreement, I am now seeing major discrepancies between Orekit and STK with these “access” times.
This is my current Orekit implementation (forgive the lack of includes, it’s part of a much larger and messier main.py):
app = FastAPI()
vm = orekit.initVM()
from orekit.pyhelpers import setup_orekit_curdir
setup_orekit_curdir()
from org.orekit.data import DataContext, DirectoryCrawler
DATA_CONTEXT = DataContext.getDefault()
orekit_data_path = '/src'
orekit_data_manager = DATA_CONTEXT.getDataProvidersManager()
from java.io import File
orekit_data_manager.addProvider(DirectoryCrawler(File(orekit_data_path)))
from org.orekit.gnss import SEMParser
sp = SEMParser(None)
sp.loadData()
from org.orekit.propagation.analytical.gnss.data import GPSAlmanac
almanacs = sp.getAlmanacs()
app.almanacs = almanacs
class EarthLocation(BaseModel):
lat: float
lon: float
alt: float
# Dictionary for TimeScale names to singleton TimeScale instance
TimeScaleByName = {
'TAI': TimeScalesFactory.getTAI(),
'UTC': TimeScalesFactory.getUTC(),
'GPS': TimeScalesFactory.getGPS(),
'TT': TimeScalesFactory.getTT()
}
def getVisTestParallel(isoStrInitialEpoch: str, isoStrFinalEpoch: str, maskRadians: float, points: list[EarthLocation]):
orekit.getVMEnv().attachCurrentThread()
# Define the earth's shape.
equatorialRadius = Constants.WGS84_EARTH_EQUATORIAL_RADIUS
flattening = Constants.WGS84_EARTH_FLATTENING
earthFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
earth = OneAxisEllipsoid(equatorialRadius, flattening, earthFrame)
initialEpoch = AbsoluteDate(isoStrInitialEpoch, TimeScaleByName["GPS"])
finalEpoch = AbsoluteDate(isoStrFinalEpoch, TimeScaleByName["GPS"]) # Timescale should be GPS, n'est-ce pas?
propagators = []
eventsLoggers = []
for alm in app.almanacs:
logger = EventsLogger() # Add new event to be detected
prop = GNSSPropagatorBuilder(alm).build() # Build GNSS propagator from individual almanac entry
i = 0
for point in points:
i += 1
stationPOV = GeodeticPoint(point.lat, point.lon, point.alt)
prn = alm.getPRN()
svn = alm.getSVN()
name = str(prn) + ":" + str(svn) + ":" + str(i)
stationFrame = TopocentricFrame(earth, stationPOV, name) # set name to PRN:SVN:point index
# Define the elevation event checker above a masking angle
stationVis = ElevationDetector(stationFrame).withConstantElevation(maskRadians).withHandler(ContinueOnEvent()) # maskingAngle is user input [rad]
# The detector should only be added to the propagator wrapped in the logger
prop.addEventDetector(logger.monitorDetector(stationVis))
propagators.append(prop)
eventsLoggers.append(logger)
step_handler = SimpleStepHandler()
parallelizer = PropagatorsParallelizer(Arrays.asList(propagators), step_handler)
parallelizer.propagate(initialEpoch, finalEpoch)
# For now we're returning the results as a comma-separated set of strings like "2:61:1,2023-09-24T17:23:21.15050731898612Z,True",
# where the form is (PRN):(SVN):(Point index), (event time), (rising/setting, where true = rising).
results = []
for eLogger in eventsLoggers:
allEvents = eLogger.getLoggedEvents()
# riseTimes = []
# setTimes = []
for logged in allEvents:
elDetector = ElevationDetector.cast_(logged.getEventDetector())
name = elDetector.getTopocentricFrame().getName()
# # eventState = logged.getState().toString()
eventTime = logged.getDate().toString()
print(name, eventTime, logged.isIncreasing())
results.append("%s,%s,%s" % (name, eventTime, logged.isIncreasing()))
return results
Here are the results for the first two weeks or so from GPS sat PRN 02/SVN 6, propagated from the epoch of the almanac included above, 2024-08-23T19:56:30.000Z, where the first column denotes a “rise” time (logged.increasing == “true”), and the second column denotes a “set” time (logged.increasing == “false”), I have changed the format a tad and truncated decimals from what’s exported by the above code to match what comes out of the STK access report:
23 Aug 2024 23:07:57.504
24 Aug 2024 14:44:35.534 24 Aug 2024 23:03:55.026
25 Aug 2024 14:40:31.737 25 Aug 2024 22:59:52.548
26 Aug 2024 14:36:27.940 26 Aug 2024 22:55:50.070
27 Aug 2024 14:32:24.143 27 Aug 2024 22:51:47.591
28 Aug 2024 14:28:20.346 28 Aug 2024 22:47:45.113
29 Aug 2024 14:24:16.549 29 Aug 2024 22:43:42.634
30 Aug 2024 14:20:12.752 30 Aug 2024 22:39:40.154
31 Aug 2024 14:16:08.955 31 Aug 2024 22:35:37.675
1 Sep 2024 14:12:05.158 1 Sep 2024 22:31:35.195
2 Sep 2024 14:08:01.362 2 Sep 2024 22:27:32.716
3 Sep 2024 14:03:57.565 3 Sep 2024 22:23:30.236
4 Sep 2024 13:59:53.768 4 Sep 2024 22:19:27.755
5 Sep 2024 13:55:49.971 5 Sep 2024 22:15:25.275
6 Sep 2024 13:51:46.174 6 Sep 2024 22:11:22.794
7 Sep 2024 13:47:42.377 7 Sep 2024 22:07:20.313
8 Sep 2024 13:43:38.580 8 Sep 2024 22:03:17.832
And here are the results from the STK scenario I ran, which I am more than happy to provide any additional context or details as to how I built the scenario, if that is the suggested suspect (but I think it’s my Orekit implementation):
23 Aug 2024 20:10:31.282 23 Aug 2024 23:47:33.332
24 Aug 2024 02:58:34.015 24 Aug 2024 09:03:59.863
24 Aug 2024 20:06:27.685 24 Aug 2024 23:43:28.122
25 Aug 2024 02:54:29.858 25 Aug 2024 08:59:56.320
25 Aug 2024 20:02:24.084 25 Aug 2024 23:39:22.910
26 Aug 2024 02:50:25.708 26 Aug 2024 08:55:52.777
26 Aug 2024 19:58:20.490 26 Aug 2024 23:35:17.696
27 Aug 2024 02:46:21.551 27 Aug 2024 08:51:49.233
27 Aug 2024 19:54:16.891 27 Aug 2024 23:31:12.483
28 Aug 2024 02:42:17.395 28 Aug 2024 08:47:45.686
28 Aug 2024 19:50:13.296 28 Aug 2024 23:27:07.265
29 Aug 2024 02:38:13.244 29 Aug 2024 08:43:42.146
29 Aug 2024 19:46:09.698 29 Aug 2024 23:23:02.055
30 Aug 2024 02:34:09.087 30 Aug 2024 08:39:38.603
30 Aug 2024 19:42:06.098 30 Aug 2024 23:18:56.843
31 Aug 2024 02:30:04.937 31 Aug 2024 08:35:35.060
31 Aug 2024 19:38:02.503 31 Aug 2024 23:14:51.628
1 Sep 2024 02:26:00.780 1 Sep 2024 08:31:31.517
1 Sep 2024 19:33:58.905 1 Sep 2024 23:10:46.415
2 Sep 2024 02:21:56.626 2 Sep 2024 08:27:27.970
2 Sep 2024 19:29:55.311 2 Sep 2024 23:06:41.201
3 Sep 2024 02:17:52.474 3 Sep 2024 08:23:24.430
3 Sep 2024 19:25:51.713 3 Sep 2024 23:02:35.985
4 Sep 2024 02:13:48.319 4 Sep 2024 08:19:20.887
4 Sep 2024 19:21:48.113 4 Sep 2024 22:58:30.771
5 Sep 2024 02:09:44.168 5 Sep 2024 08:15:17.344
5 Sep 2024 19:17:44.519 5 Sep 2024 22:54:25.557
6 Sep 2024 02:05:40.011 6 Sep 2024 08:11:13.801
6 Sep 2024 19:13:40.921 6 Sep 2024 22:50:20.344
7 Sep 2024 02:01:35.861 7 Sep 2024 08:07:10.254
7 Sep 2024 19:09:37.326 7 Sep 2024 22:46:15.129
8 Sep 2024 01:57:31.705 8 Sep 2024 08:03:06.715
8 Sep 2024 19:05:33.729 8 Sep 2024 22:42:09.912
Again, more than happy and able to provide more sample data (either longer periods of time or the results for more of the GPS sats), just don’t want to flood this post too much with numbers.
The first thing I noticed was that the STK results yielded an [expected] second window of visibility/access per day, whereas Orekit only yielded the one. On the “set” times that do align, Orekit is consistently about 40 minutes behind STK. The cadence of the Orekit times matches that of every-other STK time. All of these things lead me to believe I have an incorrect frame implementation somewhere in my Orekit code. I have double checked the EOP used in both situations, and they do match.
Is there something obvious that is just not jumping out at me? Or does anybody have any suggestions for further investigation?