Accurate prediction of relative formation position in LEO

Orekit Python Wrapper usage
1

Dear all,
I am currently trying to propagate the trajectories of two different satellites in LEO (500km), which are characterized by the same initial orbit Keplerian elements, except for the true anomaly - that is, a separation in along-track of a few hundred meters is imposed.

I am currently using this method to evaluate the relative position vector of one spacecraft with respect to the other, computed in one of the spacecraft LVLH frame. I propagate the orbits of the two spacecrafts independently, and then use the method linked to get the relative position of one spacecraft with respect to the LVLH of the other at each step. The propagation duration is set at one orbital period.

I have obtained different results when adopting for propagation:

  • a Keplerian propagator
  • an Eckestein-Hecler propagator
  • a numerical propagator (DormandPrince853 including HolmesFeatherstone model up to 10x10, and optionally drag and sun/moon attraction, min/maxSteps set at 0.001/2 seconds, position tolerance 1E-3).

I have also noticed very different trends between the relative position evolution (in LVLH, so along/across/radial directions) when I set the orbit eccentricity at 0 instead of the nominal 0.001, and I do not understand if that is correct or not.

My question is whether which result may be the most representative of reality: I have seen around the forum that the propagators above may provide mean/osculating elements and I am not really an expert about what that means.

Since, for now, I am interested in the short term evolution - that is, one orbital period - of the relative position of the two spacecrafts in LEO, which propagator would provide the most accurate results in terms of relative along-track, across-track and radial distance? It would be nice to validate my results in some way.

Thanks in advance!

2

You should use numerical propagator, perhaps increasing gravity field to 60x60 as the altitude is very low and Earth gravity field irregularities are most sensitive when you are close to Earth. You could increase the max step to much larger than 2 seconds, you could even go as far as 120 or 600 seconds if you want, as long as you let the tolerance at 1.0e-3. Beware when you compute the tolerance to be sure to use the same OrbitType that will be used by the propagator.

For short term evolution, you need osculating elements, hence numerical propagator (or DSST with short periods, but it is more difficult to set up).

The fact changing the eccentricity changes the separation is normal, the satellites will move with respect to each other. You would also notice differences if you change slightly the places, either by changing a smidge inclination or right ascension of ascending node. In fact, using eccentricity and inclination separation is one way formation flying is performed: the satellites enter what is called a elliptical relative motion, with the center of the relative ellipse being driven by the true anomaly difference, the along track and radial dimension of the ellipse being driven by eccentricity difference and the across track harmonic being driven by planes differences. If for example you use zero true anomaly difference but non-zero eccentricity/plane separation, you will end up with one satellite revolving around the other one.

1 Like
3

Dear Luc,
thank you for your reply. I have implemented your comments, I wanted to know if from code side everything seems fine or if I am incorrectly propagating the orbits / obtaining the relative positions.

Note that I am using the same “propagationType” for the propagator definition and the tolerances, but does it matter if I choose keplerian or cartesian?
Furthermore, by choosing (60,60) for the gravity model, does this mean I am including harmonics parameters up to J60?

# Initial date
year = 2026
month = 1
day = 1
hour = 12
minute = 0
seconds = 00.000

duration = 5650. # Propagation duration [s]
step_time = 0.4  # Propagation step-time [s]

#Orbital parameters (keplerian)
a = 6890. * 1000    # Orbit semi-major axis [m]
e = 0.001       # Orbit eccentricity [-]
i = radians(97.4)       # Orbit inclination [rad]
omega = radians(90.0)     # Orbit perigee argument [rad]
raan = radians(50.0)      # Orbit RAAN [rad]
true_anomaly = radians(90.0)    # Spacecraft true anomaly [rad]

at_separation = 200      # Spacecraft along-track separation [m]

### here funciton to compute delta-true anomaly based on along track separation - skipped for easier reading ###

utc = TimeScalesFactory.getUTC()
epochDate = AbsoluteDate(year, month, day, hour, minute, seconds, utc)             # Propagation start-time [s]

# Initial orbit definition, spacecraft behind in the formation
initialOrbit = KeplerianOrbit(a, e, i, omega, raan, true_anomaly,
                              PositionAngle.TRUE,
                              inertialFrame, epochDate, Constants.EIGEN5C_EARTH_MU)

