1. Orbit Propagator

This blog post introduces a high precision 6D orbit propagator that I developed for my research. This 6D orbit propagator leverages Orekit low-level spaceflight dynamics library (3D dynamics only) and CNES VTS (Celestia) visualization tool. Below is the structure of the proposed 6D flight dynamics library (Figure 1).

Figure 1: Layout of the proposed high precision 6D orbit propagator.

There are 6 classes that handles the associated tasks. State class contains all functions and variables that defines the physical shape, size, and material properties of space objects, celestial bodies, frame definitions and conversions. Force class includes functions that compute the conservative, such as non-aspherical gravitational potential and third-body mass, and non-conservative forces, such as solar radiation pressure and air drag. Torque class contains functions for computing gravity gradient torque for gravitational forces and solar radiation pressure and air drag torque for non-gravitational forces. Propagator class includes functions to propagate the orbital equations of motion in Newtonian formalism using adaptive numerical integrators such as Runge-Kutta integrator of order 8 or Dormand and Prince of order 8. Visibility class contains the functions for various methods to compute the bidirectional reflectance distribution function (BRDF) of the space objects based on visibility geometry. All reference frames use the IERS 2010 conventions, and earth orientation parameters (EOP) are updated regularly.

2. Demo Setup

For the demo, two extreme cases in terms of sensitivity of the space objects, namely massive cannonball and a small MLI-type material that is wrinkled, to the non-gravitational forces are presented in the demo. Cannonball is less sensitive to the non-gravitational forces, whereas non-gravitational forces and torques drastically affect the attitude/orbit evolution of high area-to-mass ratio wrinkled sheet, which is composed of 10 flat surfaces. Figure 2 shows the shapes of the space objects that are selected for the demo. Table 1 shows some of the physical parameters used for the demo. In Table 1, $d$ and $s$ are fractions of the diffuse and specular reflection ($s+d=1$), $\rho$ and $F_{0}$ are diffuse reflectance and specular reflectance of the surface at normal incidence, and $n_{u}$ and $n_{v}$ are phong-like exponents that control the shape of the specular lobe.

Figure 2 : Shapes of cannonball and HAMR object that are used for the demo.

Properties	Cannonball	HAMR object
Mass	7850 kg	0.0001125 kg
Area	1.209 $m^2$	0.0025 $m^2$ (each facet)
Sunlight intensity	455 $W/m^2$	455 $W/m^2$
d	1.0	0.5
s	0.0	0.5
$\rho$	0.26	0.26
$F_{0}$	-	0.5
$n_{u}$	-	0.5
$n_{v}$	-	0.5

Table 1 : Physical model parameters used for the space objects.

For the initial orbit, the orbital elemets of a satellite (OPS 5111) are used by transforming them from TEME to J2000 reference frames. First, a two-line element set is fetched from spacetrack.org, and then the state (cartesian coordinates) at epoch is converted to J2000. The TLE used is ,

0 OPS 5111

1 10684U 78020A 20110.24021950 +.00000083 +00000-0 +00000-0 0 9994

2 10684 063.0042 120.7837 0079708 206.3069 153.3306 01.98073183291879

The air drag coefficient is $2.8$ (for both objects) and aspherical gravitational potential order and degree are 20. The initial rotation rate for the HAMR object is $0.5$ degrees/s around the orbital normal axis.

3. Results

Figure 3 and Figure 4 show the light intensity observed by the ground telescope for the hamr object and the cannonball respectively. Figure 5 shows the relative orbital evolution starting from the epoch (1 day propagation starting at 2020-04-19T05:45:54.964 (UTC)). Please note that the time window for the visibility starts (for the location at lat: $57.0^{\circ}$ and long: $-5.0^{\circ}$) at 2020-04-19T21:17:04.964 (UTC) and ends at 2020-04-19T22:55:24.964 (UTC).

Figure 3 : Light intensity observed by the ground telescope for hamr object.

Figure 4 : Light intensity observed by the ground telescope for cannonball.

Figure 5 : Relative orbital evolution of space objects.

