Hi everyone!
I have been trying to build a software to track a satellite from the ground and I am trying to debug the propagation aspect. The input that I have is a TLE retrieved from space-track so I’ve used the code below the generate my satellite’s position for the on-going 2 weeks.
I thought I would use MATLAB to compare the propagations and I was actually expecting to obtain the exact same results - my thought was that MATLAB aerospace library was actually based on Orekit. But I get some deviations.
On Matlab I create my satellite using the following command and inputing the same TLE:
%% Simulation timeline
SimuTime=86400; % s
startTime = datetime(2025,7,28,4,0,0);
sampleTime = 1;
stopTime = startTime+seconds(SimuTime);
scenario = satelliteScenario(startTime,stopTime,sampleTime);
%% TLE Spacetrack
mySat = satellite(scenario,"mySat_TLE.tle");
For instance, when computing my 1st visibility I get:
- with Orekit:
- Start visibility: 2025-07-28T14:10:24
- End visibility: 2025-07-28T14:13:29
- Duration: 186s
- with Matlab:
- Start visibility: 2025-07-28T14:10:31
- End visibility: 2025-07-28T14:13:24
- Duration: 173s
It is obviously not a huge difference but I wouldn’t want to miss our satellite by only seconds
Would any of you be able to explain the differences that I observe?
Thank you very much for your help!
Benoist
Orekit code
TLE initialTLE = new TLE("1 99999U 25135AZ 25209.16473940 .00000882 00000-0 88936-4 0 9991",
"2 99999 97.7524 322.9374 0011509 294.1619 65.8397 14.93607525 5545");
AbsoluteDate initialDate = initialTLE.getDate();
Frame inertialFrame = FramesFactory.getEME2000();
// Define Earth shape
CelestialBody earth = CelestialBodyFactory.getEarth();
Frame earthFrame = FramesFactory.getITRF(IERSConventions.IERS_2010, true);
BodyShape earthShape = new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, earthFrame);
double simDuration = 14*86400; // 2 weeks in seconds
double dt = 1; // 1s
TLEPropagator tlePropagator = TLEPropagator.selectExtrapolator(initialTLE);
tlePropagator.setAttitudeProvider(new CelestialBodyPointed(inertialFrame, earth, Vector3D.PLUS_K, Vector3D.PLUS_K, Vector3D.MINUS_J));
tlePropagator.propagate(initialDate.shiftedBy(simDuration));
final EphemerisGenerator generator = tlePropagator.getEphemerisGenerator();
BoundedPropagator ephemeris = generator.getGeneratedEphemeris();
// Define ground station
double latitude = 43.6047; // Toulouse latitude in degrees
double longitude = 1.4442; // Toulouse longitude in degrees
double altitude = 150; // Altitude in meters
GeodeticPoint stationPoint = new GeodeticPoint(FastMath.toRadians(latitude), FastMath.toRadians(longitude), altitude);
TopocentricFrame stationFrame = new TopocentricFrame(earthShape, stationPoint, "Toulouse");
double visibilityTime = 0.0;
boolean visible = false;
int visibilityCount = 0;
for (AbsoluteDate propagationDate = initialDate;
propagationDate.compareTo(initialDate.shiftedBy(simDuration)) <= 0;
propagationDate = propagationDate.shiftedBy(dt)){
SpacecraftState spacecraftState = ephemeris.propagate(propagationDate);
// Propagate the orbit to the current date
PVCoordinates pvInert = spacecraftState.getPVCoordinates();
Attitude attitude = spacecraftState.getAttitude();
// Transform the satellite position to the station frame
PVCoordinates pvStationFrame = inertialFrame.getTransformTo(stationFrame, propagationDate).transformPVCoordinates(pvInert);
// Compute geodetic coordinates
GeodeticPoint satelliteGeodetic = earthShape.transform(spacecraftState.getPVCoordinates().getPosition(), inertialFrame, spacecraftState.getDate());
double lat = satelliteGeodetic.getLatitude();
double lon = satelliteGeodetic.getLongitude();
// Compute azimuth & elevation
double azimuth = stationFrame.getAzimuth(pvStationFrame.getPosition(), stationFrame, propagationDate);
double elevation = stationFrame.getElevation(pvStationFrame.getPosition(), stationFrame, propagationDate);
// Convert UTC to local time at the ground station location
ZonedDateTime localDateTime = ZonedDateTime.ofInstant(propagationDate.toInstant(), ZoneId.of("Europe/Paris"));
AbsoluteDate localDate = new AbsoluteDate( localDateTime.toLocalDate().getYear(),
localDateTime.toLocalDate().getMonthValue(),
localDateTime.toLocalDate().getDayOfMonth(),
localDateTime.toLocalTime().getHour(),
localDateTime.toLocalTime().getMinute(),
localDateTime.toLocalTime().getSecond(), TimeScalesFactory.getUTC());
if (Math.toDegrees(elevation) >= 25) {
visibilityTime += dt;
if (!visible) {
visible = true;
visibilityCount += 1;
}
} else {
visibilityTime = 0.0;
if (visible) {
visible = false;
}
}
}