Hello all,
I am building a propagator for LEO, MEO, and GEO orbits. As such, I have designed it to use Earth harmonics (120 degrees), SRP, atmospheric drag, and point masses for Earth, Moon, and Sun.
Over the course of 9 hours, I have an residual of ~168m; however, a colleague of mine using Java Orekit got an error of ~50m with the same input parameters. He states that it’s setup with the exact same perturbing forces and numerical integrator bounds. Unfortunately, he cannot share his code with me.
TLDR: I need help tracking down a small source of error that could throw a GEO satellite propagation off by 100m over the course of 9 hours. I’ve attached my code. Thank you!
import numpy as np
import pandas as pd
import time as time
from datetime import datetime, timedelta
#import ipywidgets
start_time = time.time()
#%matplotlib widget
import matplotlib.pyplot as plt
import orekit
from orekit.pyhelpers import setup_orekit_curdir
orekit.initVM()
setup_orekit_curdir()
from org.orekit.frames import FramesFactory
from org.orekit.orbits import CartesianOrbit, PositionAngleType # Import CartesianOrbit
from org.orekit.time import AbsoluteDate, TimeScalesFactory
from org.orekit.utils import Constants, DoubleArrayDictionary, PVCoordinates
from org.orekit.propagation import SpacecraftState
from org.orekit.propagation.numerical import NumericalPropagator
from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel
from org.orekit.forces.gravity.potential import GravityFieldFactory
from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator
from orekit import JArray_double
from org.orekit.orbits import OrbitType
from org.orekit.utils import IERSConventions
from org.orekit.forces.gravity import NewtonianAttraction
from org.orekit.forces.gravity import ThirdBodyAttraction
from org.orekit.bodies import CelestialBodyFactory
from org.orekit.models.earth import ReferenceEllipsoid
from org.orekit.forces.radiation import SolarRadiationPressure, IsotropicRadiationSingleCoefficient
from org.orekit.bodies import OneAxisEllipsoid
# Set up initial orbit FOR NORAD 22314
Cross_Area = 13.873
cd = 0.0
satellite_mass = 2511.4 # The models need a spacecraft mass, unit kg.
solar_coeff = 0.0345417
# 1. Define initial state as a state vector (Cartesian Position and Velocity)
utc = TimeScalesFactory.getUTC()
initial_date = AbsoluteDate(2024, 6, 17, 22, 8, 42.425880, utc)
final_date = AbsoluteDate(2024, 6, 18, 7, 27, 12.129134, utc)
delta_t_seconds = final_date.durationFrom(initial_date) # date2 - date1
print(f"Time difference: {delta_t_seconds} seconds")
inertial_frame = FramesFactory.getEME2000() # Or GCRF, J2000, etc.
wgs84Ellipsoid = ReferenceEllipsoid.getWgs84(inertial_frame)
#Rotating Reference Earth
# Preferred: ITRF for WGS84 ellipsoid
itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
Constants.WGS84_EARTH_FLATTENING,
itrf)
# Define your initial state vector components
# Example: GEO satellite
# Create Vector3D objects for position and velocity
initial_position = Vector3D(-41105.92780995852e3, -9050.561409355181e3, -1809.49929235732e3)
initial_velocity = Vector3D(0.671176665193438e3,-2.910029484561473e3, -0.741988755546521e3)
# Final Position vectors
Truth_Positions = np.array([37459.41326275722e3 ,-18756.77890229113e3 ,-5164.26162320142e3])
Truth_Velocities = np.array([1.414999564621974e3, 2.64689489090523e3,0.651570654272371e3])
# Create PVCoordinates (Position-Velocity Coordinates)
initial_pv = PVCoordinates(initial_position, initial_velocity)
# Create CartesianOrbit from PVCoordinates
mu = Constants.WGS84_EARTH_MU # Gravitational parameter of Earth
initial_orbit = CartesianOrbit(initial_pv, inertial_frame, initial_date, mu)
# Define spacecraft mass (e.g., 1000.0 kg)
initial_state = SpacecraftState(initial_orbit,satellite_mass)
# Set up numerical integrator
minStep = 0.0001
maxStep = 1000.0
initStep = 60.0
positionTolerance = 1.0
tolerances = NumericalPropagator.tolerances(positionTolerance,
initial_orbit,
initial_orbit.getType())
integrator = DormandPrince853Integrator( minStep, maxStep,
10e-14,
10e-12)
propagator_num = NumericalPropagator(integrator)
propagator_num.setOrbitType(OrbitType.CARTESIAN)
propagator_num.setInitialState(initial_state)
##### EARTH GRAVITY + HARMONICS
gravityProvider = GravityFieldFactory.getNormalizedProvider(120, 120)
propagator_num.addForceModel(HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, True), gravityProvider))
propagator_num.addForceModel(NewtonianAttraction(mu))
# Add third-body perturbations
sun = CelestialBodyFactory.getSun()
moon = CelestialBodyFactory.getMoon()
propagator_num.addForceModel(ThirdBodyAttraction(sun))
propagator_num.addForceModel(ThirdBodyAttraction(moon))
# Add Atmospheric Drag
from org.orekit.models.earth.atmosphere.data import CssiSpaceWeatherData
cswl = CssiSpaceWeatherData("SpaceWeather-All-v1.2.txt")
from org.orekit.models.earth.atmosphere import NRLMSISE00
atmosphere = NRLMSISE00(cswl, sun, earth)
from org.orekit.forces.drag import IsotropicDrag
isotropic_drag = IsotropicDrag(Cross_Area, cd)
from org.orekit.forces.drag import DragForce
drag_force = DragForce(atmosphere, isotropic_drag)
propagator_num.addForceModel(drag_force)
# ADD SRP
srp = SolarRadiationPressure(
CelestialBodyFactory.getSun(),
earth,
IsotropicRadiationSingleCoefficient(Cross_Area, solar_coeff)
)
propagator_num.addForceModel(srp)
############# PROPAGATION
duration = delta_t_seconds
step = 60.0
steps = int(duration / step)
times = np.linspace(0, duration, steps)
positions = []
velocities = []
for t in times:
state = propagator_num.propagate(initial_date.shiftedBy(float(t)))
pvT = state.getPVCoordinates()
positions.append(pvT.getPosition().toArray())
velocities.append(pvT.getVelocity().toArray())
positions = np.array(positions)
last_positions = positions[-1]
last_velocities = velocities[-1]
Poss_diff = last_positions - Truth_Positions
Poss_diff_norm = np.linalg.norm(Poss_diff)
Vel_diff = last_velocities - Truth_Velocities
Vel_diff_norm = np.linalg.norm(Vel_diff)
print(Poss_diff)
print(Poss_diff_norm)
print(Vel_diff)
print(Vel_diff_norm)