Incorrect satellite rotation rate in nadir pointing mode

Orekit Python Wrapper usage
1

Hi all,

I currently am trying to simulate a switch in pointing laws in a satellite that is overflying a ground station. Once the satellite enters the ground station field of view, it should switch from nadir to target pointing, and once it exits, it should switch back to nadir pointing.

My goal is to extract the rotation rates of the satellite about all three body-fixed axes in the body-fixed frame.

I’ve defined pointing laws, event detectors and an attitude sequence that depends on the event detectors. Finally, after propagation, I added an angular rate getter, taken from here under “Direct use of transforms : coordinates with respect to spacecraft frame”.

I’ve come to the point where I can indeed obtain the angular rates and plot them. There is however a strange behaviour that you can see in this plot here:

Screenshot (1318)
Screenshot (1318)1610×537 41 KB

The rotation rate in nadir pointing mode does not match that of Earth, as it should be doing. In fact, it appears to be pretty much half of the real value (~0.06 deg/s).

The code is here (there is lots of visualisation and data saving information, but I hope the comments are clear enough so you can see where the actual programming of the attitudes sequence/pointing laws/etc. is taking place):

# %% [markdown]
# # Elevation, range & angular rate detector for Optical Ground Station overflight simulation

# %% [markdown]
# Initialise Orekit and bring up the Python-Java interface

# %%
import orekit
vm = orekit.initVM()

# %% [markdown]
# Now set up the pointer to the orekit-data.zip file, using one of the helper files. The file should be in current directory if not specified otherwise.

# %%
from orekit.pyhelpers import setup_orekit_curdir, absolutedate_to_datetime
setup_orekit_curdir()
#from orekit.pyhelpers import download_orekit_data_curdir
#download_orekit_data_curdir()

# %% [markdown]
# Now we are set up to import and use objects from the orekit library. Packages can be imported as they were native Python packages

# %%
from org.orekit.orbits import KeplerianOrbit, PositionAngle
from org.orekit.propagation import SpacecraftState;
from org.orekit.propagation.analytical import EcksteinHechlerPropagator
from org.orekit.propagation.sampling import OrekitFixedStepHandler, PythonOrekitFixedStepHandler
from org.orekit.time import AbsoluteDate, TimeScalesFactory
from org.orekit.utils import Constants, IERSConventions, AngularDerivativesFilter
from org.orekit.frames import FramesFactory, TopocentricFrame, LOFType
from org.orekit.bodies import OneAxisEllipsoid, CelestialBodyFactory, GeodeticPoint
from org.orekit.attitudes import AttitudesSequence
from org.orekit.attitudes import LofOffset
from org.orekit.attitudes import TargetPointing, NadirPointing 
from org.orekit.bodies import BodyShape

import matplotlib.pyplot as plt

from java.io import File

from math import radians, pi, degrees
import pandas as pd
import numpy as np
from scipy import stats

# %%
import math 
import pandas as pd

# %%
utc = TimeScalesFactory.getUTC()

# %% [markdown]
# Define orbit as per requirements

# %%
 ## Calculation of inclination according to Cap14

# constant of Sun-synchronicity kh (Cap14)

kh = 10.10949

satellite_height = 550e3       # [m]

semimajsatellite = Constants.WGS84_EARTH_EQUATORIAL_RADIUS + satellite_height

def incl_calc(semimaj):
    # reduced distance eta (semi-major axis over Earth radius)
    eta = semimaj/Constants.WGS84_EARTH_EQUATORIAL_RADIUS
    inclrad = np.arccos(-eta**(7/2)/kh)
    incl = np.rad2deg(inclrad)
    return incl

# %%
ra = satellite_height         #  Apogee
rp = satellite_height         #  Perigee
ideg = incl_calc(semimajsatellite)
i = radians(ideg)      # inclination
omega = radians(90.0)   # perigee argument, value selected to achieve "frozen" orbit as per Boa04
raan = radians(112.8)  # right ascension of ascending node
lv = radians(0.0)    # True anomaly

epochDate = AbsoluteDate(2024, 1, 1, 00, 00, 00.000, utc)
initial_date = epochDate

a = (rp + ra + 2 * Constants.WGS84_EARTH_EQUATORIAL_RADIUS) / 2.0    
e = 1.0 - (rp + Constants.WGS84_EARTH_EQUATORIAL_RADIUS) / a

## Inertial frame where the satellite is defined
inertialFrame = FramesFactory.getEME2000()

## Orbit construction as Keplerian
initialOrbit = KeplerianOrbit(a, e, i, omega, raan, lv,
                              PositionAngle.TRUE,
                              inertialFrame, epochDate, Constants.WGS84_EARTH_MU)

satellite_state = SpacecraftState(initialOrbit)

# %% [markdown]
# # Pointing laws definition
# Define pointing laws for satellite. In the context of this simulation, only two are required: a nadir pointing law, for when the satellite is outside the OGS FOV, and a target pointing law, for when the satellite is inside the OGS FOV.
# 
# Nadir pointing law defined with the NadirPointing class.

# %%
ITRF = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 
                         Constants.WGS84_EARTH_FLATTENING, 
                         ITRF)

# %%
nadirPointingLaw = NadirPointing(inertialFrame, earth)

# %% [markdown]
# Target pointing law defined in terms of inertial frame where satellite is defined, geodetic point where the target lies and the shape of the target. Define it using the TargetPointing class.

# %%
# according to .java Elevation Detector, ElevationDetector requires max. three input parameters: a maximum checking interval, a maximum divergence threshold and a topo reference to a topocentric model

# definition of topo reference to a topocentric model 
longitude = radians(13.071011)
latitude  = radians(53.329364)
altitude  = 76.0
station_point = GeodeticPoint(latitude, longitude, altitude)
station_frame = TopocentricFrame(earth, station_point, "Neustrelitz")

longitude_br = radians(203.743079)
latitude_br  = radians(20.706489)
altitude_br  = 3056.272
station_point_br = GeodeticPoint(latitude_br, longitude_br, altitude_br)
station_frame_br = TopocentricFrame(earth, station_point_br, "Hawaii")

# %%
targetPointingLaw = TargetPointing(inertialFrame, station_point, earth) 
targetPointingLaw_br = TargetPointing(inertialFrame, station_point_br, earth)

# %% [markdown]
# ## Adding Event Detectors

# %%
from org.orekit.propagation.events import EventsLogger, ElevationDetector
from org.orekit.propagation.events.handlers import ContinueOnEvent
from org.orekit.propagation.events.handlers import PythonEventHandler
from org.hipparchus.ode.events import Action

