orekit_hohmann_mod.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("orekit-data/")


# 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, EquinoctialOrbit
from org.orekit.propagation.analytical import KeplerianPropagator
from org.orekit.propagation.events import PositionAngleDetector, ApsideDetector
from org.orekit.propagation.events.handlers import StopOnDecreasing, StopOnIncreasing, StopOnEvent
from org.orekit.time import AbsoluteDate, TimeScalesFactory
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
utc = TimeScalesFactory.getUTC()
# r0_vec = Vector3D(float(7200e3), float(-1000e3), float(0.0))  # [m]
r0_vec = Vector3D(float(7200.000000990574e3), float(-999.9999928676314e3), float(-0.0006883381606922735e3)) # GMAT report init values *1e3
# v0_vec = Vector3D(float(0), float(8e3), float(0))  # [m / s]
v0_vec = Vector3D( float(-7.924862975769942e-06), float(7.999999999999999e3), float(-9.335036198378994e-05)) 
rv_0 = PVCoordinates(r0_vec, v0_vec)
epoch_00 = AbsoluteDate(2000, 1, 1, 11, 59, 28.0, utc)
gcrf_frame = FramesFactory.getGCRF()

ss_0 = KeplerianOrbit(rv_0, gcrf_frame, epoch_00, k)
# ss_0 = EquinoctialOrbit(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)

# print(f"{dv_a=}, {dv_b=}")
# GMAT values
dv_a = 1.47424834733 * 1e3
dv_b = 1.44567622075 * 1e3

# print(f"{dv_a=}, {dv_b=}")

# print(f"{dv_a=} {dv_b=}")
# 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)
# )
at_apoapsis = ApsideDetector(ss_0).withHandler(StopOnDecreasing()) # default apside detector stops on increasing g function (aka at perigee)
at_periapsis = ApsideDetector(ss_0)
# 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]

# explicit tof from GMAT report timestamp
# expected_tof = AbsoluteDate(2000, 1, 1, 17, 29, 18.489, utc).durationFrom(epoch_00)

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),
]

# my run on gmat
rf_expected = [-18447.18133044379e3,-30659.53028128713e3, 0.002151253568205197e3 ]

vf_expected = [
    float(2.859891379982459e3),
    float(-1.720738617222231e3),
    float(1.600645301237219e-14),
]

# my run on gmat
vf_expected = [2.859891380100071e3,-1.720738616179917e3,-2.579688114910805e-07*1e3    ]



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