Numerical propagation, point insider ellipsoid

Orekit Python Wrapper usage
1

Hi everyone,

I’m working on a Python script to validate the accuracy of a propagator I developed. To do this, I decided to use Orekit as a reference to check whether my propagator meets the desired precision.

However, I’m encountering issues with the propagator setup. I’ve tried configuring the models and parameters to achieve the highest possible accuracy (I also plan to compare my results with the JPL Horizons dataset). Despite this, there’s a problem during propagation: the satellite falls to Earth’s surface, and I’m not sure why.

It might be related to the initial state, but I’m extracting it from the Hubble dataset in JPL Horizons at start_date. I’ve also tried disabling drag force and modifying drag and solar pressure surfaces, but the issue persists.

Note: The code is not well-structured—it’s just a quick test to compare Orekit’s results with other propagators.

Below, you’ll find the Python script and output logs.

Any questions or suggestions would be greatly appreciated!

import orekit
from orekit.pyhelpers import setup_orekit_curdir, download_orekit_data_curdir
import os
import logging
from java.io import File
from org.orekit.utils import IERSConventions, Constants
from org.orekit.propagation.numerical import NumericalPropagator
from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel
from org.orekit.forces.gravity.potential import GravityFieldFactory
from org.orekit.forces.drag import DragForce, IsotropicDrag
from org.orekit.forces.radiation import SolarRadiationPressure, IsotropicRadiationSingleCoefficient
from org.orekit.models.earth.atmosphere import NRLMSISE00
from org.orekit.models.earth.atmosphere.data import CssiSpaceWeatherData
from org.orekit.bodies import CelestialBodyFactory, OneAxisEllipsoid
from org.orekit.frames import FramesFactory
from org.orekit.orbits import CartesianOrbit, OrbitType
from org.orekit.time import AbsoluteDate, TimeScalesFactory
from org.orekit.utils import PVCoordinates
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator
from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.orekit.propagation import SpacecraftState
from org.orekit.forces.gravity import ThirdBodyAttraction
import math
import numpy as np

# Initialize Orekit
orekit.initVM()
setup_orekit_curdir()

# Download Orekit data if not already available
if not os.path.isdir('orekit-data'):
    download_orekit_data_curdir()

# Configure logging
logging.basicConfig(
    level=logging.INFO,
    format='%(asctime)s [%(levelname)s] %(message)s',
)
logger = logging.getLogger(__name__)

initial_state = [-4202834.74, -5275624.06, -1327977.87,
                 5740.649980970608, -3738.503978889843, -3330.233806655470]

mass = 11063.1  # kg
drag_coeff = 2.0
drag_area = 25.0  # m²
srp_coeff = 1.5
srp_area = 25.0  # m²
start_date = AbsoluteDate(2025, 3, 4, 23, 58, 50.815, TimeScalesFactory.getUTC())
# end_date = start_date.shiftedBy(100.0)  # 100 segundos para prueba inicial
end_date = AbsoluteDate(2025, 3, 9, 23, 58, 50.815, TimeScalesFactory.getUTC())  # 5 días
step_size = 10.0  # seconds
file_path = 'propagation_data.txt'
include_drag = False

# 1. Set up Earth parameters
earth = CelestialBodyFactory.getEarth()
iers_convention = IERSConventions.IERS_2010
earthFrame = FramesFactory.getITRF(iers_convention, True)
earthRadius = Constants.WGS84_EARTH_EQUATORIAL_RADIUS
earthFlattening = Constants.WGS84_EARTH_FLATTENING
earthShape = OneAxisEllipsoid(earthRadius, earthFlattening, earthFrame)

# 2. Set up Earth gravity model
gravity_provider = GravityFieldFactory.getNormalizedProvider(20, 20)
gravity_model = HolmesFeatherstoneAttractionModel(earthFrame, gravity_provider)

# 3. Create initial orbit in ICRF
pvCoordinates = PVCoordinates(
    Vector3D(initial_state[0], initial_state[1], initial_state[2]),
    Vector3D(initial_state[3], initial_state[4], initial_state[5])
)
initial_orbit = CartesianOrbit(pvCoordinates, FramesFactory.getICRF(), start_date, Constants.WGS84_EARTH_MU)