# %%
# define EventDetectors for the switching. Two specialized event detectors are built upon an EclipseDetector: one, invisibVisibEvent, to detect the invisible to visible transition, another, visibInvisibEvent, to detect the visible to invisible transition. This is done by overriding the eventOccurred method of the standard EclipseDetector. To override methods in a java class in Python, a special wrapped version of the class is used, called PythonElevationDetector. When subclassing Orekit classes in Python, all methods of that interface needs to be defined, if so just with a "pass" statement.

class myVisibDetector(PythonEventHandler):
    
    def init(self, initialstate, target, detector):
        pass
    
    def eventOccurred(self, s, detector, increasing):
        if increasing:
            print("\n", s.getDate()," : event occurred, entering OGS FOV\n")
        return Action.CONTINUE
    
    def resetState(self, detector, oldState):
        return oldState;

# %%
class myInvisibDetector(PythonEventHandler):
    
    def init(self, initialstate, target, detector):
        pass
    
    def eventOccurred(self, s, detector, increasing):
        if not increasing:
            print("\n", s.getDate()," : event occurred, exiting OGS FOV\n")
        return Action.CONTINUE
    
    def resetState(self, detector, oldState):
        return oldState

# %%
visibEvent = ElevationDetector(station_frame).withConstantElevation(radians(20)).withHandler(ContinueOnEvent())
visibEvent = visibEvent.withHandler(myVisibDetector().of_(ElevationDetector))

# %%
invisibEvent = ElevationDetector(station_frame).withConstantElevation(radians(20)).withHandler(ContinueOnEvent())
invisibEvent = invisibEvent.withHandler(myInvisibDetector().of_(ElevationDetector))

# %% [markdown]
# Visibility/Invisibility event for Hawaii OGS

# %%
visibEvent_br = ElevationDetector(station_frame_br).withConstantElevation(radians(20)).withHandler(ContinueOnEvent())
visibEvent_br = visibEvent_br.withHandler(myVisibDetector().of_(ElevationDetector))

# %%
invisibEvent_br = ElevationDetector(station_frame_br).withConstantElevation(radians(20)).withHandler(ContinueOnEvent())
invisibEvent_br = invisibEvent_br.withHandler(myInvisibDetector().of_(ElevationDetector))

# %% [markdown]
# # Attitude sequence definition
# Define sequence of attitude behaviours, i.e. two pointing laws, that get triggered by the occurrence of an Event:
# * the first one enables the transition from nadirPointingLaw to targetPointingLaw when an increasing visibEvent occurs
# * the second one enables the transition from targetPointingLaw to nadirPointingLaw  when a decreasing invisibEvent occurs
# 
# * see this [link](https://www.orekit.org/static/apidocs/org/orekit/attitudes/AttitudesSequence.html#addSwitchingCondition-org.orekit.attitudes.AttitudeProvider-org.orekit.attitudes.AttitudeProvider-T-boolean-boolean-double-org.orekit.utils.AngularDerivativesFilter-org.orekit.attitudes.AttitudesSequence.SwitchHandler-) under "addSwitchingCondition" for detailed explanation of parameters.
# * .USE_RRA to match rotation, rotation rate and rotation acceleration

# %%
attitudesSequence = AttitudesSequence()

attitudesSequence.addSwitchingCondition(nadirPointingLaw, targetPointingLaw,
                                        visibEvent, True, False, 10.0, 
                                        AngularDerivativesFilter.USE_RRA, None)

attitudesSequence.addSwitchingCondition(nadirPointingLaw, targetPointingLaw_br,
                                        visibEvent_br, True, False, 10.0, 
                                        AngularDerivativesFilter.USE_RRA, None)


attitudesSequence.addSwitchingCondition(targetPointingLaw, nadirPointingLaw, 
                                        invisibEvent, False, True, 10.0, 
                                        AngularDerivativesFilter.USE_RRA,  None)

attitudesSequence.addSwitchingCondition(targetPointingLaw_br, nadirPointingLaw, 
                                        invisibEvent_br, False, True, 10.0, 
                                        AngularDerivativesFilter.USE_RRA,  None)

# %%
if (visibEvent.g(SpacecraftState(initialOrbit)) >= 0):
    # initial position is outside of OGS FOV
    attitudesSequence.resetActiveProvider(nadirPointingLaw)
else:
    # initial position is in OGS FOV
    attitudesSequence.resetActiveProvider(targetPointingLaw)

# %% [markdown]
# # Orbit propagation

# %% [markdown]
# Use an EcksteinHechlerPropagator based on the analytical Eckstein-Hechler model. The propagator is built upon the initialOrbit, the attitudeSequence and physical constants for the gravity field potential.

# %%
propagator = EcksteinHechlerPropagator(initialOrbit, 
                                       attitudesSequence,
                                       Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 
                                       Constants.WGS84_EARTH_MU,Constants.WGS84_EARTH_C20, 
                                       0.0, 
                                       0.0, 
                                       0.0, 
                                       0.0)

# %% [markdown]
# Let attitudesSequence register the switching events to the propagator

# %%
attitudesSequence.registerSwitchEvents(propagator)

# %% [markdown]
# The propagation is executed in same way as in previous examples.

# %%
final_date = initial_date.shiftedBy(3600.0*13+30*60.0)
                                    #)*24.0*3.0
                                    #
# shift initial date to look at time window within same orbit propagated from initialOrbit with parameters given at initial date 2024-01-01 defined above
initial_date = initial_date.shiftedBy(3600.0*12)
state = propagator.propagate(initial_date, final_date)
state.getDate()

# %% [markdown]
# Add small loop that steps through the propagation and stores values in a Python array

# %%
pvs = []
ranges = []
ranges_br = []
extrapDate = initial_date

while (extrapDate.compareTo(final_date) <= 0.0):
    pv = propagator.getPVCoordinates(extrapDate, inertialFrame)
    range = station_frame.getRange(pv.getPosition(),
                                   inertialFrame,
                                   extrapDate)
    range_br = station_frame_br.getRange(pv.getPosition(),
                                   inertialFrame,
                                   extrapDate)
    ranges.append(range)
    ranges_br.append(range_br)
    pvs.append(pv)
    master_step = 5.0  # in seconds
    extrapDate = extrapDate.shiftedBy(master_step)

# %% [markdown]
# # Plotting results: Visibility times
# 
# Take the results and populate a Pandas DataFrame with them for managing the results.

