Modeling LEO Non-Impulsive Orbit Raise

Orekit Python Wrapper usage
1

Hello everyone,

I am brand new to Orekit and have been the exploring the wrapper docs and examples to research a task I’d like to accomplish to compare with another analysis toolkit. I am interested in raising an initial LEO Orbit 1 to a higher altitude Orbit 2 (~100km) using a Constant Thrust (Non-Impulsive) Electric burn. I hope to specify parameters for Vehicle Mass and Isp/flow rate values for various e-propulsive candidates to reach the final orbit.

I found a couple of examples (below) of the ImpulseManeuver and ConstantThrustManeuver class but from what I can see, the latter is defined as f(t0, duration, thrust, isp, direction) with duration and direction but not to a specific objective/target Orbit 2 (in my case a higher altitude in same plane). Could someone please illustrate a potential use case to non-impulsively maneuver between two pre-defined orbits perhaps using one of the event detector(s) at the final altitude instead of being constrained by a defined burn duration?

TYIA!

2

Orekit Python Users,
I attempted this query for a ConstantThrust Maneuver application to a higher altitude target orbit in the below. I’m wondering what I’m doing wrong since I’m not achieving an orbit raise with this setup. Any help would be very appreciated!

  1. What is the proper method to apply the altitude detector trigger to stop a maneuver burn
  2. Is it at all possible to omit constraining on Thrust burn duration and/or Delta-V and instead allow the propagation find a solution which reaches the intended final orbit? Could a different propagator type resolve this issue since I tried the same with Keplerian propagator with no success.
# constants
inertial_frame = FramesFactory.getEME2000()
attitude = LofOffset(inertial_frame, LOFType.LVLH)
EARTH_RADIUS = Constants.WGS84_EARTH_EQUATORIAL_RADIUS
G0_STANDARD_GRAVITY = Constants.G0_STANDARD_GRAVITY
MU = Constants.WGS84_EARTH_MU
UTC = TimeScalesFactory.getUTC()
PI = np.pi

##################################################
# INITIAL Orbit, 950km
a = 7328137.0
ex = 1.0e-6
ey = 1.0e-6
i = radians(80.0)
raan = radians(0.0)
initialOrbit = CircularOrbit(a, ex, ey, i, raan, 0.0, 
                   PositionAngle.MEAN,
                   inertial_frame,
                   epoch_00,
                   MU)
mass_wet = 5000.0
initialState = SpacecraftState(initialOrbit, 
                    attitude.getAttitude(initialOrbit, initialOrbit.getDate(), initialOrbit.getFrame()), mass_wet)

# Targeted FINAL Orbit (e.g. 10km higher)
ra = 960 * 1000 + EARTH_RADIUS # m Apogee
rp = 960 * 1000 + EARTH_RADIUS # m Perigee
a = (rp + ra) / 2.0
e = 1.0 - (rp / a)
i = radians(80.0)
argperigee = radians(0.0)
raan = radians(0.0)
anomaly = radians(0.0) #mean, eccentric or true anomaly
finalOrbit = KeplerianOrbit(a, e, i, argperigee, raan, anomaly,
                                PositionAngle.TRUE,
                                inertial_frame,
                                epoch_00,
                                MU)

# Generate a numerical propagator and add events
minStep = 0.001
maxstep = 1000.0
initStep = 10.0
positionTolerance = 1.0
orbitType = initialOrbit.getType()
tolerances = NumericalPropagator.tolerances(positionTolerance, initialOrbit, orbitType)
integrator = DormandPrince853Integrator(minStep, maxstep, 
                                        JArray_double.cast_(tolerances[0]),
                                        JArray_double.cast_(tolerances[1]))
integrator.setInitialStepSize(initStep)
propagator = NumericalPropagator(integrator)
propagator.setInitialState(initialState)
propagator.setOrbitType(orbitType)
propagator.setMu(MU)
propagator.setAttitudeProvider(attitude)

# Add force model gravity field
newattr = NewtonianAttraction(MU)
itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, True)  # International Terrestrial Reference Frame, earth fixed
earth = OneAxisEllipsoid(EARTH_RADIUS,
                        Constants.WGS84_EARTH_FLATTENING,itrf)