The 3D visualisation for the relative motion (orbits diverge because hamr object is drastically perturbed by the non-gravitational forces) can be reached from the link,

https://drive.google.com/file/d/1FeZUaFtIT5W8D5UL9v2Mr_PMCZL0Nqz8/view?usp=sharing

And the visible pass in 3D (visible pass occurs when the ground observer is in dark and space object is in sun light) can be reached from the link,

https://drive.google.com/file/d/14F6Q6MH2FCJmfes40f9W3-kjYU7av1Ee/view?usp=sharing.

4. References

[1] Space-track.org (fetched TLE for the initial orbit), https://www.space-track.org/auth/login

[2] Orekit (used for frame transformations and gravitational potential computations), https://www.orekit.org/

[3] CNES VTS (used for 3D visualisation), https://timeloop.fr/vts/

[4] Ashikhmin, Michael and Shirley, Peter, An Anisotropic Phong BRDF Model, Journal of Graphics Tools, 2000.

5. Code

The code used to generate above analysis is provided below.

# Call the required libraries
import seaborn as sns
import numpy as np
from matplotlib import pyplot
from state import State
from visibility import LambertianSphere # Replace ashikhminPremoze with LambertianSphere for cannonball
from forces import NeutralDrag
from forces import HolmesFeatherstoneGravity
from forces import ThirdBodyForce
from forces import ImprovedRadiation
from torques import Torques
from torques import GravityGradient
from torques import MagneticDipole
from propagator import NumericalPropagator
from sgp4.io import twoline2rv
from sgp4.earth_gravity import wgs84
import utils as utl

# Initial orbit is transformed from TLE mean elements in TEME to cartesian states in J2000.
TLE=[]
TLE.append('0 OPS 5111')
TLE.append('1 10684U 78020A   20110.24021950 +.00000083 +00000-0 +00000-0 0  9994')
TLE.append('2 10684 063.0042 120.7837 0079708 206.3069 153.3306 01.98073183291879')#Eccentricity modified
satellite = twoline2rv( TLE[1],TLE[2],wgs84)
tlepos = satellite.propagatemin(0)[0]
tlevel = satellite.propagatemin(0)[1]
tleTime = str(satellite.epoch.strftime("%Y-%m-%dT%H:%M:%S.%f")[0:-3])
state = np.array([tlepos[0],tlepos[1],tlepos[2],tlevel[0],tlevel[1],tlevel[2]])*1000.0
stateNew = utl.rotateFrame(tleTime,state,'TEME','J2000')

# The configuration for Steel cannonball all diffusive and HAMR object with diffuse and specular reflection
spacecraftConfig = {}
spacecraftConfig['name'] = 'Steel'
spacecraftConfig['type'] = '3DoF'#6DoF or 3DoF
spacecraftConfig['mass'] = 7850#kg
spacecraftConfig['date'] = tleTime#UTC time
spacecraftConfig['orbit'] = np.array([stateNew[0] ,stateNew[1] ,stateNew[2],\
                               stateNew[3] ,stateNew[4] ,stateNew[5]])#PV in m and m/s (J2000)
