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()