Dear all,

I am currently trying to measure the effects of higher order perturbations and drag on the nominal trajectory of a satellite over 3 days. The method I am using is:

- initialize the spacecraft state with its osculating keplerian orbit parameters
- propagate the spacecraft state considering only gravitational perturbations up to a certain degree
- propagate the spacecraft state with a different propagator, which considers higher order gravitational terms and drag
- evaluate the position difference between the two satellites, centered in the LVLH frame of the nominal spacecraft.

In this way, I should be able to evaluate the relative position between the nominal trajectory and the perturbed one, and the choice of coordinates allows me to evaluate whether the real trajectory is within the boundaries of an orbital tube, centered in the nominal trajectory, with components on the across-track and radial direction.

The method I am using to compute at each time the LVLH frame of the nominal trajectory is taken from this thread.

What I expect to see is a slight difference on the along-track direction, an oscillation for what regards the across-track direction and a constant increase for the radial direction, due to drag.

What I observe is that the component in the along-track direction rapidly diverges and reaches thousands of kilometers in a few days! The across-track component oscillates around 0 and the radial increases, which is something that I would expect instead (but I am not certain about the magnitude)

Here is the code that I am using, did I implement something wrongly?

```
year = 2026
month = 1
day = 1
hour = 6
minute = 0
seconds = 03.146
date_propagation = '2026-01-01T06:00:03.416'
utc = TimeScalesFactory.getUTC()
epochDate = AbsoluteDate(year, month, day, hour, minute, seconds, utc)
epoch = Time(date_propagation)
duration = 3*24*60*60. # Propagation duration [s]
J2 = -Constants.EIGEN5C_EARTH_C20
h = 500e3 # Orbit altitude [m]
R_e = Constants.WGS84_EARTH_EQUATORIAL_RADIUS
a = Constants.WGS84_EARTH_EQUATORIAL_RADIUS + h
e = 0.0
mu_e = Constants.EIGEN5C_EARTH_MU
i = radians(60)
omega = radians(90.0)
raan = radians(0.0)
true_anomaly = radians(0.0)
inertialFrame = FramesFactory.getEME2000()
# Initial orbit definition
initialOrbit = KeplerianOrbit(a, e, i, omega, raan, true_anomaly,
PositionAngle.TRUE,
inertialFrame, epochDate, Constants.EIGEN5C_EARTH_MU)
period = initialOrbit.getKeplerianPeriod()
step_time = period/300 # Propagation step-time [s]
satellite_mass = 70. # Spacecraft mass [kg]
# Numerical propagator definition - nominal
minStep = 0.001
maxStep = step_time
initStep = step_time
positionTolerance = 1E-2
propagationType = OrbitType.CARTESIAN
tolerances = NumericalPropagator.tolerances(positionTolerance, initialOrbit, propagationType)
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)
gravityProvider1 = GravityFieldFactory.getNormalizedProvider(20, 20)
numericalPropagator.addForceModel(HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravityProvider1))
numericalPropagator.setInitialState(SpacecraftState(initialOrbit, satellite_mass))
# Numerical propagator definition - perturbed
tolerances2 = NumericalPropagator.tolerances(positionTolerance, initialOrbit, 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)
gravityProvider2 = GravityFieldFactory.getNormalizedProvider(60, 60)
numericalPropagator2.addForceModel(HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravityProvider2))
numericalPropagator2.setInitialState(SpacecraftState(initialOrbit, satellite_mass))
crossSection_drag = float ( 2 ) # Cross section area used for drag
dragCoeff = 2.2 # Spacecraft drag coefficient
# Initialize drag perturbation
sun = CelestialBodyFactory.getSun() # Create sun in simulation
strenghtLevel = MarshallSolarActivityFutureEstimation.StrengthLevel.STRONG
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)
numericalPropagator2.addForceModel(drag)
initialDate = epochDate
t = [initialDate.shiftedBy(float(dt)) for dt in np.arange(0, duration+step_time, step_time)]
spacecraft1 = [numericalPropagator.propagate(tt) for tt in t] # Nominal
spacecraft2 = [numericalPropagator2.propagate(tt) for tt in t] # Perturbed
pv = [x.getPVCoordinates() for x in spacecraft1]
lvlh = [LofOffset(sc.getFrame(), LOFType.LVLH_CCSDS) for sc in spacecraft1]
converted_sc1 = [SpacecraftState(sc1.getOrbit(), lvlh_list.getAttitude(sc1.getOrbit(), sc1.getDate(), sc1.getFrame())) for sc1, lvlh_list in zip(spacecraft1, lvlh)]
p1_along = [x.getPosition().getX() for x in pv]
p1_normal = [y.getPosition().getY() for y in pv]
p1_radial = [z.getPosition().getZ() for z in pv]
pv2 = [y.getPVCoordinates() for y in spacecraft2]
pv2t = [a.toTransform().transformPVCoordinates(y.getPVCoordinates(x.getFrame())) for a,y,x in zip(converted_sc1, spacecraft2, spacecraft1)]
pv2t_along = [x.getPosition().getX() for x in pv2t]
pv2t_normal = [y.getPosition().getY() for y in pv2t]
pv2t_radial = [z.getPosition().getZ() for z in pv2t]
plt.plot(np.arange(0, duration+step_time, step_time)/24/60/60, pv2t_along)
plt.title('Along track difference [m]')
plt.show()
plt.plot(np.arange(0, duration+step_time, step_time)/24/60/60, pv2t_normal)
plt.title('Cross track difference [m]')
plt.show()
plt.plot(np.arange(0, duration+step_time, step_time)/24/60/60, pv2t_radial)
plt.title('Radial difference [m]')
plt.show()
```

Another way I could perform my activity is to evaluate the variation in the keplerian parameters and derive the position difference from these, but I would like to maintain this approach, if possible.