# 4. Set up numerical propagator
min_step = 1e-6
max_step = 60.0
position_tolerance = 10.0
integrator = DormandPrince853Integrator(min_step, max_step, 1e-12, 1e-12)
propagator = NumericalPropagator(integrator)
propagator.setOrbitType(OrbitType.CARTESIAN)
propagator.setInitialState(SpacecraftState(initial_orbit, mass))
propagator.clearEventsDetectors()  # Desactivar detectores automáticos

# 5. Add force models
propagator.addForceModel(gravity_model)
cssiSpaceWeatherData = CssiSpaceWeatherData("SpaceWeather-All-v1.2.txt")
atmosphere = NRLMSISE00(cssiSpaceWeatherData, CelestialBodyFactory.getSun(), earthShape)
dragSC = IsotropicDrag(drag_area, drag_coeff)
drag_force = DragForce(atmosphere, dragSC)
propagator.addForceModel(drag_force)
sun = CelestialBodyFactory.getSun()
earth = OneAxisEllipsoid(earthRadius, Constants.WGS84_EARTH_FLATTENING, FramesFactory.getICRF())
isotropicRadiationSingleCoeff = IsotropicRadiationSingleCoefficient(srp_area, srp_coeff)
srpForce = SolarRadiationPressure(sun, earth, isotropicRadiationSingleCoeff)
propagator.addForceModel(srpForce)
propagator.addForceModel(ThirdBodyAttraction(CelestialBodyFactory.getSun()))
propagator.addForceModel(ThirdBodyAttraction(CelestialBodyFactory.getMoon()))

# 6. Propagation and data output
logger.info(f"Starting propagation at {start_date}")
try:
    with open(file_path, 'w') as file:
        # Escribir el estado inicial
        initial_utc = start_date.toString(TimeScalesFactory.getUTC()).replace('T', ' ')
        initial_line = f"{initial_utc},{initial_state[0]},{initial_state[1]},{initial_state[2]},{initial_state[3]},{initial_state[4]},{initial_state[5]}\n"
        file.write(initial_line)

        logger.info(f"Initial altitude: {pvCoordinates.getPosition().getNorm() - earthRadius} m")

        # Propagación en intervalos de 10 segundos
        sample_date = start_date
        while sample_date.compareTo(end_date) <= 0:
            state = propagator.propagate(sample_date)
            pv = state.getPVCoordinates(FramesFactory.getICRF())
            pos = pv.getPosition()
            vel = pv.getVelocity()
            altitude = pos.getNorm() - earthRadius

            logger.info(f"Date: {state.getDate()}, Altitude: {altitude} m")

            if altitude < 0:
                logger.error(f"Altitude below surface at {state.getDate()}: {altitude} m")
                break

            # Guardar en el archivo
            utc_date = state.getDate().toString(TimeScalesFactory.getUTC()).replace('T', ' ')
            line = f"{utc_date},{pos.getX()},{pos.getY()},{pos.getZ()},{vel.getX()},{vel.getY()},{vel.getZ()}\n"
            file.write(line)

            # Avanzar 10 segundos
            sample_date = sample_date.shiftedBy(10.0)
except Exception as e:
    logger.error(f"An error occurred: {e}")
finally:
    logger.info(f"Propagation finished at {end_date}")
    logger.info(f"Propagation data saved to {file_path}")

