The drag force is not effective, why?

Dear all,

I am so new in working with Orekit and so happy to be a member of Orekit community. As a starter, I tried to simulate an object propagation considering only the Earth gravitational force and the atmospheric drag force. I chose the orbital element in a way that it is in LEO region, so the drag should be effective. However, when I comment and uncomment the line regarding the drag force, the output is the same. I’d be grateful if anyone can give me any clue about my mistakes in this regard.

I appreciate your time in advance.
You can take a look into my code as follows:

import orekit
vm = orekit.initVM()
# Now we have populated the namespace with the orekit classes

from orekit.pyhelpers import setup_orekit_curdir

from java.util import Arrays
from orekit import JArray_double

from org.orekit.bodies import  OneAxisEllipsoid
from org.orekit.frames import  FramesFactory
from import DataProvidersManager, ZipJarCrawler
from org.orekit.time import TimeScalesFactory, AbsoluteDate
from org.orekit.orbits import KeplerianOrbit
from org.orekit.utils import Constants
from org.orekit.propagation.analytical import EcksteinHechlerPropagator
from org.orekit.propagation.analytical.tle import TLEPropagator
from org.orekit.propagation.conversion import FiniteDifferencePropagatorConverter
from org.orekit.propagation.conversion import TLEPropagatorBuilder
from datetime import datetime
from org.orekit.propagation import SpacecraftState
from org.orekit.orbits import OrbitType, PositionAngle
from org.orekit.propagation.numerical import NumericalPropagator
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator
from org.orekit.forces.gravity.potential import GravityFieldFactory
from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel
from org.orekit.utils import IERSConventions

#Cartesian Elements
from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.orekit.utils import PVCoordinates

#for drag force
from import CssiSpaceWeatherData, MarshallSolarActivityFutureEstimation
from import NRLMSISE00, DTM2000
from org.orekit.forces.drag import IsotropicDrag
from org.orekit.forces.drag import DragForce

#For third bodies
from org.orekit.bodies import CelestialBodyFactory

import matplotlib.dates as mdates
from math import radians
import numpy as np
import matplotlib.pyplot as plt

# Some Constants
mu = Constants.WGS84_EARTH_MU
utc = TimeScalesFactory.getUTC()

ra = 10290 * 1000         # km Apogee
rp = 7284 * 1000         # km Perigee
i = radians(13.2)      # inclination
omega = radians(20.07)   # perigee argument
raan = radians(255.03)  # right ascension of ascending node
lv = radians(28.45)    # True anomaly
epochDate = AbsoluteDate(2012, 1, 26, 16, 0, 00.000, utc)

a = (rp + ra) / 2.0
e = (ra - rp)/(rp + ra)

#Inertial frame where the satellite is defined
inertialFrame = FramesFactory.getEME2000()

#Orbit construction as Keplerian
initialOrbit = KeplerianOrbit(a, e, i, omega, raan, lv,
                              inertialFrame, epochDate, mu)
satellite_mass = 10.0 # kg

initialState = SpacecraftState(initialOrbit, satellite_mass)

minStep = 0.001
maxstep = 1000.0
initStep = 60.0

positionTolerance = 10.0 
orbitType = OrbitType.KEPLERIAN
tol = NumericalPropagator.tolerances(positionTolerance, initialOrbit, orbitType)

integrator = DormandPrince853Integrator(minStep, maxstep, 
    JArray_double.cast_(tol[0]),  # Double array of doubles needs to be casted in Python

propagator_num = NumericalPropagator(integrator)

itrf    = FramesFactory.getITRF(IERSConventions.IERS_2010, True) # International Terrestrial Reference Frame, earth fixed
earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
gravityProvider = GravityFieldFactory.getNormalizedProvider(8, 8)

##Define Drag Force
sun = CelestialBodyFactory.getSun()
cswl = CssiSpaceWeatherData("SpaceWeather-All-v1.2.txt")
msafe = MarshallSolarActivityFutureEstimation(
#atmosphere = NRLMSISE00(msafe, sun, earth)
atmosphere = DTM2000(msafe, sun, earth)
cd = 1.1
cross_section = 9.6
isotropic_drag = IsotropicDrag(cross_section, cd)
drag_force = DragForce(atmosphere, isotropic_drag)

propagator_num.addForceModel(HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravityProvider))

## Create time vector
startDate = AbsoluteDate(2012, 1, 26, 11, 0, 00.000, utc)

#Overall duration in seconds for extrapolation
duration = 86400 * 10 
step_time = 60 * 3

# Time array in orekit AbsoluteDate format
t = [startDate.shiftedBy(float(dt)) \
        for dt in np.arange(0, duration, step_time)]

state = [propagator_num.propagate(tt) for tt in t]
pos = [propagator_num.propagate(tt).getPVCoordinates().getPosition() for tt in t]

#for i in pos:
#    print(i.x)

plt.plot([n.x for n in pos])

Hi @Sajji

Welcome to the Orekit forum!

The perigee altitude is 906 km and the apogee altitude is 3882 km.
The drag effect is close to negligeable for this orbit.

Just a test to see if there is no mistakes in the code, could you set:

ra = 6828* 1000         # km Apogee
rp = 6778 * 1000         # km Perigee

In that case Perigee altitude is 400 km and apogee altitude is 450 km. The purpose is to see if drag as an effect (it should have) to confirm there is no bug in your code.

If drag has no effect, could you check the epoch of the position you print. You can check it using the state parameter. Looking at the code, they should be at the same epoch.


Dear @bcazabonne

I appreciate your time and consideration. It worked perfectly. Just as a quick question, which model is working better for modeling drag force in LEO region according to your experiences?

Again I thank you very much.

I appreciate the NRLMSISE00 atmosphere model.
For input data, I always wonder which ones are the best. Personally, I prefer to use the CssiSpaceWeatherData

Best regards,