Hey Vincent,
sure I can. Running the following code leads to the infinite loop.
# === INITIALISATION AND IMPORTS ===
import orekit
from orekit.pyhelpers import setup_orekit_curdir
# Start the Java Virtual Machine (JVM) for Orekit
orekit.initVM()
# Load Orekit data from the current working directory
setup_orekit_curdir()
# Import Orekit data
from org.orekit.utils import Constants, IERSConventions
from org.orekit.bodies import OneAxisEllipsoid, GeodeticPoint, CelestialBodyFactory
from org.orekit.time import AbsoluteDate, TimeScalesFactory
from org.orekit.orbits import KeplerianOrbit, PositionAngleType, OrbitType
from org.orekit.frames import FramesFactory, TopocentricFrame
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.generation import FixedPointTleGenerationAlgorithm
from org.orekit.propagation.events import ElevationDetector, EventsLogger, ApsideDetector, DateDetector
from org.orekit.propagation.events.handlers import ContinueOnEvent
from org.orekit.forces.gravity import HolmesFeatherstoneAttractionModel
from org.orekit.forces.gravity.potential import GravityFieldFactory
from org.orekit.forces.maneuvers import ImpulseManeuver
from org.orekit.models.earth.atmosphere import NRLMSISE00, NRLMSISE00InputParameters
from org.orekit.models.earth.atmosphere.data import MarshallSolarActivityFutureEstimation, CssiSpaceWeatherData
from org.orekit.forces.drag import DragForce, IsotropicDrag
from org.orekit.data import DataProvidersManager, DirectoryCrawler
from orekit import JArray_double
# Import hipparchus Data
from org.hipparchus.util import FastMath
from org.hipparchus.ode.nonstiff import DormandPrince853Integrator
from org.hipparchus.geometry.euclidean.threed import Vector3D, RotationConvention, Rotation
# Import standard Python libraries
import math
import pandas as pd
# === FUNCTIONS ===
# Elevation/Contact detector
def elevation_detector(prop):
elevation_detector = (ElevationDetector(station_frame).withConstantElevation(min_elevation)
.withHandler(ContinueOnEvent()))
logger = EventsLogger()
logged_detector = logger.monitorDetector(elevation_detector)
prop.addEventDetector(logged_detector) # Add Event to Propagator
state = prop.propagate(initDate, finalDate) # Detection occurs during this propagation
visible_times = logger.getLoggedEvents()
current_TLE = TLE.stateToTLE(state, my_TLE, FixedPointTleGenerationAlgorithm()) # Convert propagated state to TLE
return visible_times, current_TLE
# Contact duration Calculator
def contact_duration_calc(visible_times):
vis_result = [] # List to store the results of contact duration calculations
start_time = None # initialise start_time
for time in visible_times:
if time.isIncreasing():
start_time = time.getState().getDate() # time.isIncreasing means contact duration starts
# When the event ends (isIncreasing is False), calculate the contact duration
elif start_time: # Only process if start_time has been set
stop_time = time.getState().getDate() # Contact duration ends, so set the stop time
duration = stop_time.durationFrom(start_time) / 60 # Convert seconds to minutes
vis_result.append({'Start': start_time,
'Stop': stop_time,
'Contact duration [min]': duration})
start_time = None # Reset start_time for next contact duration period
return vis_result
# Creates contact duration table
def contact_duration_print(vis_result):
if not vis_result: # If no result logged
print("No events were logged. Check the event detection logic.")
else:
# Convert list of results to a DataFrame for better formatting
result_df = pd.DataFrame.from_dict(vis_result)
# Limit 'Start' and 'Stop' columns to the first 21 characters for better readability
result_df['Start'] = result_df['Start'].astype(str).str[:21] # Limits column width to 21 characters
result_df['Stop'] = result_df['Stop'].astype(str).str[:21]
pd.set_option('display.max_columns', None) # Display alle columns in dataframe
pd.set_option('display.max_colwidth', 25) # Set limit on column width
print(f"{result_df}")
# === FURTHER DATA ===
# CONSTANTS
r_earth = Constants.WGS84_EARTH_EQUATORIAL_RADIUS
mu_earth = Constants.WGS84_EARTH_MU
g_0 = Constants.G0_STANDARD_GRAVITY
# BODIES
sun = CelestialBodyFactory.getSun()
earth = OneAxisEllipsoid(r_earth,
Constants.WGS84_EARTH_FLATTENING,
FramesFactory.getITRF(IERSConventions.IERS_2010, True))
# FRAME & TIMEZONE
frame = FramesFactory.getTEME()
utc = TimeScalesFactory.getUTC()
# GROUND STATION
longitude = FastMath.toRadians(9.1038)
latitude = FastMath.toRadians(48.7495)
altitude = 0.0
station = GeodeticPoint(latitude, longitude, altitude) # Position of ground station
station_frame = TopocentricFrame(earth, station, "Station") # Frame for further calculations
min_elevation = FastMath.toRadians(5.0) # Minimum ground station elevation angle
# ENGINE
bd = 6 # burn duration
#initMass = float(input("Enter initial satellite mass [kg]:"))
initMass = 82.5
mass = initMass # initial mass
mf = 0.0046 # mass flow
isp_ideal = 360.0 # ideal specific impulse
# DRAG
cd = 2.2
area = 0.64
# THRUST DEVIATION
isp = []
isp.append(isp_ideal * 0.97)
isp.append(isp_ideal) # Thrust is directly proportional to Isp, that's why thrust deviations are adjusted via Isp
isp.append(isp_ideal * 1.03)
# THRUST VECTOR MISALIGNMENT
opening_angle_deg = 3.53 # maximum deviation of thrust direction
step_angle_deg = 45.0 # iteration step size for conical vectors (minimum 5°)
opening_angle_rad = FastMath.toRadians(opening_angle_deg) # convert to rad
step_angle_rad = FastMath.toRadians(step_angle_deg) # convert to rad
# DEFINE TIME SPAN TO PROPAGATE
t_shifted = 60.0 * 60 * 24 # Elapsed time in seconds, default value: 24 hours
# reference date
year = 2026
month = 10
day = 1
hour = 0
minute = 0
second = 0.0
# Initial and final Date of propagation
initDate = AbsoluteDate(year, month, day, hour, minute, second, utc)
finalDate = initDate.shiftedBy(t_shifted)
a = 7000 * 1000.0
e = 0.13
incl = 98.6
raan = 90.0
omega = 100.0
ta = 20.0
# Create orbit
orbit = KeplerianOrbit(a, e, incl, omega, raan, ta, PositionAngleType.TRUE, frame,
initDate, mu_earth)
# Set up propagator structure with default settings
minStep = 0.001
maxStep = 500.0
initStep = 60.0
positionTolerance = 0.001 # [m]
# tolerances[0]: maximum absolute error, tolerances[1]: maximum relative error
tolerances = NumericalPropagator.tolerances(positionTolerance, orbit, orbit.getType())
# DormandPrince is an 8th-order integrator with variable step size
integrator = DormandPrince853Integrator(minStep, maxStep, JArray_double.cast_(tolerances[0]),
JArray_double.cast_(tolerances[1]))
initialState = SpacecraftState(orbit, initMass) # Create initial spacecraft state using orbit configuration
integrator.setInitialStepSize(initStep)
propagator = NumericalPropagator(integrator) # Create numerical propagator with selected integrator
propagator.setOrbitType(OrbitType.CARTESIAN) # Required because ForceModel operates with cartesian coordinates
propagator.setInitialState(initialState)
# gravity force model
gravityProvider = GravityFieldFactory.getNormalizedProvider(10, 10)
ForceModel_gravity = HolmesFeatherstoneAttractionModel(FramesFactory.getITRF(IERSConventions.IERS_2010, True),
gravityProvider)
propagator.addForceModel(ForceModel_gravity)
#atmosphere force model
atmosphere = NRLMSISE00(CssiSpaceWeatherData(CssiSpaceWeatherData.DEFAULT_SUPPORTED_NAMES), sun, earth)
drag_sensitive = IsotropicDrag(area, cd)
ForceModel_drag = DragForce(atmosphere, drag_sensitive)
propagator.addForceModel(ForceModel_drag)
# Convert initial orbit representation to Keplerian coordinates
OrbitType.KEPLERIAN.convertType(propagator.getInitialState().getOrbit())
# The TLE is used here to enable a stateToTLE conversion later
# Therefore, this TLE is only needed to gather information about the satellite ID
tle_file = "/Users/maxgantze/Desktop/BA/flp17016.tle" # add file path of TLE
with open(tle_file, "r") as f: # open file in read mode
lines = f.readlines() # Read all lines into a list
tle_line1 = lines[1].strip() # First line of TLE (remove extra spaces/newlines)
tle_line2 = lines[2].strip() # Second line of TLE
my_TLE = TLE(tle_line1, tle_line2) # Create TLE object
# Print dates
print(f"\nInitial Date: {initDate}")
print(f"Final Date: {finalDate}")
first_event = None # Used only for type identification of first event (Peri/Apo)
first_detection = None # Ensures that the Peri/Apo list is only printed once
# Ideal propagation
visible_times_ideal, current_TLE_ideal = elevation_detector(propagator)
vis_result_ideal = contact_duration_calc(visible_times_ideal)
print("\n\nContact duration table, ideal:\n")
contact_duration_print(vis_result_ideal)
# Display final TLE after propagation in ideal scenario
print(f"\n\nTLE at final Date, ideal:\n{current_TLE_ideal}")
Thank you once again!
Cheers,
Max