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
setup_orekit_curdir("resources")
from java.util import Arrays
from orekit import JArray_double
from org.orekit.bodies import OneAxisEllipsoid
from org.orekit.frames import FramesFactory
from org.orekit.data 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 org.orekit.models.earth.atmosphere.data import CssiSpaceWeatherData, MarshallSolarActivityFutureEstimation
from org.orekit.models.earth.atmosphere 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
ae = Constants.WGS84_EARTH_EQUATORIAL_RADIUS
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,
PositionAngle.TRUE,
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
JArray_double.cast_(tol[1]))
integrator.setInitialStepSize(initStep)
propagator_num = NumericalPropagator(integrator)
propagator_num.setOrbitType(orbitType)
propagator_num.setInitialState(initialState)
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(8, 8)
##Define Drag Force
sun = CelestialBodyFactory.getSun()
cswl = CssiSpaceWeatherData("SpaceWeather-All-v1.2.txt")
msafe = MarshallSolarActivityFutureEstimation(
MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES,
MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE
)
#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))
#propagator_num.addForceModel(drag_force)
## 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]
print(pos[0:5])
#for i in pos:
# print(i.x)
plt.plot([n.x for n in pos])
plt.show()