2025-03-12 15:38:56,699 [INFO] Starting propagation at 2025-03-04T23:58:50.815Z
2025-03-12 15:38:56,699 [INFO] Initial altitude: 496421.48017976247 m
2025-03-12 15:38:56,699 [INFO] Date: 2025-03-04T23:58:50.815Z, Altitude: 496421.48017976247 m
2025-03-12 15:38:57,972 [INFO] Date: 2025-03-04T23:59:00.815Z, Altitude: 496514.56888198294 m
2025-03-12 15:38:57,995 [INFO] Date: 2025-03-04T23:59:10.815Z, Altitude: 496739.52475478593 m
2025-03-12 15:38:58,029 [INFO] Date: 2025-03-04T23:59:20.815Z, Altitude: 497095.1619379893 m
2025-03-12 15:38:58,048 [INFO] Date: 2025-03-04T23:59:30.815Z, Altitude: 497580.1849763403 m
2025-03-12 15:38:58,081 [INFO] Date: 2025-03-04T23:59:40.815Z, Altitude: 498193.18983183615 m
2025-03-12 15:38:58,097 [INFO] Date: 2025-03-04T23:59:50.815Z, Altitude: 498932.6649817908 m
2025-03-12 15:38:58,104 [INFO] Date: 2025-03-05T00:00:00.815Z, Altitude: 499796.99260077346 m
2025-03-12 15:38:58,125 [INFO] Date: 2025-03-05T00:00:10.815Z, Altitude: 500784.44982437603 m
2025-03-12 15:38:58,143 [INFO] Date: 2025-03-05T00:00:20.815Z, Altitude: 501893.2100926293 m
2025-03-12 15:38:58,165 [INFO] Date: 2025-03-05T00:00:30.815Z, Altitude: 503121.3445707932 m
...
2025-03-12 15:38:59,201 [INFO] Date: 2025-03-05T00:24:20.815Z, Altitude: 169329.2349564666 m
2025-03-12 15:38:59,206 [INFO] Date: 2025-03-05T00:24:30.815Z, Altitude: 152559.3680682471 m
2025-03-12 15:38:59,210 [INFO] Date: 2025-03-05T00:24:40.815Z, Altitude: 135526.8768665064 m
2025-03-12 15:38:59,216 [INFO] Date: 2025-03-05T00:24:50.815Z, Altitude: 118233.5797212068 m
2025-03-12 15:38:59,218 [INFO] Date: 2025-03-05T00:25:00.815Z, Altitude: 100681.42329558823 m
2025-03-12 15:38:59,218 [INFO] Date: 2025-03-05T00:25:10.815Z, Altitude: 82872.48486507777 m
2025-03-12 15:38:59,228 [INFO] Date: 2025-03-05T00:25:20.815Z, Altitude: 64808.97468862124 m
2025-03-12 15:38:59,232 [INFO] Date: 2025-03-05T00:25:30.815Z, Altitude: 46493.238432354294 m
2025-03-12 15:38:59,235 [INFO] Date: 2025-03-05T00:25:40.815Z, Altitude: 27927.759645403363 m
2025-03-12 15:38:59,239 [INFO] Date: 2025-03-05T00:25:50.815Z, Altitude: 9115.1622874178 m
2025-03-12 15:38:59,250 [ERROR] An error occurred: <super: <class 'JavaError'>, <JavaError object>>
    Java stacktrace:
org.orekit.errors.OrekitException: el punto se encuentra en el interior del elipsoide
        at org.orekit.bodies.Ellipsoid.pointOnLimb(Ellipsoid.java:407)
        at org.orekit.utils.OccultationEngine.angles(OccultationEngine.java:83)
        at org.orekit.propagation.events.EclipseDetector.g(EclipseDetector.java:214)
        at org.orekit.propagation.integration.AbstractIntegratedPropagator$AdaptedEventDetector.g(AbstractIntegratedPropagator.java:985)
        at org.hipparchus.ode.events.DetectorBasedEventState.g(DetectorBasedEventState.java:154)
        at org.hipparchus.ode.events.DetectorBasedEventState.evaluateStep(DetectorBasedEventState.java:226)
        at org.hipparchus.ode.AbstractIntegrator.lambda$acceptStep$6(AbstractIntegrator.java:327)
        at java.util.ArrayList$ArrayListSpliterator.forEachRemaining(ArrayList.java:1384)
        at java.util.stream.Streams$ConcatSpliterator.forEachRemaining(Streams.java:742)
        at java.util.stream.ReferencePipeline$Head.forEach(ReferencePipeline.java:647)
        at org.hipparchus.ode.AbstractIntegrator.acceptStep(AbstractIntegrator.java:327)
        at org.hipparchus.ode.nonstiff.EmbeddedRungeKuttaIntegrator.integrate(EmbeddedRungeKuttaIntegrator.java:267)
        at org.orekit.propagation.integration.AbstractIntegratedPropagator.integrateDynamics(AbstractIntegratedPropagator.java:509)
        at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:440)
        at org.orekit.propagation.integration.AbstractIntegratedPropagator.propagate(AbstractIntegratedPropagator.java:400)

2025-03-12 15:38:59,250 [INFO] Propagation finished at 2025-03-09T23:58:50.815Z
2025-03-12 15:38:59,251 [INFO] Propagation data saved to propagation_data.txt
2

I do not know it this is what is causing the error, but 10m for positionTolerance seems very high: I have seen 1E-3 to be typically recommended.

3

I changed it, but the output remains the same. However, I agree that 1e-3 is better.

4

Are you sure the coordinates are defined in ICRF?
Since ICRF is centered on the solar system’s barycenter, your satellite seems far from the Earth.

2 Likes
5

Perfect, thank you so much! It’s working great now.

2 Likes