# %%
prop_data = pd.DataFrame(data=pvs, columns=['pv'])

prop_data['Position'] = prop_data['pv'].apply(lambda x: x.getPosition())
prop_data['datetime'] = prop_data['pv'].apply(lambda x: absolutedate_to_datetime(x.getDate()))
prop_data['elevation'] = prop_data['pv'].apply(lambda x: station_frame.getElevation(x.getPosition(), inertialFrame, x.getDate())*180.0/pi )
prop_data['azimuth'] = prop_data['pv'].apply(lambda x: station_frame.getAzimuth(x.getPosition(), inertialFrame, x.getDate())*180.0/pi )
prop_data['elevation_br'] = prop_data['pv'].apply(lambda x: station_frame_br.getElevation(x.getPosition(), inertialFrame, x.getDate())*180.0/pi )
prop_data['azimuth_br'] = prop_data['pv'].apply(lambda x: station_frame_br.getAzimuth(x.getPosition(), inertialFrame, x.getDate())*180.0/pi )
prop_data.set_index('datetime', inplace=True, drop=False)
prop_data.index.name = 'Timestamp'
prop_data.head()

# %% [markdown]
# Calculate the satellite subpoint geographical coordinates and store them in the pandas DataFrame

# %%
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

# %% [markdown]
# Merge both data frames

# %%
# follow tips under https://stackoverflow.com/questions/12850345/how-do-i-combine-two-dataframes#12850453
# https://stackoverflow.com/a/40339932

range_data = pd.DataFrame(data=ranges, columns=['range'])
range_data_br = pd.DataFrame(data=ranges_br, columns=['range_br'])
prop_data.reset_index(drop=True, inplace=True)
range_data.reset_index(drop=True, inplace=True)
range_data_br.reset_index(drop=True, inplace=True)
prop_range_merged = pd.concat([prop_data, range_data, range_data_br], axis=1)
#prop_range_merged.set_index('datetime', inplace=True, drop=False)
#prop_range_merged.index.name = 'Timestamp'

# %% [markdown]
# Make separate vectors for combined contact constraint (range < 1500 km && elevation > 20 deg) for all ground stations

# %%
prop_range_merged['visible'] = prop_range_merged.elevation.apply(lambda el: 'Yes' if el>20 else 'No')
prop_range_merged['inrange'] = prop_range_merged.range.apply(lambda rng: 'Yes' if rng<1.5e6 else 'No')
prop_range_merged['visible_br'] = prop_range_merged.elevation_br.apply(lambda el: 'Yes' if el>20 else 'No')
prop_range_merged['inrange_br'] = prop_range_merged.range_br.apply(lambda rng: 'Yes' if rng<1.5e6 else 'No')
prop_range_merged.head()

# %% [markdown]
# Select contact data based on combined contact constraint (range < 1500 km and elevation > 20 deg) for all ground stations

# %%
# as per https://pandas.pydata.org/pandas-docs/stable/getting_started/intro_tutorials/03_subset_data.html#how-do-i-select-a-subset-of-a-dataframe
contact = prop_range_merged[((prop_range_merged["range"] < 150e6) & (prop_range_merged["elevation"] > 20)) | ((prop_range_merged["range_br"] < 150e6) & (prop_range_merged["elevation_br"] > 20))]
contact.head()

# %% [markdown]
# Plot range, elevation of entire overflight, and contact time, interactively with Plotly.

# %%
import plotly.express as px

# %%
px.line(prop_range_merged, y='elevation', x='datetime', hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation',
                    'latitude',
                    'longitude', 
                    'range'])

# %%
px.line(prop_range_merged, y='elevation_br', x='datetime', hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation_br',
                    'latitude',
                    'longitude', 
                    'range_br'])

# %%
px.line(contact, y='elevation', x='datetime', hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation',
                    'latitude',
                    'longitude', 
                    'range'])

# %%
px.line(prop_range_merged, y='range', x='datetime', hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation',
                    'latitude',
                    'longitude', 
                    'range'])

# %%
px.line(prop_range_merged, y='range_br', x='datetime', hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation_br',
                    'latitude',
                    'longitude', 
                    'range_br'])

# %%
px.line(contact, y='range', x='datetime', hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation',
                    'latitude',
                    'longitude', 
                    'range'])

# %% [markdown]
# # Angular rate extraction
# 
# The propagator operating mode is set to master mode with fixed step. The implementation of the interface OrekitFixedStepHandler aims to define the handleStep method called within the loop. 
# 
# The handleStep method will print at the current date the rotation rate of the spacecraft as a 3D vector, similar to [here](https://www.orekit.org/site-orekit-tutorials-11.3/tutorials/frames.html).

# %%
class mystephandler(PythonOrekitFixedStepHandler): 
    
    spins = []  # holds lists of Vector3D types (without conversion to array with toArray!)
    dates = []
    
    def init(self, s0, t, step):
        pass
        
    def handleStep(self, currentState):
        # SpacecraftState.toTransform(): Compute the transform from state defining frame to spacecraft frame as per https://www.orekit.org/static/apidocs/org/orekit/propagation/SpacecraftState.html#toTransform--
        inertToSpacecraft = currentState.toTransform()
        # getRotationRate(): "the result of the “getRotationRate()” call should return the angular velocity of refB with respect to refA expressed in the “refB”." as per https://forum.orekit.org/t/transform-between-two-reference-frames/1265. Assume here: spacecraft frame is B, inertial frame is A.
        spin = np.array(inertToSpacecraft.getRotationRate().toArray())
        self.spins.append(spin)
        self.dates.append(currentState.getDate())
        print("The current date is " + currentState.getDate().toString())

    def __str__(self):
        return str(self.spins)

    def finish(self, finalState):
        pass

# %%
handler = mystephandler()
propagator.cast_(propagator).getMultiplexer().add(master_step, handler)
#propagator.cast_(propagator).setStepHandler(180.0, handler)

# %%
#final_date = initial_date.shiftedBy(3600.0*11.5)
# shift initial date to look at time window within same orbit propagated from initialOrbit with parameters given at initial date 2024-01-01 defined above
#initial_date = initial_date.shiftedBy(3600.0*11.25)
finalState = propagator.propagate(initial_date, final_date)
#print(mystephandler.__str__(handler))


# %% [markdown]
# # Plotting results: angular rates

# %%
print("Propagation ended at " + finalState.getDate().toString())

# %%
from orekit.pyhelpers import absolutedate_to_datetime

