I have an ephemeris SP3 file with coords of one satellite in ITRF with time in TAI. I’ve constructed form it an sp3 propagator like this (code in python but I accept responses in both Python and Java):
data = DataSource("ssaja320.b23130.e23140.DG_.sp3.001") sp3_file = sp3.SP3Parser().parse(data) sp3_eph = sp3_file.getSatellites().get("L39") sp3_prop = sp3_eph.getPropagator()
next I defined a little function that will help me convert python datetimes to AbsoluteDate:
def datetime_to_AbsoluteDate(row, colname): tmp_date = row[colname] y = tmp_date.year m = tmp_date.month d = tmp_date.day h = tmp_date.hour mi = tmp_date.minute s = tmp_date.second ms = float(tmp_date.microsecond/1e6) s = s + ms return AbsoluteDate(y,m,d,h,mi,s, TimeScalesFactory.getUTC())
I have a dataframe in which one of the column is named “UTC” and it contains datetimes of observations. With it I can construct new column with AbsoluteDates:
df["AbsDate"] = df.apply(lambda row: datetime_to_AbsoluteDate(row, "UTC"), axis = 1)
Then I define frame and observing station:
itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, True) earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, itrf) longitude = math.radians(lon) latitude = math.radians(lat) altitude = el station = GeodeticPoint(latitude, longitude, altitude) topo = TopocentricFrame(earth, station, "Observer")
Next I’m getting PV coordinates for my AbsDates and converting them to topocentric AZEL:
df["timestamped"] = df.apply(lambda row: sp3_prop.getPVCoordinates(row["AbsDate"], itrf), axis = 1) df["SP3_alt_noref"] = df.apply(lambda row: math.degrees(topo.getElevation(row["timestamped"].getPosition(), itrf, row["timestamped"].getDate())), axis = 1) df["SP3_az"] = df.apply(lambda row: math.degrees(topo.getAzimuth(row["timestamped"].getPosition(), itrf, row["timestamped"].getDate())), axis = 1)
Later I apply refraction to elevation.
I have a set of data generated with STK for comparison. In ITRF data from SP3 propagated with Orekit are consistent with data from STK (there are around 2-3 cm differences which is negligible for satellite of altitude ~1330km). Unfortunately data converted to AZEL using Orekit differs from data from STK. Residuals Orekit-STK look like this:
Does anybody has any idea if there is a mistake in my code? Parameters like latitude, longitude, elevation are the same in both STK and my code. Refraction model should be roughly the same as well, but even if there are small differences in parameters used for refraction it should only affect elevation, but the biggest differences are in azimuth.