# Below is requred for 3DoF propagation - change type to 3DoF as well
spacecraftConfig['radius'] = 0.6204#m
spacecraftState = State(spacecraftConfig)
# # Below is the configuration for a HAMR object
# spacecraftConfig['name'] = 'Hamr'
# spacecraftConfig['type'] = '6DoF'#6DoF or 3DoF
# spacecraftConfig['mass'] = 0.0001125#kg
# spacecraftConfig['date'] = tleTime#UTC time
# spacecraftConfig['orbit'] = np.array([stateNew[0] ,stateNew[1] ,stateNew[2],\
#                                stateNew[3] ,stateNew[4] ,stateNew[5]])#PV in m and m/s (J2000)
# # Below is requred for 3DoF propagation - change type to 3DoF as well
# spacecraftConfig['radius'] = 0.3#m
# # Below is requred for 6DoF propagation - change type to 6DoF as well
# spacecraftConfig['attitude'] = np.array([20.0,20.0,0.0,0.0,0.0,0.5])#omega and omega dot w.r.t orbit frame
# spacecraftConfig['moi'] = np.array([[1.54180631e-07, 0.0, -1.08223763e-07],
#                                     [0.0, 2.31268163e-07, 0.0],
#                                     [-1.08223763e-07, 0.0, 8.41187813e-08]])#m^2 w.r.t body frame
# spacecraftConfig['magDipole'] = np.array([0.01,0.01,0.01])#Nm/Tesla = Am^2
# spacecraftConfig['comOffset'] = np.array([0.0,0.0,0.0])#m w.r.t body frame
# spacecraftConfig['centerOfPressure'] = np.array([[0.0,0.0,0.0558 ,0.0558 ,0.02165 ,0.02165  , \
#                                                   -0.02165,-0.02165,-0.0558  ,-0.0558],
#                                                  [0.0,0.0,0.0    ,0.0    ,0.0     ,0.0      , \
#                                                   0.0     ,0.0     ,0.0      ,0.0],
#                                                  [0.0,0.0,0.07165,0.07165,0.0375  ,0.0375   , \
#                                                   -0.0375 ,-0.0375 ,-0.07165 ,-0.07165]])\
#                                        # location of facet center w.r.t body frame  - meters
# spacecraftConfig['areaOfFacets'] = np.array([0.0025,0.0025,0.0025,0.0025,0.0025,0.0025,0.0025,0.0025, \
#                                              0.0025,0.0025])\
#                                    # area of each facet - meters^2
# spacecraftConfig['orientationOfFacets'] = np.array([[0.0,0.0  ,0.0  ,0.0  ,0.0 ,0.0  ,0.0 ,0.0  ,0.0 ,0.0  ],
#                                                     [0.0,180.0,30.0 ,210.0,60.0,240.0,60.0,240.0,30.0,210.0],
#                                                     [0.0,0.0  ,0.0  ,0.0  ,0.0 ,0.0  ,0.0 ,0.0  ,0.0 ,0.0  ]])\
#                                         # orientation of the facet w.r.t body frame - degrees
# spacecraftState = State(spacecraftConfig)

# Defines the airdrag parameters
dragCoeff = 2.8
dragNeutral = NeutralDrag(spacecraftState,dragCoeff)

# Defines the aspherical gravitational potential parameters
order = 20
degree = 20
gravity = HolmesFeatherstoneGravity(spacecraftState,order,degree)

# Defines the third-bosy mass interaction parameters 
sunBody = spacecraftState.sun
moonBody = spacecraftState.moon
thirdBodies = ThirdBodyForce(spacecraftState,[sunBody,moonBody])

# Defines the visibility parameters
gs = [57.0, -5.0, 606.0]#latitude, longitude, altitude
brdfCoeff = {}
brdfCoeff['cSunVis'] = 455#W/m^2
brdfCoeff['d'] = 0.5
brdfCoeff['rho'] = 0.26
brdfCoeff['s'] = 0.5
brdfCoeff['Fo'] = 0.5
brdfCoeff['nu'] = 0.5
brdfCoeff['nv'] = 0.5
# AshikhminPremoze = ashikhminPremoze(spacecraftState,gs,brdfCoeff)
# Below is when cannonball option
LambertianSphere = LambertianSphere(spacecraftState,gs,brdfCoeff)

# SRPClassical = ClassicalRadiation(spacecraftState,brdfCoeff)

# Initializes the SRP force class
emissivity = 0.05
SRPImproved = ImprovedRadiation(spacecraftState,brdfCoeff,emissivity)

# Initializes the torque classes
# aeroTorque = Torques(spacecraftState,dragNeutral)
# radiationTorque = Torques(spacecraftState,SRPImproved)
# ggTorque = GravityGradient(spacecraftState,[])#No need for a force Model
# magneticTorque = MagneticDipole(spacecraftState,[])#No need for a force Model

