Hey guys,
I’m working on my first python script and having some issues implementing different propagators. I set up three different propagators (Keplerian, Numerical with DormandPrice843 and SGP4). Therefore I tried to build a TLE by using Keplerian parameters. Keplerian and Numerical Propagator seem to fit very well, as they’re always on a similar attitude & position. But when I choose the SGP4 Propagator, something is not working out as intended.
For example when I shift the time 100.0s from the starting point I get the following value differences, using numerical as reference:
attitude +2km for SGP4
position x +5km
position y -38km
position z +32km
Maybe the TLE is wrong? Maybe there’s an issue when converting from TEME to GCRF frame? I didn’t manage to find out…
Thank you for helping me out!
Code:
frame = FramesFactory.getGCRF()
#date of mission launch (reference)
year = 2026
month = 10
day = 1
hour= 0
minute = 0
second = 0.0
utc = TimeScalesFactory.getUTC()
initDate = AbsoluteDate(year,month,day,hour,minute,second,utc)
t_shifted = 100.0 #elapsed seconds
finalDate = initDate.shiftedBy(t_shifted)
#Phase 1 - SSO Data - Circular - CircularOrbit because of performance reasons
a_c = r_earth + 450e3 #semi-major axis
i_c = FastMath.toRadians(98.6) #inclination
raan_c = FastMath.toRadians(0.0) #argument of ascending node
omega_c = FastMath.toRadians(0.0) #argument of perigee
nu_c = FastMath.toRadians(0.0) #true anomaly
u_c = omega_c+nu_c #argument of latitude
#create Orbit und propagate, circular because of performance reasons
create_orbit = CircularOrbit(a_c,0.0,0.0,i_c,raan_c,u_c,PositionAngle.TRUE,
FramesFactory.getGCRF(),initDate,mu_earth)
new_orbit = create_orbit.shiftedBy(t_shifted)
new_orbit_K = OrbitType.KEPLERIAN.convertType(new_orbit) #convert circular to keplerian for better understanding
#choose Propagator
strategy = int(input(“Choose Propagator and enter number\nKeplerian: 1\nNumerical: 2\nSGP4: 3”))
if strategy == 1:
propagator = KeplerianPropagator(create_orbit)
propagate_orbit = propagator.propagate(finalDate)
elif strategy == 2:
minStep = 0.001
maxStep = 1000.0
initStep = 60.0
positionTolerance = 1.0
tolerances = NumericalPropagator.tolerances(positionTolerance, create_orbit, create_orbit.getType())
#tolerance array: position 0 allowed absolute error, position 1 allowed relative error
integrator = DormandPrince853Integrator(minStep, maxStep, JArray_double.cast_(tolerances[0]),
JArray_double.cast_(tolerances[1]))
#DormanPrince is a 8th-order integrator with variable step size
integrator.setInitialStepSize(initStep)
satellite_mass = 82.5 #model needs a spacecraft mass [kg]
initialState = SpacecraftState(create_orbit, satellite_mass)
propagator_num = NumericalPropagator(integrator)
propagator_num.setOrbitType(OrbitType.CARTESIAN)
propagator_num.setInitialState(initialState)
gravityProvider = GravityFieldFactory.getNormalizedProvider(10, 10)
propagator_num.addForceModel(
HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, True), gravityProvider))
orbit = OrbitType.KEPLERIAN.convertType(propagator_num.getInitialState().getOrbit())
propagate_orbit = propagator_num.propagate(initDate, finalDate)
elif strategy == 3:
tle_line1 = “1 99999U 26001A 26274.00000000 .00000000 00000-0 00000-0 0 9991”
tle_line2 = “2 99999 98.6000 0.0000 0000000 00.0000 0.0000 15.38683946 01”
my_TLE= TLE(tle_line1, tle_line2)
print(my_TLE)
propagate_orbit = TLEPropagator.selectExtrapolator(my_TLE)
if strategy in [1, 2]: # SpacecraftState
position = propagate_orbit.getPVCoordinates(frame).getPosition().scalarMultiply(1/1000)
velocity = propagate_orbit.getPVCoordinates(frame).getVelocity().scalarMultiply(1/1000)
acceleration = propagate_orbit.getPVCoordinates(frame).getAcceleration()
elif strategy == 3: # TLEPropagator
position = propagate_orbit.getPVCoordinates(finalDate, frame).getPosition().scalarMultiply(1/1000)
velocity = propagate_orbit.getPVCoordinates(finalDate, frame).getVelocity().scalarMultiply(1/1000)
acceleration = propagate_orbit.getPVCoordinates(finalDate, frame).getAcceleration()
#calculate absolute values
r = math.sqrt(position.getX()**2 + position.getY()**2 + position.getZ()**2)
v_abs = math.sqrt(velocity.getX()**2 + velocity.getY()**2 + velocity.getZ()**2)
a_abs = math.sqrt(acceleration.getX()**2 + acceleration.getY()**2 + acceleration.getZ()**2)