Negative eccentricity error when propagating numerically

Orekit Python Wrapper usage
1

Hello,

I’m building a numerical propagator class with the Python wrapper, and facing a very strange error:

org.orekit.errors.OrekitException: invalid parameter eccentricity: -0 not in range [0, ∞]

I have already tried some different techniques, but the problem persists. It propagates until a certain point, and then throws the above error. My code below:

import orekit
vm = orekit.initVM()
from orekit.pyhelpers import setup_orekit_curdir
setup_orekit_curdir()

from orekit.pyhelpers import absolutedate_to_datetime
from orekit.pyhelpers import datetime_to_absolutedate

from orekit import JArray_double

from org.orekit.orbits import KeplerianOrbit
from org.orekit.orbits import CartesianOrbit
from org.orekit.orbits import OrbitType, PositionAngle

from org.orekit.propagation import SpacecraftState
from org.orekit.propagation.numerical import NumericalPropagator
from org.orekit.propagation.analytical.tle import TLE, TLEPropagator
from org.orekit.propagation.analytical.tle import SGP4
from org.orekit.propagation.analytical import KeplerianPropagator
from org.orekit.propagation.analytical import EcksteinHechlerPropagator

from org.orekit.attitudes import NadirPointing

from org.hipparchus.ode.nonstiff import DormandPrince853Integrator
from org.hipparchus.ode.nonstiff import ClassicalRungeKuttaIntegrator
from org.hipparchus.geometry.euclidean.threed import Vector3D

from org.orekit.models.earth.atmosphere import NRLMSISE00
from org.orekit.models.earth.atmosphere.data import MarshallSolarActivityFutureEstimation

from org.orekit.forces.gravity.potential import GravityFieldFactory
from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel
from org.orekit.forces.gravity import ThirdBodyAttraction
from org.orekit.forces.gravity import Relativity

from org.orekit.forces.radiation import IsotropicRadiationSingleCoefficient
from org.orekit.forces.radiation import SolarRadiationPressure

from org.orekit.forces.drag import IsotropicDrag
from org.orekit.forces.drag import DragForce

from org.orekit.frames import FramesFactory, TopocentricFrame
from org.orekit.bodies import OneAxisEllipsoid, GeodeticPoint, CelestialBodyFactory
from org.orekit.utils import Constants
from org.orekit.utils import IERSConventions
from org.orekit.utils import PVCoordinates

from org.orekit.time import TimeScalesFactory, AbsoluteDate

import time
from datetime import datetime
from math import degrees, radians

utc = TimeScalesFactory.getUTC()
gcrf = FramesFactory.getGCRF()
itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, False)
eci  = gcrf  # Inertial frame
ecef = itrf  # Non-inertial frame
earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, Constants.WGS84_EARTH_FLATTENING, ecef)
sun  = CelestialBodyFactory.getSun()
moon = CelestialBodyFactory.getMoon()

class Spacecraft:
	def __init__(self, name, mass, cross_section, cd, cr, identification, norad_id):
		self.name           = name
		self.mass           = mass
		self.cross_section  = cross_section
		self.cd             = cd
		self.cr             = cr
		self.identification = identification
		self.norad_id       = norad_id

class KeplerOrbit:
    def __init__(self, apogee, perigee, inclination, aop, raan, lv):
	    self.apogee         = apogee*1000
	    self.perigee        = perigee*1000
	    self.inclination    = radians(inclination)
	    self.aop            = radians(aop)
	    self.raan           = radians(raan)
	    self.lv             = radians(lv)

    def sma(self):
	    return (self.perigee + self.apogee + 2*Constants.WGS84_EARTH_EQUATORIAL_RADIUS)/2.0

    def eccentricity(self):
	    return abs((self.apogee - self.perigee)/(self.apogee + self.perigee + 2*Constants.WGS84_EARTH_EQUATORIAL_RADIUS))