# Initial orbit definition, first spacecraft of the formation
initialOrbit2 = KeplerianOrbit(a, e, i, omega, raan, true_anomaly+dtheta,
                              PositionAngle.TRUE,
                              inertialFrame, epochDate, Constants.EIGEN5C_EARTH_MU)

satellite_mass = 20.

# Numerical propagator definition
# Step size
minStep = 0.001
maxStep = 10.
initStep = 0.4
positionTolerance = 1E-3
propagationType = OrbitType.CARTESIAN 
tolerances = NumericalPropagator.tolerances(positionTolerance, initialOrbit, propagationType)

# Spacecraft behind - numerical propagator
integrator = DormandPrince853Integrator(minStep, maxStep, orekit.JArray('double').cast_(tolerances[0]), orekit.JArray('double').cast_(tolerances[1]))
integrator.setInitialStepSize((initStep))

numericalPropagator = NumericalPropagator(integrator)
numericalPropagator.setOrbitType(propagationType)

itrf    = FramesFactory.getITRF(IERSConventions.IERS_2010, True) # International Terrestrial Reference Frame, earth fixed

earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
                         Constants.WGS84_EARTH_FLATTENING,
                         itrf)
gravityProvider = GravityFieldFactory.getNormalizedProvider(60, 60)
numericalPropagator.addForceModel(HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravityProvider))
numericalPropagator.setInitialState(SpacecraftState(initialOrbit, satellite_mass))

# First spacecraft in the formation - numerical propagator
tolerances2 = NumericalPropagator.tolerances(positionTolerance, initialOrbit2, propagationType)

integrator2 = DormandPrince853Integrator(minStep, maxStep, orekit.JArray('double').cast_(tolerances2[0]), orekit.JArray('double').cast_(tolerances2[1]))
integrator2.setInitialStepSize((initStep))

numericalPropagator2 = NumericalPropagator(integrator2)
numericalPropagator2.setOrbitType(propagationType)

numericalPropagator2.addForceModel(HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravityProvider))
numericalPropagator2.setInitialState(SpacecraftState(initialOrbit2, satellite_mass))

# Drag
crossSection_drag = 0.1
dragCoeff = 2.2

strenghtLevel = MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE
supportedNames = MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES
inputParameters_Drag = MarshallSolarActivityFutureEstimation(supportedNames, strenghtLevel)
atmosphere = NRLMSISE00(inputParameters_Drag, sun, earth)
spacecraftDrag = IsotropicDrag(crossSection_drag, dragCoeff)
drag = DragForce(atmosphere, spacecraftDrag)
numericalPropagator.addForceModel(drag)
numericalPropagator2.addForceModel(drag)

# SRP
crossSection_SRP = 0.1
reflectivityCoeff = 1.3

spacecraftSRP = IsotropicRadiationSingleCoefficient(crossSection_SRP, reflectivityCoeff)
solarRadiationPressure = SolarRadiationPressure(sun, R_e, spacecraftSRP)
numericalPropagator.addForceModel(solarRadiationPressure)
numericalPropagator2.addForceModel(solarRadiationPressure)

#Moon and Sun perturbations
moon = CelestialBodyFactory.getMoon()
sun = CelestialBodyFactory.getSun()
moon_3dbodyattraction = ThirdBodyAttraction(moon)
numericalPropagator.addForceModel(moon_3dbodyattraction)
numericalPropagator2.addForceModel(moon_3dbodyattraction)

sun_3dbodyattraction = ThirdBodyAttraction(sun)
numericalPropagator.addForceModel(sun_3dbodyattraction)
numericalPropagator2.addForceModel(sun_3dbodyattraction)

# Start propagations
initialDate = epochDate
t = [initialDate.shiftedBy(float(dt)) for dt in np.arange(0, duration, step_time)]

# Propagated orbits
spacecraft1 = [numericalPropagator.propagate(tt) for tt in t]
spacecraft2 = [numericalPropagator2.propagate(tt) for tt in t]