# %%
pydates = [absolutedate_to_datetime(t) for t in handler.dates]

# %%
# store stepHandler angular rate values in list
spins = handler.spins

# from https://stackoverflow.com/a/903867
def column(matrix, i):
    return [row[i] for row in matrix]

# split up multidimensional list in three lists

spinsxrad = column(spins, 0)
spinsyrad = column(spins, 1)
spinszrad = column(spins, 2)

# convert rotation values to degrees

spinsx = np.degrees(spinsxrad)
spinsy = np.degrees(spinsyrad)
spinsz = np.degrees(spinszrad)

# %% [markdown]
# Make data frames for angular rate values

# %%
angular_rates_x = pd.DataFrame(data=spinsx, columns=['spinsx'])
angular_rates_y = pd.DataFrame(data=spinsy, columns=['spinsy'])
angular_rates_z = pd.DataFrame(data=spinsz, columns=['spinsz'])

# %% [markdown]
# Merge angular rate data frames with merged propagation and range data frame

# %%
prop_range_merged.reset_index(drop=True, inplace=True)
angular_rates_x.reset_index(drop=True, inplace=True)
angular_rates_y.reset_index(drop=True, inplace=True)
angular_rates_z.reset_index(drop=True, inplace=True)
complete_data = pd.concat([prop_range_merged, angular_rates_x, angular_rates_y, angular_rates_z], axis=1)
complete_data.set_index('datetime', inplace=True, drop=False)
complete_data.index.name = 'Timestamp'
complete_data.head()

# %% [markdown]
# Apply combined contact constraint on complete data DataFrame.

# %%
contact_rates = complete_data[((complete_data["range"] < 150e6) & (complete_data["elevation"] > 20)) | ((complete_data["range_br"] < 150e6) & (complete_data["elevation_br"] > 20))]
contact_rates.head()

# %% [markdown]
# Display angular rates, both for entire overflight and during contact window, visually with Plotly

# %%
figspinsxcomplete = px.line(complete_data, 
        y='spinsx', 
        x='datetime', 
        hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation',
                    'latitude',
                    'longitude', 
                    'range', 
                    'spinsx'], 
        labels=dict(spinsx="Angular rate about body-frame x-axis, deg/s", datetime="Time"))
figspinsxcomplete.show()

# %%
px.line(contact_rates, 
        y='spinsx', 
        x='datetime', 
        hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation',
                    'latitude',
                    'longitude', 
                    'range', 
                    'spinsx'],
        labels=dict(spinsx="Angular rate about body-frame x-axis, deg/s", datetime="Time"))

# %%
figspinsycomplete = px.line(complete_data, 
        y='spinsy', 
        x='datetime', 
        hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation',
                    'latitude',
                    'longitude', 
                    'range', 
                    'spinsy'], 
        labels=dict(spinsy="Angular rate about body-frame y-axis, deg/s", datetime="Time"))
figspinsycomplete.show()

# %%
px.line(contact_rates, 
        y='spinsy', 
        x='datetime', 
        hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation',
                    'latitude',
                    'longitude', 
                    'range', 
                    'spinsy'],
        labels=dict(spinsy="Angular rate about body-frame y-axis, deg/s", datetime="Time"))

# %%
figspinszcomplete = px.line(complete_data, 
        y='spinsz', 
        x='datetime', 
        hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation',
                    'latitude',
                    'longitude', 
                    'range', 
                    'spinsz'], 
        labels=dict(spinsz="Angular rate about body-frame z-axis, deg/s", datetime="Time"))
figspinszcomplete.show()

# %%
px.line(contact_rates, 
        y='spinsz', 
        x='datetime', 
        hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation',
                    'latitude',
                    'longitude', 
                    'range', 
                    'spinsz'],
        labels=dict(spinsz="Angular rate about body-frame z-axis, deg/s", datetime="Time"))

# %% [markdown]
# Calculate absolute value of angular rate in body-fixed frame

# %%
angular_rate_abs = np.sqrt(np.square(complete_data["spinsx"]) + np.square(complete_data["spinsy"]) + np.square(complete_data["spinsz"]))
angular_rate_abs_df = pd.DataFrame(data=angular_rate_abs, columns=['angular_rate_abs'])

# %% [markdown]
# Pot all three angular rate graphs, along with absolute value of angular rate, in a single plot

# %%
complete_data.reset_index(drop=True, inplace=True)
angular_rate_abs_df.reset_index(drop=True, inplace=True)
angular_rates_complete = pd.concat([complete_data, angular_rate_abs_df], axis=1)
angular_rates_complete.set_index('datetime', inplace=True, drop=False)
angular_rates_complete.index.name = 'Timestamp'
all_rates_fig = px.line(angular_rates_complete, x="datetime", y=["spinsx", "spinsy", "spinsz", "angular_rate_abs"],)
all_rates_fig.show()

Has anybody encountered a similar problem before? I am not sure what is occurring there and I’d be very glad to have your input on this.

Cheers,

Alex

2

Looking at the range you get, I wonder if you are observing the station from the opposite side of the Earth, through it.

3

Hi @der-a.bfuehler,

Welcome to the forum !

I must say I don’t get it neither. I’ve run your code and found the same issue.

At initial date, both attitudes and spins should be the same but it’s not the case

attitudesSequence.getAttitude(propagator, initial_date, inertialFrame).getSpin().scalarMultiply(180/np.pi)
<Vector3D: {0,001257726; -0,0326834274; -0,0038981296}>
nadirPointingLaw.getAttitude(propagator, initial_date, inertialFrame).getSpin().scalarMultiply(180/np.pi)
<Vector3D: {0,0000007366; -0,0619089561; -0,0000268643}>

The attitude sequence is definitely the cause of the problem here but I don’t know why.
I will have to translate your code to Java and debug more thoroughly to understand what’s happening.
If anyone has a better idea I’ll be glad not to go through this :sweat_smile:

Cheers,
Maxime

4

Thank you both for your replies.
@Maxime, I think your reply put me on the right track by getting the rotation rates at the outset of the propagation.

Beforehand, the problem with the rotation rate being ~0.03 deg/s while not in FOV of the station was due to the fact that it was in target pointing mode the entire time. I cross-checked the rotation rate in target pointing mode at propagation outset, and indeed, ~0.03 deg/s is the value if the propagation starts in target pointing mode.

