Not realistic EclipseDetector events for SC on a SSO

Hello, everyone!
I am a new user to the Orekit library. I have been trying to determine the eclipse timings for a SC on a SSO orbit using EclipseDetector. However, the eclipse periods seem to vary with DOY from 35 minutes in duration to almost no eclipses. I have looked at the tutorials and the posts here, but my eclipse crossings still look unnatural. I am attaching the code that I was used and the produced plots. I would appreciate your help. Thanks!

import orekit
vm = orekit.initVM()
print ('Java version:',vm.java_version)
print ('Orekit version:', orekit.VERSION)

from orekit.pyhelpers import setup_orekit_curdir, absolutedate_to_datetime
setup_orekit_curdir()
from org.hipparchus.geometry.euclidean.threed import Vector3D
from org.orekit.data import DataProvidersManager, ZipJarCrawler
from org.orekit.propagation import SpacecraftState 
from org.orekit.propagation.analytical import KeplerianPropagator, EcksteinHechlerPropagator
from org.orekit.orbits import KeplerianOrbit, CartesianOrbit, PositionAngle, OrbitType
from org.orekit.frames import FramesFactory, TopocentricFrame
from org.orekit.bodies import OneAxisEllipsoid, GeodeticPoint, CelestialBodyFactory
from org.orekit.time import TimeScalesFactory, AbsoluteDate, DateComponents, TimeComponents
from org.orekit.utils import IERSConventions, Constants, PVCoordinates, PVCoordinatesProvider
from org.orekit.propagation.analytical.tle import TLE, TLEPropagator
from org.orekit.propagation.events import EclipseDetector, EventsLogger, GeographicZoneDetector
from org.orekit.propagation.events.handlers import ContinueOnEvent, RecordAndContinue, StopOnDecreasing, StopOnEvent, StopOnIncreasing
from math import radians, pi
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import datetime
import cartopy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER
import matplotlib.pyplot as plt
import plotly
import plotly.express as px
import plotly.io as pio

pio.renderers.default='browser'


tle_line1 = "1 43783U 18099AB  23093.86347468  .00002220  00000-0  20937-3 0  9995"
tle_line2 = "2 43783  97.5888 155.3721 0008405  30.4825 329.6883 14.94739059236174"
mytle = TLE(tle_line1,tle_line2)
utc = TimeScalesFactory.getUTC()
inertialFrame = FramesFactory.getEME2000()

#Overall duration in seconds for extrapolation
duration = 24 * 60 * 60
step_time = 1
h=600000.0
Re=Constants.WGS84_EARTH_EQUATORIAL_RADIUS
initialDate = AbsoluteDate(2023, 4, 1, 20, 43, 24.000, utc)
finalDate=initialDate.shiftedBy(float(duration))
mu=Constants.WGS84_EARTH_MU
ITRF = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
earth = OneAxisEllipsoid(Re, 
                          Constants.WGS84_EARTH_FLATTENING, 
                         ITRF)


sun = CelestialBodyFactory.getSun()     # Here we get it as an CelestialBody
sun = PVCoordinatesProvider.cast_(sun)  # But we want the PVCoord interface
R_sun=Constants.SUN_RADIUS
p0=Vector3D(1367598.478,6842810.7,-3480.816)
v0=Vector3D(1022.636, 200.575,-7485.674)

c_orbit=CartesianOrbit(PVCoordinates(p0, v0), inertialFrame,initialDate,mu)
tle2=TLE.stateToTLE(SpacecraftState(c_orbit),mytle)
print(tle2)
print ('Epoch :',tle2.getDate())
propagator = TLEPropagator.selectExtrapolator(tle2)
eclipse_detector = EclipseDetector(sun, R_sun, earth).withPenumbra().withHandler(ContinueOnEvent()).withMaxCheck(60.0) 

logger = EventsLogger()
logged_detector = logger.monitorDetector(eclipse_detector)
propagator.addEventDetector(logged_detector)
propagator.propagate(initialDate,finalDate)
events = logger.getLoggedEvents()
logger.clearLoggedEvents()
print('# of eclipse events:', events.size())
#Time array in orekit AbsoluteDate format

t = [initialDate.shiftedBy(float(dt)) \
        for dt in np.arange(0, duration, step_time)]