# Behind spacecraft coordinates
pv = [x.getPVCoordinates() for x in spacecraft1]
p = [tpv.getPosition() for tpv in pv]
v1 = [tpv.getVelocity() for tpv in pv]

# LVLH of behind spacecraft
lvlh = [LofOffset(sc1.getFrame(), LOFType.LVLH_CCSDS) for sc1 in spacecraft1]

# Conversion function
converted_sc1 = [SpacecraftState(placeholder1.getOrbit(),
                                      placeholder2.getAttitude(placeholder1.getOrbit(), placeholder1.getDate(), placeholder1.getFrame())) for placeholder1, placeholder2in zip(spacecraft1, lvlh)]

# Positions of spacecraft 2 (first) with respect to LVLH of spacecraft 1 (behind)
pv2t = [a.toTransform().transformPVCoordinates(y.getPVCoordinates(x.getFrame())) for a,y,x in zip(converted_sc1, spacecraft2, spacecraft1)]
4

This is the part where consistency is important.

No, you can use any representation you want, at numerical noise level it will give comparable results.

Yes, and you also includes the tesseral terms (the J_n terms are the zonal terms only, i.e. terms to degree n and order 0, whereas tesseral terms have order > 0).

I would recommend to rebuild force models for each propagator instead of reusing them in both propagators. It is cleaner as some force models may keep some internal state, so they really should not be shared. Sharing force models may work, but it’s better not to do so.

5

Thanks for your reply @luc !

I have one last question regarding the numerical propagator definition, since I believe it is really important that both orbits are propagated at the same instants in time, following the time-step separation which is 0.4 s.

By looking at the code above, is this the case? Because I am worried that selecting a range for minStep and maxStep may cause the propagator of spacecraft1 to derive its PV in different instants in time with respect to the propagator2.
This would mean that I would be extracting the LVLH orbital frame from spacecraft1 and comparing it to the PV of spacecraft2 in a different instant in time.

I have printed out the time-stamped values of each propagation step by using the following commands and they are exactly equal, but still I wanted to have confirmation (and a better understanding) on the matter:

dates_1 = [xx.getDate() for xx in spacecraft1]
dates_2 = [xx.getDate() for xx in spacecraft2]
6

Yes and no.
Yes because if you do a full long range propagation on both satellites, they can loose synchronization on intermediate steps because each integrator bases its step control decision on an estimate of the current error, which is based on the current state. Since the states will be different for the two satellites, error estimation will be different too.
No because if you don’t care about the intermediate steps and just use the result of the propagate method at final time tt, then you have already fixed this time by yourself independently of underlying integrator choices.

This leads me to one advice: don’t do the propagation loop as you did, by setting a series of time and calling the propagator in a loop, in particular if your various times tt are regularly spaced (resulting from t = [initialDate.shiftedBy(float(dt)) for dt in np.arange(0, duration, step_time)]). This is really not efficient, in particular with a 0.4s long propagation. You should rather perform a single propagation for the full duration (5650s in your case), and register a step handler to the propagator that will be called every step_time. With these settings, you let the integrator choose step sized that may be far larger than 0.4s (you have set a max step to 10s, but I would even suggest to let it go to 120s or more, as the real size will really be driven by the positionTolerance). The sampling time at which your handler will be called (here step_time) is completely independent of the step the integrator uses for its computation, the integrator knows how to interpolate within the step using its own internal model, which will be consistent with the dynamics since it is really the model the integrator itself uses.

As you want to have several satellites propagated at once, you will need to use PropagatorParallelizer, together with a MultiSatFixedStepHandler. Beware however that the fixed step handler for propagator parallelizer is a new feature introduced in Orekit 12.0, that should be released probably tomorrow, and the Python wrapper will be published soon after (see [VOTE] Releasing Orekit 12.0 from release candidate 2 - #13 by petrus.hyvonen).

7

Dear @luc , thank you for your reply. I’ll try with the solution you proposed with the parallelizer when the new release becomes available.

Since I am interested in short-term variations - that is, analyzing the distances along one single orbit only - I don’t know whether increasing the step time could be applied to my analysis.

Do you believe that the current approach shown in the previous codes may differ largely with respect to a parallelized propagator? Is there a way to validate my results, to your knowledge?