I reworked the EventDetectors and the attitudesSequence according to this example (I wrote down the correspondences that I made just above the attitudesSequence definition) and there is a rotation rate of ~0.06 deg/s now, only that it lasts throughout the entire propagation. So now, the propagation starts with the correct rotation rate value, however the problem of the EventDetectors appearing to register switching (according to the print statements), yet the issue of the actual switching between the pointing laws not occurring still remains.

newplot
newplot1157×450 27.8 KB

I am still not sure what is happening here and would be glad to hear your suggestions :slight_smile:

Here is the updated code:

# %% [markdown]
# # Elevation, range & angular rate detector for Optical Ground Station overflight simulation

# %% [markdown]
# Initialise Orekit and bring up the Python-Java interface

# %%
import orekit
vm = orekit.initVM()

# %% [markdown]
# Now set up the pointer to the orekit-data.zip file, using one of the helper files. The file should be in current directory if not specified otherwise.

# %%
from orekit.pyhelpers import setup_orekit_curdir, absolutedate_to_datetime
setup_orekit_curdir()
#from orekit.pyhelpers import download_orekit_data_curdir
#download_orekit_data_curdir()

# %% [markdown]
# Now we are set up to import and use objects from the orekit library. Packages can be imported as they were native Python packages

# %%
from org.orekit.orbits import KeplerianOrbit, PositionAngle
from org.orekit.propagation import SpacecraftState;
from org.orekit.propagation.analytical import EcksteinHechlerPropagator
from org.orekit.propagation.sampling import OrekitFixedStepHandler, PythonOrekitFixedStepHandler
from org.orekit.time import AbsoluteDate, TimeScalesFactory
from org.orekit.utils import Constants, IERSConventions, AngularDerivativesFilter
from org.orekit.frames import FramesFactory, TopocentricFrame, LOFType
from org.orekit.bodies import OneAxisEllipsoid, CelestialBodyFactory, GeodeticPoint
from org.orekit.attitudes import AttitudesSequence
from org.orekit.attitudes import LofOffset
from org.orekit.attitudes import TargetPointing, NadirPointing 
from org.orekit.bodies import BodyShape

import matplotlib.pyplot as plt

from java.io import File

from math import radians, pi, degrees
import pandas as pd
import numpy as np
from scipy import stats

# %%
import math 
import pandas as pd

# %%
utc = TimeScalesFactory.getUTC()

# %% [markdown]
# Define orbit as per requirements

# %%
 ## Calculation of inclination according to Cap14

# constant of Sun-synchronicity kh (Cap14)

kh = 10.10949

satellite_height = 550e3       # [m]

semimajsatellite = Constants.WGS84_EARTH_EQUATORIAL_RADIUS + satellite_height

def incl_calc(semimaj):
    # reduced distance eta (semi-major axis over Earth radius)
    eta = semimaj/Constants.WGS84_EARTH_EQUATORIAL_RADIUS
    inclrad = np.arccos(-eta**(7/2)/kh)
    incl = np.rad2deg(inclrad)
    return incl

# %%
ra = satellite_height         #  Apogee
rp = satellite_height         #  Perigee
ideg = incl_calc(semimajsatellite)
i = radians(ideg)      # inclination
omega = radians(90.0)   # perigee argument, value selected to achieve "frozen" orbit as per Boa04
raan = radians(112.8)  # right ascension of ascending node
lv = radians(0.0)    # True anomaly

epochDate = AbsoluteDate(2024, 1, 1, 00, 00, 00.000, utc)
initial_date = epochDate

a = (rp + ra + 2 * Constants.WGS84_EARTH_EQUATORIAL_RADIUS) / 2.0    
e = 1.0 - (rp + Constants.WGS84_EARTH_EQUATORIAL_RADIUS) / a

## Inertial frame where the satellite is defined
inertialFrame = FramesFactory.getEME2000()

## Orbit construction as Keplerian
initialOrbit = KeplerianOrbit(a, e, i, omega, raan, lv,
                              PositionAngle.TRUE,
                              inertialFrame, epochDate, Constants.WGS84_EARTH_MU)

satellite_state = SpacecraftState(initialOrbit)

# %% [markdown]
# # Pointing laws definition
# Define pointing laws for satellite. In the context of this simulation, only two are required: a nadir pointing law, for when the satellite is outside the OGS FOV, and a target pointing law, for when the satellite is inside the OGS FOV.
# 
# Nadir pointing law defined with the NadirPointing class.

# %%
ITRF = FramesFactory.getITRF(IERSConventions.IERS_2010, True)
earth = OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 
                         Constants.WGS84_EARTH_FLATTENING, 
                         ITRF)

# %%
nadirPointingLaw = NadirPointing(inertialFrame, earth)

# %% [markdown]
# Target pointing law defined in terms of inertial frame where satellite is defined, geodetic point where the target lies and the shape of the target. Define it using the TargetPointing class.

# %%
# according to .java Elevation Detector, ElevationDetector requires max. three input parameters: a maximum checking interval, a maximum divergence threshold and a topo reference to a topocentric model

# definition of topo reference to a topocentric model 
longitude = radians(13.071011)
latitude  = radians(53.329364)
altitude  = 76.0
station_point = GeodeticPoint(latitude, longitude, altitude)
station_frame = TopocentricFrame(earth, station_point, "Neustrelitz")

longitude_br = radians(203.743079)
latitude_br  = radians(20.706489)
altitude_br  = 3056.272
station_point_br = GeodeticPoint(latitude_br, longitude_br, altitude_br)
station_frame_br = TopocentricFrame(earth, station_point_br, "Hawaii")

# %%
targetPointingLaw = TargetPointing(inertialFrame, station_point, earth) 

# %% [markdown]
# ## Adding Event Detectors

# %%
from org.orekit.propagation.events import EventsLogger, ElevationDetector
from org.orekit.propagation.events.handlers import ContinueOnEvent
from org.orekit.propagation.events.handlers import PythonEventHandler
from org.hipparchus.ode.events import Action

# %%
# define EventDetectors for the switching. Two specialized event detectors are built upon an EclipseDetector: one, invisToVisEvent, to detect the invisible to visible transition, another, visToInvisEvent, to detect the visible to invisible transition. This is done by overriding the eventOccurred method of the standard EclipseDetector. To override methods in a java class in Python, a special wrapped version of the class is used, called PythonElevationDetector. When subclassing Orekit classes in Python, all methods of that interface needs to be defined, if so just with a "pass" statement.