tle_pv=[propagator.propagate(tt).getPVCoordinates(inertialFrame) for tt in t]


p_tle = [tpv.getPosition() for tpv in tle_pv]
v_tle= [tpv.getVelocity() for tpv in tle_pv]

prop_data = pd.DataFrame(data=tle_pv,columns=['pv'])
prop_data['Position'] = prop_data['pv'].apply(lambda x: x.getPosition())
prop_data['x_ECI'] = prop_data['Position'].apply(lambda x: x.getX())
prop_data['y_ECI'] = prop_data['Position'].apply(lambda x: x.getY())
prop_data['z_ECI'] = prop_data['Position'].apply(lambda x: x.getZ())
prop_data['Velocity'] = prop_data['pv'].apply(lambda x: x.getVelocity())
prop_data['vx_ECI'] = prop_data['Velocity'].apply(lambda x: x.getX())
prop_data['vy_ECI'] = prop_data['Velocity'].apply(lambda x: x.getY())
prop_data['vz_ECI'] = prop_data['Velocity'].apply(lambda x: x.getZ())
prop_data['datetime'] = prop_data['pv'].apply(lambda x: absolutedate_to_datetime(x.getDate()))
prop_data.set_index('datetime', inplace=True, drop=False)
prop_data.index.name = 'Timestamp'
prop_data['groundpoint'] = prop_data['pv'].apply(lambda pv: earth.transform(pv.position, inertialFrame, pv.date))
prop_data['latitude'] = np.degrees(prop_data.groundpoint.apply(lambda gp: gp.latitude))
prop_data['longitude'] = np.degrees(prop_data.groundpoint.apply(lambda gp: gp.longitude))
prop_data['day'] = prop_data.datetime.dt.dayofyear
prop_data['hour'] = prop_data.datetime.dt.hour
prop_data['minute'] = prop_data.datetime.dt.minute

subpoint_tle = [earth.transform(tp, inertialFrame, tt)  for tt, tp in zip(t,p_tle)]
lat_tle = np.degrees([gp.getLatitude()  for gp in subpoint_tle])
lon_tle = np.degrees([gp.getLongitude() for gp in subpoint_tle])

start_time = None
result = []

for event in events:
    if not event.isIncreasing():
        start_time = event.getState().getDate()
    elif start_time:
        stop_time = event.getState().getDate()
        # result.append({"Start_iso":absolutedate_to_datetime(start_time).isoformat(' ', 'seconds'), 
        #                "Start_none":absolutedate_to_datetime(start_time),
        #                'Start_micros':absolutedate_to_datetime(start_time).replace(tzinfo=None,microsecond=0),})
        result.append({"Start":absolutedate_to_datetime(start_time).replace(tzinfo=None,microsecond=0), 
                        "Stop":absolutedate_to_datetime(stop_time).replace(tzinfo=None,microsecond=0),
                        'Start_abs':start_time,
                        'Stop_abs':stop_time,
                        "EclipseDuration": stop_time.durationFrom(start_time)})
                        # "EclipseDuration": stop_time.durationFrom(start_time)})
        start_time = None
eclipse_df = pd.DataFrame.from_dict(result)
result = []

for event in events:
    sign_event=event.isIncreasing()
    start_time = event.getState().getDate()
    diff=absolutedate_to_datetime(start_time).replace(tzinfo=None,microsecond=0)-absolutedate_to_datetime(initialDate)
    result.append({"Start":absolutedate_to_datetime(start_time).replace(tzinfo=None,microsecond=0), 
                        "isincreasing":sign_event,
                        'seconds':round(diff.total_seconds()),
                        })

    start_time = None
eclipse_event = pd.DataFrame.from_dict(result)

plt.figure(8)

ax = plt.axes(projection=ccrs.PlateCarree())
ax.stock_img()
ax.coastlines()
plt.scatter(prop_data.longitude,prop_data.latitude,
          color='b',s=1,
          # transform=ccrs.Geodetic(),
          transform=ccrs.PlateCarree(),
          )
