Error during DragForce.acceleration | Clarity on input parameter

Orekit usage
1

I need to extract the drag acceleration during propagation. For that I am using the acceleration method associated with the drag ForceModel. But in doing so, I am running into error saying: No matching overloads found for org.orekit.forces.drag.DragForce.acceleration(org.orekit.propagation.SpacecraftState,JDouble)

I understand that the inputs need to be of the form: acceleration(Spacecraft s, double[] parameters) and not JDouble. Can someone clarify what quantities this array of doubles (“parameters”) should comprise?

I have set up my drag model as follows:

atmosphere = NRLMSISE00(nrlmsis_prams,sun,earth)
spacecraftDrag = IsotropicDrag(crossSectionArea, dragCoefficient)
dragForceModel = DragForce(atmosphere, spacecraftDrag)
propagator_num.addForceModel(dragForceModel)

During propagation i am doing the following that leads to the error:

drag_acc_vec = dragForceModel.acceleration(s, dragForceModel.getParametersDrivers().get(0).getValue()).toArray()[:]

During debugging I notice that self.drag_model.getParametersDrivers() gives
<UnmodifiableRandomAccessList: [global drag factor = 1.0, drag coefficient = 2.2]>

Out of these parameters I am unsure which one to use as input to the acceleration method.

Any advice will be much appreciated. Thanks.

2

Hello @Stella28,

Compute acceleration

You can do something much simpler with the following code (in Java for now but don’t worry) :

forceModel.acceleration(spacecraftState, forceModel.getParameters(spacecraftState.getDate()))

This will directly give you the array needed to compute the acceleration i.e. the cross section and drag coefficient parameters in your case.

Store acceleration

You did not mention it so i will mention it just in case but you can use an OrekitFixedStepHandler to save the drag acceleration during propagation. For that you will need to implement your own class that will take your drag force as a parameter and then call its acceleration methods in the handleStep method of your step handler.

Code example

For this example, i’m using the classic Orekit Labs / Orekit Python Wrapper · GitLab so it will be simpler if you are using orekit jpype.

Don’t forget to set your own path to orekit data !

import orekit
from orekit import JArray_double
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator
from org.orekit.bodies import CelestialBodyFactory
from org.orekit.forces import ForceModel
from org.orekit.forces.drag import DragForce
from org.orekit.forces.drag import IsotropicDrag
from org.orekit.frames import FramesFactory
from org.orekit.models.earth import ReferenceEllipsoid
from org.orekit.models.earth.atmosphere import NRLMSISE00
from org.orekit.models.earth.atmosphere.data import CssiSpaceWeatherData
from org.orekit.orbits import EquinoctialOrbit
from org.orekit.orbits import PositionAngleType
from org.orekit.propagation import Propagator
from org.orekit.propagation import SpacecraftState
from org.orekit.propagation.numerical import NumericalPropagator
from org.orekit.propagation.sampling import PythonOrekitFixedStepHandler
from org.orekit.time import AbsoluteDate
from org.orekit.time import TimeScalesFactory
from org.orekit.utils import Constants
from org.orekit.utils import IERSConventions
from org.orekit.utils import ParameterDriversProvider

# Initialize JVM
vm = orekit.initVM()

from orekit.pyhelpers import setup_orekit_curdir

# Set orekit data
setup_orekit_curdir('YOUR PATH TO OREKIT DATA')

### CONSTANTS
UTC = TimeScalesFactory.getUTC()
INERTIAL = FramesFactory.getEME2000()
ITRF = FramesFactory.getITRF(IERSConventions.IERS_2010, False)
EARTH_BODYSHAPE = ReferenceEllipsoid.getIers2010(ITRF)
EARTH_RADIUS = Constants.WGS84_EARTH_EQUATORIAL_RADIUS
SUN = CelestialBodyFactory.getSun()
INITIAL_DATE = AbsoluteDate(2024, 8, 30, 15, 0, 00.000, UTC)


### CLASSES
class AccelerationHandler(PythonOrekitFixedStepHandler):

    def __init__(self, force: ForceModel):
        super().__init__()
        self.force = force
        self.accelerations = []

    def finish(self, spacecraftState):
        pass

    def init(self, spacecraftState, absoluteDate, double):
        self.accelerations = []

    def handleStep(self, spacecraftState):
        casted_force = ParameterDriversProvider.cast_(self.force)
        acceleration = self.force.acceleration(spacecraftState, casted_force.getParameters(spacecraftState.getDate()))
        self.accelerations.append(list(acceleration.toArray()))


### METHODS
def get_numerical_propagator(initial_orbit):
    min_step = 0.001
    maxstep = 1000.0
    init_step = 60.0
    position_tolerance = 1.0

    # Compute tolerances
    tolerances = NumericalPropagator.tolerances(position_tolerance,
                                                initial_orbit,
                                                initial_orbit.getType())
    # Get integrator
    integrator = DormandPrince853Integrator(min_step, maxstep,
                                            JArray_double.cast_(tolerances[0]),
                                            # Double array of doubles needs to be casted in Python
                                            JArray_double.cast_(tolerances[1]))
    integrator.setInitialStepSize(init_step)

    # Create numerical propagator
    return NumericalPropagator(integrator)


def get_orbit():
    # Define equatorial orbit
    a = EARTH_RADIUS + 450e3
    return EquinoctialOrbit(a, 0., 0., 0., 0., 0., PositionAngleType.TRUE, INERTIAL,
                            INITIAL_DATE, Constants.WGS84_EARTH_MU)


### MAIN
if __name__ == '__main__':
    # Get initial orbit
    initial_orbit = get_orbit()

    # Get numerical propagator
    propagator_num = get_numerical_propagator(initial_orbit)

    # Create initial state
    satellite_mass = 1000.0  # The models need a spacecraft mass, unit kg.
    initial_state = SpacecraftState(initial_orbit, satellite_mass)

    # Set initial state
    propagator_num.setInitialState(initial_state)

    # Add drag
    cross_section = 1.
    drag_coeff = 2.2
    drag_sensitive = IsotropicDrag(cross_section, drag_coeff)
    atmosphere = NRLMSISE00(CssiSpaceWeatherData(CssiSpaceWeatherData.DEFAULT_SUPPORTED_NAMES), SUN, EARTH_BODYSHAPE)
    drag_force = DragForce(atmosphere, drag_sensitive)

    propagator_num.addForceModel(drag_force)

    # Set step handler
    force_handler = AccelerationHandler(drag_force)

    Propagator.cast_(propagator_num).setStepHandler(60., force_handler)

    # Propagate for one day
    propagator_num.propagate(initial_orbit.getDate().shiftedBy(3600.))

    # Display forces
    print(force_handler.accelerations)

Hope this helps,

Cheers,
Vincent

3

@Vincent Thank you very much for the detailed answer. That was immensely helpful.

1 Like