class visToInvisDetector(PythonEventHandler):
    
    def init(self, initialstate, target, detector):
        pass
    
    def eventOccurred(self, s, detector, increasing):
        if not increasing:
            print("\n", s.getDate()," : event occurred, exiting OGS FOV, switching to nadir pointing law\n")
        return Action.CONTINUE
    
    def resetState(self, detector, oldState):
        return oldState;

# %%
class invisToVisDetector(PythonEventHandler):
    
    def init(self, initialstate, target, detector):
        pass
    
    def eventOccurred(self, s, detector, increasing):
        if increasing:
            print("\n", s.getDate()," : event occurred, entering OGS FOV, switching to target pointing law\n")
        return Action.CONTINUE
    
    def resetState(self, detector, oldState):
        return oldState

# %%
visToInvisEvent = ElevationDetector(station_frame).withConstantElevation(radians(20)).withHandler(ContinueOnEvent())
visToInvisEvent = visToInvisEvent.withHandler(visToInvisDetector().of_(ElevationDetector))

# %%
invisToVisEvent = ElevationDetector(station_frame).withConstantElevation(radians(20)).withHandler(ContinueOnEvent())
invisToVisEvent = invisToVisEvent.withHandler(invisToVisDetector().of_(ElevationDetector))

# %% [markdown]
# # Attitude sequence definition
# Define sequence of attitude behaviours, i.e. two pointing laws, that get triggered by the occurrence of an Event:
# * the first one enables the transition from targetPointingLaw to nadirPointingLaw when a decreasing visToInvisEvent occurs
# * the second one enables the transition from nadirPointingLaw to targetPointingLaw when an increasing invisToVisEvent occurs
# 
# * see this [link](https://www.orekit.org/static/apidocs/org/orekit/attitudes/AttitudesSequence.html#addSwitchingCondition-org.orekit.attitudes.AttitudeProvider-org.orekit.attitudes.AttitudeProvider-T-boolean-boolean-double-org.orekit.utils.AngularDerivativesFilter-org.orekit.attitudes.AttitudesSequence.SwitchHandler-) under "addSwitchingCondition" for detailed explanation of parameters.
# * .USE_RRA to match rotation, rotation rate and rotation acceleration

# %%
# invis == night == nadir, vis == day == target
# decreasing day to night == decreasing visToInvisEvent
# increasing night to day == increasing invisToVisEvent

attitudesSequence = AttitudesSequence()

# boolean parameters: switchOnIncrease, switchOnDecrease, in that order

attitudesSequence.addSwitchingCondition(targetPointingLaw, nadirPointingLaw,
                                        visToInvisEvent, False, True, 10.0, 
                                        AngularDerivativesFilter.USE_RRA, None)


attitudesSequence.addSwitchingCondition(nadirPointingLaw, targetPointingLaw, 
                                        invisToVisEvent, True, False, 10.0, 
                                        AngularDerivativesFilter.USE_RRA,  None)

# %%
if (visToInvisEvent.g(SpacecraftState(initialOrbit)) >= 0):
    # initial position is inside of OGS FOV
    attitudesSequence.resetActiveProvider(targetPointingLaw);
else:
    # initial position is outside of OGS FOV
    attitudesSequence.resetActiveProvider(nadirPointingLaw);

# %% [markdown]
# # Orbit propagation

# %% [markdown]
# Use an EcksteinHechlerPropagator based on the analytical Eckstein-Hechler model. The propagator is built upon the initialOrbit, the attitudeSequence and physical constants for the gravity field potential.

# %%
propagator = EcksteinHechlerPropagator(initialOrbit, 
                                       attitudesSequence,
                                       Constants.WGS84_EARTH_EQUATORIAL_RADIUS, 
                                       Constants.WGS84_EARTH_MU,Constants.WGS84_EARTH_C20, 
                                       0.0, 
                                       0.0, 
                                       0.0, 
                                       0.0)

# %% [markdown]
# Let attitudesSequence register the switching events to the propagator

# %%
attitudesSequence.registerSwitchEvents(propagator)

# %% [markdown]
# The propagation is executed in same way as in previous examples.

# %%
final_date = initial_date.shiftedBy(3600.0*13+30*60.0)
                                    #)*24.0*3.0
                                    #
# shift initial date to look at time window within same orbit propagated from initialOrbit with parameters given at initial date 2024-01-01 defined above
initial_date = initial_date.shiftedBy(3600.0*10+30*60.0)
state = propagator.propagate(initial_date, final_date)
state.getDate()

# %% [markdown]
# Add small loop that steps through the propagation and stores values in a Python array

# %%
pvs = []
ranges = []
ranges_br = []
extrapDate = initial_date

while (extrapDate.compareTo(final_date) <= 0.0):
    pv = propagator.getPVCoordinates(extrapDate, inertialFrame)
    range = station_frame.getRange(pv.getPosition(),
                                   inertialFrame,
                                   extrapDate)
    range_br = station_frame_br.getRange(pv.getPosition(),
                                   inertialFrame,
                                   extrapDate)
    ranges.append(range)
    ranges_br.append(range_br)
    pvs.append(pv)
    master_step = 5.0  # in seconds
    extrapDate = extrapDate.shiftedBy(master_step)

# %% [markdown]
# # Plotting results: Visibility times
# 
# Take the results and populate a Pandas DataFrame with them for managing the results.

# %%
prop_data = pd.DataFrame(data=pvs, columns=['pv'])

prop_data['Position'] = prop_data['pv'].apply(lambda x: x.getPosition())
prop_data['datetime'] = prop_data['pv'].apply(lambda x: absolutedate_to_datetime(x.getDate()))
prop_data['elevation'] = prop_data['pv'].apply(lambda x: station_frame.getElevation(x.getPosition(), inertialFrame, x.getDate())*180.0/pi )
prop_data['azimuth'] = prop_data['pv'].apply(lambda x: station_frame.getAzimuth(x.getPosition(), inertialFrame, x.getDate())*180.0/pi )
prop_data['elevation_br'] = prop_data['pv'].apply(lambda x: station_frame_br.getElevation(x.getPosition(), inertialFrame, x.getDate())*180.0/pi )
prop_data['azimuth_br'] = prop_data['pv'].apply(lambda x: station_frame_br.getAzimuth(x.getPosition(), inertialFrame, x.getDate())*180.0/pi )
prop_data.set_index('datetime', inplace=True, drop=False)
prop_data.index.name = 'Timestamp'
prop_data.head()

# %% [markdown]
# Calculate the satellite subpoint geographical coordinates and store them in the pandas DataFrame

# %%
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

# %% [markdown]
# Merge both data frames