# Defines the propagator
#prop = NumericalPropagator(8640.0,100.0,spacecraftState,'dop853',\
#                           [dragNeutral,gravity,thirdBodies,SRPClassical],\
#                            [aeroTorque,radiationTorque,ggTorque,magneticTorque])
# prop = NumericalPropagator(86400.0,10.0,spacecraftState,'dop853',\
#                            [dragNeutral,gravity,thirdBodies,SRPImproved],\
#                             [aeroTorque,radiationTorque],ashikhminPremoze)
# prop = NumericalPropagator(86400.0,10.0,spacecraftState,'dop853',\
#                            [dragNeutral,gravity,thirdBodies,SRPImproved],\
#                             [],AshikhminPremoze)
prop = NumericalPropagator(86400.0,10.0,spacecraftState,'dop853',\
                           [dragNeutral,gravity,thirdBodies,SRPImproved],\
                            [],LambertianSphere)

# Propagate the orbit/attitude
prop.Propagate()

# Keep steel data
dataSteel = prop.sol['state'][:,:6]
# dataHamr  = prop.sol['state'][:,:6]

# Plots the flux seen by the ground observer
selected = prop.sol['visibility'][:,2]>np.zeros(len(prop.sol['visibility'][:]))
desired = np.where(selected)
startIndex = desired[0][0]
endIndex = desired[0][-1]
sns.set()
print('Below plot is flux when start time : ' +str(prop.sol['time'][startIndex]) + \
      ' and end time : '+str(prop.sol['time'][endIndex])+' (Cannonball object)')
pyplot.plot(prop.sol['visibility'][startIndex:endIndex,1])
pyplot.xlabel('Time (x10 seconds)')
pyplot.ylabel('Flux (W/m^2)')
pyplot.show()

Below plot is flux when start time : 2020-04-19T21:17:04.964 and end time : 2020-04-19T22:55:24.964 (Cannonball object)

# Plots the flux seen by the ground observer
selected = prop.sol['visibility'][:,2]>np.zeros(len(prop.sol['visibility'][:]))
desired = np.where(selected)
startIndex = desired[0][0]
endIndex = desired[0][-1]
sns.set()
print('Below plot is flux when start time : ' +str(prop.sol['time'][startIndex]) + \
      ' and end time : '+str(prop.sol['time'][endIndex])+' (HAMR object)')
pyplot.plot(prop.sol['visibility'][startIndex:endIndex,1])
pyplot.xlabel('Time (x10 seconds)')
pyplot.ylabel('Flux (W/m^2)')
pyplot.show()

Below plot is flux when start time : 2020-04-19T21:17:04.964 and end time : 2020-04-19T22:55:24.964 (HAMR object)

# Compute the relative orbit trajectory
lenData = dataHamr.shape[0]
relTrajectory = np.zeros((lenData,3))
posvel = np.zeros(6)
for cnt1 in range(lenData):
    posvel[0:3] = dataSteel[cnt1,:3]
    posvel[3:6] = dataSteel[cnt1,3:]
    rotMatIne2Orb   = utl.OrbitFrameAxes(posvel).T
    relTrajectory[cnt1,:] = np.dot(rotMatIne2Orb,dataHamr[cnt1,:3]-dataSteel[cnt1,:3])/1000.0

# Plots the relative orbit trajectory
pyplot.plot(relTrajectory[:,0],'r',markersize=0.1, label="Along-Track")
pyplot.plot(relTrajectory[:,1],'g',markersize=0.1, label="Nadir")
pyplot.plot(relTrajectory[:,2],'b',markersize=0.1, label="Neg. Orbit Normal")
pyplot.xlabel('Time (x10 seconds)')
pyplot.ylabel('Reative Distance (km)')
pyplot.legend(loc='center left', bbox_to_anchor=(1, 0.5))
pyplot.show()

# utl.createVTSSim('hamr', prop.sol['time'], prop.sol['state'][:,0:10], '6DoF')
utl.createVTSSim('steel', prop.sol['time'], prop.sol['state'][:,0:10], '3DoF')