for i in range(len(eclipse_event.Start)-1):
    if eclipse_event.isincreasing[i]:
        plt.scatter(lon_tle[eclipse_event.seconds[i]], lat_tle[eclipse_event.seconds[i]],
          color='r',s=3,
          # transform=ccrs.Geodetic(),
          transform=ccrs.PlateCarree(),
          )
    else:
        plt.scatter(lon_tle[eclipse_event.seconds[i]:eclipse_event.seconds[i+1]], lat_tle[eclipse_event.seconds[i]:eclipse_event.seconds[i+1]],
              color='y',s=3,
              # transform=ccrs.Geodetic(),
              transform=ccrs.PlateCarree(),
              )
if not eclipse_event.isincreasing.iloc[-1]:
    plt.scatter(lon_tle[eclipse_event.seconds.iloc[-1]:-1], lat_tle[eclipse_event.seconds.iloc[-1]:-1],
          color='y',s=3,
          # transform=ccrs.Geodetic(),
          transform=ccrs.PlateCarree(),
          )
else:
    plt.scatter(lon_tle[eclipse_event.seconds.iloc[-2]:eclipse_event.seconds.iloc[-1]], lat_tle[eclipse_event.seconds.iloc[-2]:eclipse_event.seconds.iloc[-1]],
          color='y',s=3,
          # transform=ccrs.Geodetic(),
          transform=ccrs.PlateCarree(),
          )
fig = plt.gcf()
fig.set_size_inches(18.5, 10.5)
plt.show()

Hi @spaceman,

Welcome to the forum !!

I’ve quickly looked at the Orekit part of your code and didn’t see anything that seemed wrong.
Still I’m wondering why you’re using a Cartesian orbit then convert it using TLE.stateToTLE (which is a not very accurate) instead of using the original TLE at the top of your code ?

Could you please elaborate a bit on that? Do you know the local time at ascending node of your SSO orbit ?
Can you see a seasonal pattern on when there is no eclipse and when it is the longest?
Maybe the above conversion introduced errors that made your orbit drift from the SSO (a,i) combination?

1 Like

Hello and thank you, @MaximeJ!

I wanted to use SGP4 propagator in my simulations, hence the generation of custom TLE. The original TLE (mytle) is used as a template for the generation of a new custom TLE (tle2) that is based on the [p0,v0] state vector. I have generated the initial state vector using STK orbit wizard’s SSO orbit (I believe the local time was either 22:30 or 10:30) and altitude 600km.

I do see seasonal patterns:

  1. closer to March 22 the eclipse duration shortens,
  2. closer to June 22/December 22 it is at max.

Starting from a position/velocity pair (which are osculating parameters) to generate a TLE is error-prone.
You should perform a fitting on several orbits to capture the proper dynamics. This is just the dual part of avoiding to convert from TLE to osculating elements. By the way, you can also use Orekit phasing tutorial to generate your orbit (but again, it will be an osculating one, except if you modify the code to append a fitting at the end).

Using a TLE for long term prediction (several months it seems as you cover both March and December) is also not accurate. If you need to simulate a truly SSO orbit, you will need to consider station keeping, and this is a huge development. Unfortunately, we cannot provide the corresponding code. I know a program based on Orekit that does that perfectly throughout satellite lifetime (Skat, which I wrote together with Pascal), but it is proprietary program we did for Eumetsat so it belong to them.

One thing you may try is just putting a few almost arbitrary station-keeping dates every few months and just compensate semi-major axis and node drift to synchronize your orbit again with respect to SSO constraints, it will be more realistic than a single propagation.

Hello @luc !
My goal was to simulate the SSO orbit and propagate it forward up to a week for a simple analysis (eclipse detection, geographiczonedetection, etc). The initialization was to be either real TLE or a custom TLE (based on initial state/keplerian elements). I wanted to use the same propagator (SGP4) for both. But I guess I should drop the idea of creating a custom TLE based on single point (state/keplerians) and use a different propagator (EH, or numerical) for them.

You can look at the tutorials about TLE conversion with proper fitting, I think there is one.
You can also look at the conversion unit tests, I am sure there is one example there.

1 Like

thank you! I guess the problem with the eclipsedetector arises from the fact that my orbit is not really a SSO orbit (even though I get the starting inertial state vector from STK). When I plugin real TLE or proper orbital elements, I get better eclipsedetections.
off topic: is there a quick (simple) way to get an angle between two vectors that are defined using PVCoordinates in the same frame?

Hi,

The class Vector3D from Hipparchus has a static method called angle that does that.

Cheers,
Romain.

1 Like