# %%
# follow tips under https://stackoverflow.com/questions/12850345/how-do-i-combine-two-dataframes#12850453
# https://stackoverflow.com/a/40339932

range_data = pd.DataFrame(data=ranges, columns=['range'])
range_data_br = pd.DataFrame(data=ranges_br, columns=['range_br'])
prop_data.reset_index(drop=True, inplace=True)
range_data.reset_index(drop=True, inplace=True)
range_data_br.reset_index(drop=True, inplace=True)
prop_range_merged = pd.concat([prop_data, range_data, range_data_br], axis=1)
#prop_range_merged.set_index('datetime', inplace=True, drop=False)
#prop_range_merged.index.name = 'Timestamp'

# %% [markdown]
# Make separate vectors for combined contact constraint (range < 1500 km && elevation > 20 deg) for all ground stations

# %%
prop_range_merged['visible'] = prop_range_merged.elevation.apply(lambda el: 'Yes' if el>20 else 'No')
prop_range_merged['inrange'] = prop_range_merged.range.apply(lambda rng: 'Yes' if rng<1.5e6 else 'No')
prop_range_merged['visible_br'] = prop_range_merged.elevation_br.apply(lambda el: 'Yes' if el>20 else 'No')
prop_range_merged['inrange_br'] = prop_range_merged.range_br.apply(lambda rng: 'Yes' if rng<1.5e6 else 'No')
prop_range_merged.head()

# %% [markdown]
# Select contact data based on combined contact constraint (range < 1500 km and elevation > 20 deg) for all ground stations

# %%
# as per https://pandas.pydata.org/pandas-docs/stable/getting_started/intro_tutorials/03_subset_data.html#how-do-i-select-a-subset-of-a-dataframe
contact = prop_range_merged[((prop_range_merged["range"] < 150e6) & (prop_range_merged["elevation"] > 20)) | ((prop_range_merged["range_br"] < 150e6) & (prop_range_merged["elevation_br"] > 20))]
contact.head()

# %% [markdown]
# Plot range, elevation of entire overflight, and contact time, interactively with Plotly.

# %%
import plotly.express as px

# %%
px.line(prop_range_merged, y='elevation', x='datetime', hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation',
                    'latitude',
                    'longitude', 
                    'range'])

# %%
px.line(prop_range_merged, y='elevation_br', x='datetime', hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation_br',
                    'latitude',
                    'longitude', 
                    'range_br'])

# %%
px.line(contact, y='elevation', x='datetime', hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation',
                    'latitude',
                    'longitude', 
                    'range'])

# %%
px.line(prop_range_merged, y='range', x='datetime', hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation',
                    'latitude',
                    'longitude', 
                    'range'])

# %%
px.line(prop_range_merged, y='range_br', x='datetime', hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation_br',
                    'latitude',
                    'longitude', 
                    'range_br'])

# %%
px.line(contact, y='range', x='datetime', hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation',
                    'latitude',
                    'longitude', 
                    'range'])

# %% [markdown]
# # Angular rate extraction
# 
# The propagator operating mode is set to master mode with fixed step. The implementation of the interface OrekitFixedStepHandler aims to define the handleStep method called within the loop. 
# 
# The handleStep method will print at the current date the rotation rate of the spacecraft as a 3D vector, similar to [here](https://www.orekit.org/site-orekit-tutorials-11.3/tutorials/frames.html).

# %%
class mystephandler(PythonOrekitFixedStepHandler): 
    
    spins = []  # holds lists of Vector3D types (without conversion to array with toArray!)
    dates = []
    
    def init(self, s0, t, step):
        pass
        
    def handleStep(self, currentState):
        # SpacecraftState.toTransform(): Compute the transform from state defining frame to spacecraft frame as per https://www.orekit.org/static/apidocs/org/orekit/propagation/SpacecraftState.html#toTransform--
        inertToSpacecraft = currentState.toTransform()
        # getRotationRate(): "the result of the “getRotationRate()” call should return the angular velocity of refB with respect to refA expressed in the “refB”." as per https://forum.orekit.org/t/transform-between-two-reference-frames/1265. Assume here: spacecraft frame is B, inertial frame is A.
        spin = np.array(inertToSpacecraft.getRotationRate().toArray())
        self.spins.append(spin)
        self.dates.append(currentState.getDate())
        print("The current date is " + currentState.getDate().toString())

    def __str__(self):
        return str(self.spins)

    def finish(self, finalState):
        pass

# %%
handler = mystephandler()
propagator.cast_(propagator).getMultiplexer().add(master_step, handler)
#propagator.cast_(propagator).setStepHandler(180.0, handler)

# %%
#final_date = initial_date.shiftedBy(3600.0*11.5)
# shift initial date to look at time window within same orbit propagated from initialOrbit with parameters given at initial date 2024-01-01 defined above
#initial_date = initial_date.shiftedBy(3600.0*11.25)
finalState = propagator.propagate(initial_date, final_date)
#print(mystephandler.__str__(handler))


# %% [markdown]
# # Plotting results: angular rates

# %%
print("Propagation ended at " + finalState.getDate().toString())

# %%
from orekit.pyhelpers import absolutedate_to_datetime

# %%
pydates = [absolutedate_to_datetime(t) for t in handler.dates]

# %%
# store stepHandler angular rate values in list
spins = handler.spins

# from https://stackoverflow.com/a/903867
def column(matrix, i):
    return [row[i] for row in matrix]

# split up multidimensional list in three lists

spinsxrad = column(spins, 0)
spinsyrad = column(spins, 1)
spinszrad = column(spins, 2)

# convert rotation values to degrees

spinsx = np.degrees(spinsxrad)
spinsy = np.degrees(spinsyrad)
spinsz = np.degrees(spinszrad)

# %% [markdown]
# Make data frames for angular rate values

# %%
angular_rates_x = pd.DataFrame(data=spinsx, columns=['spinsx'])
angular_rates_y = pd.DataFrame(data=spinsy, columns=['spinsy'])
angular_rates_z = pd.DataFrame(data=spinsz, columns=['spinsz'])

# %% [markdown]
# Merge angular rate data frames with merged propagation and range data frame

# %%
prop_range_merged.reset_index(drop=True, inplace=True)
angular_rates_x.reset_index(drop=True, inplace=True)
angular_rates_y.reset_index(drop=True, inplace=True)
angular_rates_z.reset_index(drop=True, inplace=True)
complete_data = pd.concat([prop_range_merged, angular_rates_x, angular_rates_y, angular_rates_z], axis=1)
complete_data.set_index('datetime', inplace=True, drop=False)
complete_data.index.name = 'Timestamp'
complete_data.head()