class KeplerPropagator:
	def __init__(self, spacecraft, orbit, epoch):
		self.spacecraft = spacecraft
		self.orbit   = orbit
		self.epoch = epoch

	def numericalPropagatorBuilder(self, duration, step, degree, order,
		anomaly = 'TRUE',
		thirdBody = True,
		solarPressure = True,
		atmosphericDrag = True,
		relativity = True):

	    startTime = AbsoluteDate(self.epoch.year,self.epoch.month,self.epoch.day,
		self.epoch.hour,self.epoch.minute,float(self.epoch.second),utc)
	    endTime = startTime.shiftedBy(duration)

	    # Orbit construction as Keplerian
	    if anomaly == 'TRUE':
		    initialOrbit = KeplerianOrbit(self.orbit.sma(), self.orbit.eccentricity(),
			self.orbit.inclination, self.orbit.aop, self.orbit.raan, self.orbit.lv,
			PositionAngle.TRUE, eci, startTime, Constants.WGS84_EARTH_MU)
	    elif anomaly == 'MEAN':
		    initialOrbit = KeplerianOrbit(self.orbit.sma(), self.orbit.eccentricity(),
			self.orbit.inclination, self.orbit.aop, self.orbit.raan, self.orbit.lv,
			PositionAngle.MEAN, eci, startTime, Constants.WGS84_EARTH_MU)
	    else:
		    raise ValueError('Unrecognized anomaly type')

	    initialState = SpacecraftState(initialOrbit, self.spacecraft.mass)

	    # Setup numerical propagator
	    minStep = 1e-3
	    maxStep = 1e+3
	    positionTolerance = 1.0
	    orbitType = OrbitType.KEPLERIAN
	    tolerance = NumericalPropagator.tolerances(positionTolerance, initialOrbit, orbitType)

	    # Integrator
	    integrator = DormandPrince853Integrator(minStep, maxStep,
		    JArray_double.cast_(tolerance[0]),
		    JArray_double.cast_(tolerance[1]))
	    integrator.setInitialStepSize(step)

	    propagator = NumericalPropagator(integrator)
	    propagator.setOrbitType(orbitType)

	    # GravityModel
	    gravityModel = GravityFieldFactory.getNormalizedProvider(degree, order)
	    propagator.addForceModel(HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravityModel))

	    ## Perturbations
	    # Third-body Perturbations
	    if (thirdBody):
	    	sun3rdBodyAttraction  = ThirdBodyAttraction(sun)
	    	moon3rdBodyAttraction = ThirdBodyAttraction(moon)
	    	propagator.addForceModel(sun3rdBodyAttraction)
	    	propagator.addForceModel(moon3rdBodyAttraction)
	    # Solar Radiation Pressure
	    if (solarPressure):
		    isotropicRadiationSingleCoeff = IsotropicRadiationSingleCoefficient(self.spacecraft.cross_section, self.spacecraft.cr)
		    solarRadiationPressure = SolarRadiationPressure(sun, Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
			isotropicRadiationSingleCoeff)
		    propagator.addForceModel(solarRadiationPressure)
	    # Atmospheric Drag
	    if (atmosphericDrag):
		    parameters = MarshallSolarActivityFutureEstimation(
			MarshallSolarActivityFutureEstimation.DEFAULT_SUPPORTED_NAMES,
			MarshallSolarActivityFutureEstimation.StrengthLevel.AVERAGE)
		    atmosphere = NRLMSISE00(parameters, sun, earth)
		    isotropicDrag = IsotropicDrag(self.spacecraft.cross_section, self.spacecraft.cd)
		    dragForce = DragForce(atmosphere, isotropicDrag)
		    propagator.addForceModel(dragForce)
	    #Relativity
	    if (relativity):
		    relativityEffect = Relativity(Constants.WGS84_EARTH_MU)
		    propagator.addForceModel(relativityEffect)

	    propagator.setInitialState(initialState)

	    return propagator, startTime, endTime

	def propagate(self, propagator, startTime, endTime, step):#, type):

		while(startTime.compareTo(endTime) <= 0.0):
			pvPosition = propagator.propagate(startTime).getPVCoordinates().getPosition()
			pvVelocity = propagator.propagate(startTime).getPVCoordinates().getVelocity()

			print('DATA pos \"{} {} {}\"'.format(pvPosition.getX()/1e3,
				pvPosition.getY()/1e3,pvPosition.getZ()/1e3))

			startTime = startTime.shiftedBy(step)

sat_list = {
     'SAT': {
         'norad_id': 43226,                     # SpaceTrack ID for TLE
         'cospar_id': '0200901',                # [N/A]
         'sic_id': '6179',                      # [N/A]
         'mass': 5192.0,                        # [kg]
         'cross_section': 13.25,                # [m2]
         'cd': 2.2,                             # Drag coefficient
         'cr': 1.8,                             # SRP coefficient
         'identification': 11111,               # Spacecraft ID
         'origin': 'USA'                        # Country
     }
}
sat_name   = 'SAT'
spacecraft = Spacecraft(sat_name, sat_list[sat_name]['mass'], sat_list[sat_name]['cross_section'],
        sat_list[sat_name]['cd'],sat_list[sat_name]['cr'],
        sat_list[sat_name]['identification'],sat_list[sat_name]['norad_id'])

epoch = datetime(2022, 1, 1, 10, 30, 0)
ra = 814
rp = 786
i =  98.55
omega = 90.0
raan = 5.1917
lv = 0.0567634 
orbit = KeplerOrbit(ra, rp, i, omega, raan, lv)

degree = 70
order  = 70
step   = 60.0   # seconds
duration = 365*24*3600.0 # seconds

keplerianProp = KeplerPropagator(spacecraft, orbit, epoch)
propagator, startTime, endTime = keplerianProp.numericalPropagatorBuilder(duration, step, degree, order,
    anomaly='TRUE',
    thirdBody=True,
    solarPressure = True,
    atmosphericDrag = True,
    relativity = True)
keplerianProp.propagate(propagator, startTime, endTime, step)

Would anyone have an idea of why does this happen?

Many many thanks in advance.

Best regards,
/Leonardo.

2

Hi @lghizoni,

Since your initial eccentricity is small and you propagate on a very long time it’s possible that inside the steps calculation of the integrator the eccentricity becomes negative and the exception gets thrown.
After how long does it happen?
Try propagating in circular elements (use setOrbitType(OrbitType.CIRCULAR) before running the propagator in your propagate function).
And if it works, print the value of e at each step and see if it does not get to close to 0.

Best regards,
Max

1 Like
3

Dear @MaximeJ ,

Thanks a lot for the input! Indeed the propagation works now. I guess I’ll have to add an IF loop to compute eccentricities and change the OrbitType depending on the initial parameters.

Thank you again for your help.

Best regards,
/Leonardo.

4

Hey @MaximeJ , I know this is old, but I am having a similar issue. However if I switch to circular, the estimation portion just goes forever. Any advice?

5

Hi @tru3anomaly,

Hard to say. How big is your measurements’ arc and how many measurements do you have ?
Do you use an OD observer to watch the iterations ?

6

Thanks Maxime. The measurements are somewhat close together, some separated by an hour or so and others by mins. I have a days worth of measurements. Around 30 or so. I was not aware of the OD observer until you told me, I will try it out.