Orbit Determination Example¶
In this tutorial we will perform a basic orbit determination, using user defined data.
First, load all the godot modules we will need:
from godot.core import util, tempo, autodif as ad
from godot import cosmos
import godot.core.astro as astro
import godot.core.events as events
import godot.model.common as common
import godot.model.obs as obs
Load the universe configuration and create the universe object.
We print out all the parameters and their status. Notice all are fixed at this time.
# avoid verbose output
util.suppressLogger()
uni_cfg = cosmos.util.load_yaml("universe.yml")
uni = cosmos.Universe(uni_cfg)
print(uni.parameters)
We have created a simple trajectory for a Lunar orbiter with a single control point. We load this trajectory configuration and create the trajectory object.
Again, we print out all the parameters and their status. Notice the trajectory has created more parameters, but all are still fixed.
tra_cfg = cosmos.util.load_yaml("trajectory.yml")
tra = cosmos.Trajectory(uni, tra_cfg, True)
print(uni.parameters)
We have created a simple od problem for our Lunar orbiter. We load this problem configuration and create the problem object.
pro_cfg = cosmos.util.load_yaml("problem.yml")
prob = cosmos.orb.Problem(uni,pro_cfg)
Once the problem is created, the track-status of the parameters may be changed to free if they are to be estimated or considered.
print(uni.parameters)
Now we can propagate the trajectory, since we have established - via the problem configuration - which parameters to track
tra.compute(True)
First, we will perform a very simple OD solution, using only the a priori information which is specified in the problem configuration:
# create the solver from the problem
solv=cosmos.orb.Solver(prob)
# add equations from the problem configuration to the solver
solv.processProblemEquations()
# solve
solv.solve()
# generate a filter solution
fs = solv.filterSolution()
# print the solution and covariance
print(f"The solution {fs.solution()}")
print(f"The solution covariance {fs.covariance()}")
Adding tracking data¶
We can add equations to the filter, which means adding either a priori information or observations.
The user can add arbitrary equations by creating their own classes.
We use the OneWayLightTime class to build a simple Range observable:
class Range(common.VectorTimeEvaluable):
"""
Range Data
"""
def __init__(self, uni, station, spacecraft, ul_bias=None, dl_bias=None, xpr_bias=None, noise=None):
common.VectorTimeEvaluable.__init__(self)
self.__uni = uni
self.__station_point = self.__uni.frames.pointId(station)
self.__spacecraft_point = self.__uni.frames.pointId(spacecraft)
# optional biases and noise - not implemented below
self.__ul_bias = ul_bias
self.__dl_bias = dl_bias
self.__xpr_bias = xpr_bias
self.__noise = noise
self.__owlt = obs.OneWayLightTime(
self.__uni.frames,
self.__uni.frames.pointId("SSB"), # Light Time center
self.__uni.frames.axesId("ICRF"), # Light Time axes
maxIter = 5,
thdIter = 1e-6
)
def eval(self, epoch):
lt_downlink = self.__owlt.eval(
epoch,
self.__station_point, # receiver
self.__spacecraft_point, # transmitter
obs.Backward
)
if self.__dl_bias is not None:
lt_downlink += self.__dl_bias.eval(epoch)
if self.__xpr_bias is not None:
lt_downlink += self.__xpr_bias.eval(epoch - lt_downlink)
lt_uplink = self.__owlt.eval(
epoch - lt_downlink,
self.__spacecraft_point, # transmitter
self.__station_point, # receiver
obs.Backward
)
if self.__ul_bias is not None:
lt_uplink += self.__ul_bias.eval(epoch)
noise = self.__noise.eval(epoch) if self.__noise is not None else 0
return ad.concatenate([lt_downlink + lt_uplink + noise])
We can evaluate this as follows:
nno_rng = Range(uni, 'New_Norcia', 'SC_center')
o = nno_rng.eval(tempo.XEpoch("2026-01-07T16:42:13.815276 TDB"))
print(o)
print(common.resolve(o))
We can also compute random noise as follows:
import numpy as np
class GaussianNoise(common.ScalarTimeEvaluable):
"""Simple Gaussian noise"""
def __init__(self, std_dev, bias=0.0):
"""Constructor"""
super().__init__()
self.__bias = bias
self.__std_dev = std_dev
def eval(self, e) -> np.ndarray:
"""Eval method"""
return np.random.normal(self.__bias, self.__std_dev, 1)[0]
Observation Schedules¶
We can use the events system to compute when the observations can be made
def sc_nno_elev( epo ):
x = uni.frames.vector3('New_Norcia', 'SC_center', 'New_Norcia', epo)
return astro.sphericalFromCart(x)[2]
grid = tra.range().createGrid(86400/2)
sc_visible = events.generateEventIntervalSet( sc_nno_elev, 1e-6, grid, 1e-6 )
# TBD add occultation events and merge with sc_visible
# see https://gitlab.esa.int/midas/midas/-/blob/master/midas/od/_scheduler_filters.py#L307
# ...
# sc_occulted = events.generateEventIntervalSet( func, 1e-6, grid, 1e-6 )
# sc_not_occulted = events.complement(sc_occulted)
# get overlap of all
#sc_visible = events.overlap([sc_occulted, sc_visible])
# create observation times during visibility
time_tags = []
for pss in sc_visible:
for tt in pss.createGrid(600):
time_tags.append(tt.value())
Next we compute the True trajectory (in this case a perturbed initial state from the pre-loaded trajectory)
tra_true_cfg = cosmos.util.load_yaml("trajectory_true.yml")
tra_true = cosmos.Trajectory(uni, tra_true_cfg, True)
tra_true.compute(False)
Now we can process these observations, weighted by a noise value:
# create a new solver
solv=cosmos.orb.Solver(prob)
# add apriori information from the problem
solv.processProblemEquations()
# range noise std dev in seconds
range_noise = 0.33e-8
# create a modelled range object
rng_comp = Range(uni, "New_Norcia", "SC_center")
# create a observed range object
rng_obs = Range(uni, "New_Norcia", "SC_true_center", noise=GaussianNoise(range_noise))
# create residuals
# define a reference epoch
epoch = tempo.Epoch("2026-01-07T16:42:13.815276 TDB")
# process one hour of range data with 60 second time step
residuals0 = []
for ep in time_tags:
residual = rng_obs.eval(ep) - rng_comp.eval((tempo.XEpoch(ep)))
# process weighted residuals
solv.process(residual/range_noise)
residuals0.append([ep, residual.value()[0]])
# solve the od
solv.solve()
# create a solution object
fs = solv.filterSolution()
print(fs.covariance())
solv.update()
tra.compute(True)
residuals1 = []
for ep in time_tags:
residual = rng_obs.eval(ep) - rng_comp.eval(ep)
residuals1.append([ep, residual.value()[0]])
GODOT has powerful functionality for transforming and propagating the solution and covariance, for example:
class Spherical(common.VectorTimeEvaluable):
"""Class for converting Cartesian to Spherical Polar Coordinates
"""
def __init__(self, uni, spacecraft):
common.VectorTimeEvaluable.__init__(self)
self.__uni = uni
self.__spacecraft_point = spacecraft
def eval(self, epoch):
x = self.__uni.frames.vector3('Earth', self.__spacecraft_point, 'ITRF', epoch)
return astro.sphericalFromCart(x)
spherical = Spherical(uni, 'SC_center')
# now compute the covariance in spherical coordinates
epoch = tempo.Epoch("2026-01-07T16:42:13.815276 TDB")
for ep in tempo.EpochRange(epoch,epoch+86400/2).createGrid(3600):
print(fs.covariance(spherical, tempo.XEpoch(ep)))
Load the universe configuration and create the universe object
%matplotlib inline
import matplotlib.pyplot as plt
plt.plot([e.mjd() for e,x in residuals0], [x for e,x in residuals0],'.',label='pre-fit')
plt.plot([e.mjd() for e,x in residuals1], [x for e,x in residuals1],'.',label='post-fit')
plt.legend()