# %% [markdown]
# Apply combined contact constraint on complete data DataFrame.

# %%
contact_rates = complete_data[((complete_data["range"] < 150e6) & (complete_data["elevation"] > 20)) | ((complete_data["range_br"] < 150e6) & (complete_data["elevation_br"] > 20))]
contact_rates.head()

# %% [markdown]
# Get angular rates both from `attitudesSequence.getAttitude` and `nadirPointingLaw.getAttitude` as per [this](https://forum.orekit.org/t/incorrect-satellite-rotation-rate-in-nadir-pointing-mode/2554/3) recommendation in the Orekit help forum.

# %%
# from attitudesSequence.getAttitude
attitudesSequence.getAttitude(propagator, initial_date, inertialFrame).getSpin().scalarMultiply(180/np.pi)

# %%
# from nadirPointingLaw.getAttitude
nadirPointingLaw.getAttitude(propagator, initial_date, inertialFrame).getSpin().scalarMultiply(180/np.pi)

# %% [markdown]
# Get initial angular rate in target pointing for elimination cross-checking

# %%
targetPointingLaw.getAttitude(propagator, initial_date, inertialFrame).getSpin().scalarMultiply(180/np.pi)

# %% [markdown]
# Display angular rates, both for entire overflight and during contact window, visually with Plotly

# %%
figspinsxcomplete = px.line(complete_data, 
        y='spinsx', 
        x='datetime', 
        hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation',
                    'latitude',
                    'longitude', 
                    'range', 
                    'spinsx'], 
        labels=dict(spinsx="Angular rate about body-frame x-axis, deg/s", datetime="Time"))
figspinsxcomplete.show()

# %%
px.line(contact_rates, 
        y='spinsx', 
        x='datetime', 
        hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation',
                    'latitude',
                    'longitude', 
                    'range', 
                    'spinsx'],
        labels=dict(spinsx="Angular rate about body-frame x-axis, deg/s", datetime="Time"))

# %%
figspinsycomplete = px.line(complete_data, 
        y='spinsy', 
        x='datetime', 
        hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation',
                    'latitude',
                    'longitude', 
                    'range', 
                    'spinsy', 'visible'], 
        labels=dict(spinsy="Angular rate about body-frame y-axis, deg/s", datetime="Time"))
figspinsycomplete.show()

# %%
px.line(contact_rates, 
        y='spinsy', 
        x='datetime', 
        hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation',
                    'latitude',
                    'longitude', 
                    'range', 
                    'spinsy'],
        labels=dict(spinsy="Angular rate about body-frame y-axis, deg/s", datetime="Time"))

# %%
figspinszcomplete = px.line(complete_data, 
        y='spinsz', 
        x='datetime', 
        hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation',
                    'latitude',
                    'longitude', 
                    'range', 
                    'spinsz'], 
        labels=dict(spinsz="Angular rate about body-frame z-axis, deg/s", datetime="Time"))
figspinszcomplete.show()

# %%
px.line(contact_rates, 
        y='spinsz', 
        x='datetime', 
        hover_name="datetime", 
        hover_data=['azimuth',
                    'elevation',
                    'latitude',
                    'longitude', 
                    'range', 
                    'spinsz'],
        labels=dict(spinsz="Angular rate about body-frame z-axis, deg/s", datetime="Time"))

# %% [markdown]
# Calculate absolute value of angular rate in body-fixed frame

# %%
angular_rate_abs = np.sqrt(np.square(complete_data["spinsx"]) + np.square(complete_data["spinsy"]) + np.square(complete_data["spinsz"]))
angular_rate_abs_df = pd.DataFrame(data=angular_rate_abs, columns=['angular_rate_abs'])

# %% [markdown]
# Pot all three angular rate graphs, along with absolute value of angular rate, in a single plot

# %%
complete_data.reset_index(drop=True, inplace=True)
angular_rate_abs_df.reset_index(drop=True, inplace=True)
angular_rates_complete = pd.concat([complete_data, angular_rate_abs_df], axis=1)
angular_rates_complete.set_index('datetime', inplace=True, drop=False)
angular_rates_complete.index.name = 'Timestamp'
all_rates_fig = px.line(angular_rates_complete, x="datetime", y=["spinsx", "spinsy", "spinsz", "angular_rate_abs"],)
all_rates_fig.show()
5

:smiling_face_with_tear: ah ah, so the problem was that we were never in Nadir pointing… sorry I missed that.

Looks like the attitude switch never occurs… I don’t know why. I’ll try to take some time to look into it.
Just some side remarks:

  • I don’t think you need to create two identical detectors for the switches
  • Maybe you could try to implement the SwitchHandler in AttitudeSequence and print the event occurrence there so we are sure that the switch is actually triggered.
  • I saw that you perform the propagation twice. I think you could do it just once. Know that you can have several StepHandler with different sampling steps added to a single propagator.

Maxime

6

Hi @der-a.bfuehler,

I managed to make the switch work by forcing the maxCheckInterval of the event detector.
It’s 600s by default so it seems the propagator doesn’t “see” the event in some cases.
Setting it to 60s (or event 300s) fixed the problem.

This means changing this:

visToInvisEvent = ElevationDetector(station_frame).withConstantElevation(radians(20)).withHandler(ContinueOnEvent())

By adding a max check:

visToInvisEvent = ElevationDetector(station_frame).withConstantElevation(radians(20)).withHandler(ContinueOnEvent()).withMaxCheck(60.)

Now, the switch occurs but 10s of transition is very small and it looks like the satellite is rotating in the wrong direction (doing a full round-house turn instead of going for the shortest path), see the figure below.
I’m not sure, you should probably print the Euler angles of the attitude to study this further.

image
image1369×520 10 KB

Also, I tried to implement the PythonSwitchHandler class but I couldn’t make it work (I haven’t battled a lot though).

Anyway, good luck with your investigations, and don’t hesitate to ask if you have further questions.

Maxime

7

Hi Maxime,

Thanks a lot for your answers, I had encountered the same behaviour as you had in this image by modifying the maxCheckInterval as well, unfortunate those peaks in angular rate indeed. Thank you for your effort!

About point 3 you had mentioned in your first message: I had thought so too, that performing the propagation twice could maybe play a role in this problem… I will try to modify that so that I can have the propagation just once and see what that does.