orekit_hohmann.py

""" Modeling Hohmann impulsive transfer with Orekit Python wrapper """

import numpy as np

# Setup the virtual machine and data location
import orekit
from orekit.pyhelpers import setup_orekit_curdir

vm = orekit.initVM()

# Data orekit-data.zip must be located at same dir/ level that this script
setup_orekit_curdir("data/orekit-data.zip")

# Import different utilities from Orekit API
from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.orekit.attitudes import LofOffset
from org.orekit.forces.maneuvers import ImpulseManeuver
from org.orekit.frames import FramesFactory, LOFType
from org.orekit.orbits import KeplerianOrbit, OrbitType, PositionAngle
from org.orekit.propagation.analytical import KeplerianPropagator
from org.orekit.propagation.events import PositionAngleDetector, ApsideDetector
from org.orekit.time import AbsoluteDate
from org.orekit.utils import Constants as C
from org.orekit.utils import PVCoordinates

# Build the initial orbit
k = C.IAU_2015_NOMINAL_EARTH_GM
r0_vec = Vector3D(float(7200e3), float(-1000e3), float(0))  # [m]
v0_vec = Vector3D(float(0), float(8e3), float(0))  # [m / s]
rv_0 = PVCoordinates(r0_vec, v0_vec)
epoch_00 = AbsoluteDate.J2000_EPOCH
gcrf_frame = FramesFactory.getGCRF()

ss_0 = KeplerianOrbit(rv_0, gcrf_frame, epoch_00, k)
print(ss_0, end="\n\n")

# Final desired orbit radius and auxiliary variables
rf_norm = float(35781.34857e3)  # [m]
r0_norm, v0_norm = r0_vec.getNorm(), v0_vec.getNorm()
a_trans = (r0_norm + rf_norm) / 2

# Compute the magnitude of Hohmann's deltaV
dv_a = np.sqrt(2 * k / r0_norm - k / a_trans) - v0_norm
dv_b = np.sqrt(k / rf_norm) - np.sqrt(2 * k / rf_norm - k / a_trans)

# Convert previous magnitudes into vectors within the TNW frame, which
# is aligned with local velocity vector
dVa_vec, dVb_vec = [Vector3D(float(dV), float(0), float(0)) for dV in (dv_a, dv_b)]

# Local orbit frame: X-axis aligned with velocity, Z-axis towards momentum
attitude_provider = LofOffset(gcrf_frame, LOFType.TNW)

# Define the events when that trigger the impulses: first one happens at
# periapsis (true anomaly = 0 [deg]) while the second does at apoapsis (true anomaly =
# 180 [deg])
at_periapsis = PositionAngleDetector(OrbitType.KEPLERIAN, PositionAngle.TRUE, float(0))
at_apoapsis = PositionAngleDetector(
    OrbitType.KEPLERIAN, PositionAngle.TRUE, float(np.pi)
)

# Build the impulsive maneuvers
ISP = float(300)
imp_A = ImpulseManeuver(at_periapsis, attitude_provider, dVa_vec, ISP)
imp_B = ImpulseManeuver(at_apoapsis, attitude_provider, dVb_vec, ISP)

# Generate the propagator and add the maneuvers
propagator = KeplerianPropagator(ss_0, attitude_provider)
propagator.addEventDetector(imp_A)
propagator.addEventDetector(imp_B)

# Propagate the orbit and retrieve final state vectors. Force vectors to be in
# original frame (GCRF), not in (TNW)
expected_tof = float(0.229056589 * 24 * 3600)  # [day * h/day * s/h = s]
expected_tof = float(19790.48929209821)  # [day * h/day * s/h = s]
epoch_ff = epoch_00.shiftedBy(expected_tof)
rv_f = propagator.propagate(epoch_ff).getPVCoordinates(gcrf_frame)
rf_vec, vf_vec = (
    (list(rv_f.getPosition().toArray())),
    (list(rv_f.getVelocity().toArray())),
)

# Check if final vectors match expected ones predicted by GMAT-NASA 2020a
rf_expected = [
    float(-18447.18134824541e3),
    float(-30659.53026513587e3),
    float(2.116148304903275e-10),
]
vf_expected = [
    float(2.859891379982459e3),
    float(-1.720738617222231e3),
    float(1.600645301237219e-14),
]

print(f"Assert final position\n{rf_vec = }\n{rf_expected = }\n")
print(f"Assert final velocity\n{vf_vec = }\n{vf_expected = }")