gravityProvider = GravityFieldFactory.getNormalizedProvider(8, 8)
propagator.addForceModel(HolmesFeatherstoneAttractionModel(earth.getBodyFrame(), gravityProvider))

# Define ConstantThrustManeuver and addForceModel
# date - maneuver date
# duration - the duration of the thrust (s) (if negative, the date is considered to be the stop date)
# thrust - the thrust force (N)
# isp - engine specific impulse (s)
# direction - the acceleration direction in satellite frame.
#? Do we have to specify propagation time or can it forced to stop at target altitude?
#? Same as above for omit providing a dV and stop at altitude?
t_maneuver = initialOrbit.getDate().shiftedBy(1000.0)
t_prop = t_maneuver.shiftedBy(10 * initialOrbit.getKeplerianPeriod())  
dV = Vector3D(10.0, Vector3D.PLUS_I)

f_thrust = 20.0  #N
isp = 300.0   # seconds
vExhaust = G0_STANDARD_GRAVITY * isp
mass = initialState.getMass()
dt = -(mass * vExhaust / f_thrust) * FastMath.expm1(-dV.getNorm() / vExhaust)
maneuver = ConstantThrustManeuver(t_maneuver, dt, f_thrust, isp, dV.normalize())
print(f"FiniteThrustManeuver: {dt = :.2f} seconds, f_thrust = {f_thrust:.1f} N, {isp = :.1f} seconds, dV = {dV.getNorm():.3f} m/s")
propagator.addForceModel(maneuver)

# Define AltitudeDetector at the Target/Final Orbit
earth = OneAxisEllipsoid(EARTH_RADIUS, 
                        Constants.WGS84_EARTH_FLATTENING, 
                        inertial_frame)
target_alt_m = finalOrbit.getPVCoordinates().getPosition().getNorm() - EARTH_RADIUS # (m)
altDetector = AltitudeDetector(target_alt_m,  earth).withHandler(StopOnIncreasing().of_(AltitudeDetector))  #stop propagation on ascending altitude
print(f"target altitude: {altDetector.getAltitude()} meters")
propagator.addEventDetector(altDetector)
# propagator.setEphemerisMode()

finalState = propagator.propagate(t_prop)
finalPropTime = finalState.getDate()
print(f"propagation limited to {t_prop}\n")
print(finalState)
print(f"{finalPropTime = }")
finalAlt = finalState.getPVCoordinates().getPosition().getNorm() - EARTH_RADIUS
finalStateEqualsAlt = abs(finalAlt - target_alt_m) < 1e-5
print(f"{finalStateEqualsAlt = } , {finalAlt = }")
# finalStateEqualsAlt = False , finalAlt = 946938.5746795079
3

Hi @conz1,

Welcome to Orekit and its forum! And sorry for the late answer.

I haven’t tried your code but will (or its equivalent in Java at least) once I have time.
I’ve read it very quickly and, unless I missed something, you’re maneuvering along the X axis of the LVLH LOF frame which is actually the radial axis.
That may explain why you never reach the target altitude.
Perhaps you could try with the Y axis instead (velocity axis), and also print the Keplerian orbit at regular steps to see how it evolves with time during the maneuver.

Perhaps you could try your luck with the EventBasedManeuverTriggers instead of using the DateBasedManeuverTriggers that is built in the ConstantThrustManeuver constructor you are using. The starting event would be a DateDetector and the ending event would be the altitude detector.
I’ve never personally tried it though.

Not yet, unfortunately there is no maneuver computation/optimization implementation at the moment in Orekit.

It would be a very good addition to the library and we would like to have it; but it’s quite an important piece of work.

Contributions are always welcome if you feel like getting involved in the project :wink:

I wonder if we couldn’t use a “Fielded” version of the Maneuver model (i.e. a version where we would compute the nth-order derivatives along the propagation) and invert the Taylor algebra (i.e. find the thrust, direction and duration that allow reaching the desired orbit) to compute the maneuver needed to reach a target from a starting point. But it would be a two-fold process anyway, a first propagation intended to find the proper maneuver, and a second one to actually realize it.

Have a good day,
Max