Adding Action/Events to a Propagator¶
We will perform a propagation in a number of different ways, first using the trajectory and then using the BallisticPropagator and finally the standard Propagator.
We will show how to use events with triggers and actions to perform user defined actions on the propagation.
First perform all necessary imports
from godot import cosmos
from godot.model import prop
from godot.model import common
from godot.model import interface
from godot.model import geometry
from godot.core.autodif import bridge as br
from godot.core import autodif as ad
from godot.core import tempo, util, constants, integ
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
# avoid verbose output
util.suppressLogger()
next load the universe.
uni_config = cosmos.util.load_yaml('universe.yml')
uni = cosmos.Universe(uni_config)
Note that in the universe configuration the dynamics model named
SunEarthMarsDynamicsis of type Combined. Only Combined dynamics models can have more dynamics models added to them:
Now load the trajectory configuration and create the trajectory object using the universe object
tra_config = cosmos.util.load_yaml('trajectory.yml')
tra = cosmos.Trajectory(uni, tra_config)
Now we plot the trajectory
def plot(points):
ran = tempo.EpochRange(tempo.Epoch("2020-01-01 TDB"),
tempo.Epoch("2022-01-01 TDB"))
earth = np.asarray([uni.frames.vector3('Sun', 'Earth', 'EMC', e) /
constants.AU for e in ran.createGrid(1.0 * tempo.SecondsInDay)])
mars = np.asarray([uni.frames.vector3('Sun', 'Mars', 'EMC', e) /
constants.AU for e in ran.createGrid(1.0 * tempo.SecondsInDay)])
def pos(e,point):
try:
return uni.frames.vector3('Sun', point, 'EMC', e)
except:
return [np.nan] * 3
plt.figure(figsize=(8, 8))
grid = tra.range().createGrid(0.01 * tempo.SecondsInDay)
t = [e.mjd() for e in grid]
for point in points:
x = np.asarray([pos(e, point) for e in grid])/constants.AU
plt.plot(x[:, 0], x[:, 1], '-', linewidth=2, label=point)
plt.plot(x[0, 0], x[0, 1], 'o')
plt.plot(x[-1, 0], x[-1, 1], 'o')
plt.xlabel('EMC X (AU)')
plt.ylabel('EMC Y (AU)')
plt.plot(earth[:, 0], earth[:, 1], '--', label="Earth")
plt.plot(mars[:, 0], mars[:, 1], '--', label="Mars")
plt.xlim([-2, 2])
plt.ylim([-2, 2])
plt.legend()
plt.grid()
if len(points) > 1:
plt.figure(figsize=(8, 4))
ref = points[0]
for point in points[1:]:
x = np.asarray([np.linalg.norm(pos(e, point)-pos(e,ref)) for e in grid])
plt.plot(t, x, '-', linewidth=2, label=f"|{point}-{ref}|")
plt.legend()
plt.grid()
tra.compute(False)
plot(['SC_center'])
Ballistic Propagator¶
We create a cosmos.BallisticPropagator to propagate the spacecraft again with the same dynamic model. The difference with the trajectory is that in this case only the state is propagated and in the trajectory, state mass and deltav ar propagated.
We can create a simple action to print the epoch of the trigger time, but with this propagator no action can be taken to adjust the state.
The difference in the propagator is
epoch0 = tempo.Epoch("9750.0 TDB")
epoch1 = tempo.Epoch("9755.0 TDB")
# create the BallisticPropagator
pro = cosmos.BallisticPropagator(
uni, 'SC', 'SunEarthMarsDynamics', epoch0, 'Earth', 1e-12)
# Action to print output only
class ExampleAction(interface.EventAction):
def __init__(self):
# We have to call C++ trampoline class constructor explicitly
interface.EventAction.__init__(self)
# Assign member variables
def eval(self, epoch, crossDir, integDir):
def actionFunc(): return integ.Action.Continue
print(f"epoch={epoch}")
return interface.create_action(epoch, actionFunc, dict())
action = ExampleAction()
# create the trigger
trigger = common.SimpleFixedEventTrigger(
[tempo.Epoch("9751 TDB"), tempo.Epoch("9752 TDB")])
# create the event and add it to the propagator
event = common.SimpleEvent(trigger, action)
pro.addEvent(event)
# propagate state from trajectory initial state
x0 = uni.frames.vector6('Earth', 'SC_center', 'ICRF', epoch0)
pro.compute(x0, 0, epoch1)
# plot results
plot(['SC_center','SC'])
Propagator¶
Lets now add a spacecraft propagator to the universe. This time we will use a different dynamic model for the propagation and compare the results.
First get variables for the points we need and create a point for the propagator to use:
# Get frames object and variables for points to use
fra = uni.frames
icrf = fra.axesId('ICRF')
earth = fra.pointId('Earth')
# Create point for spacecraft
try:
scPoint = fra.addPoint('Spacecraft', tempo.TimeScale.TDB)
except:
scPoint = fra.pointId('Spacecraft')
Now we create a propagator
propPoint = prop.PropagatorPoint(fra, icrf, scPoint)
And define the dynamics (GravityOnly, so without SRP). We need to set some aliases to use the dynamics with the newly created spacecraft point:
dyn = uni.dynamics.get("GravityOnly")
fra.setAlias(dyn.point, scPoint)
fra.setAlias(dyn.coi, earth)
fra.setAlias(dyn.axes, icrf)
for key, alias in dyn.inputs.items():
alias.set(uni.parameters.get("{}_{}".format("Spacecraft", key)))
Now we create a propagator and propagate and plot the difference with the propagation with SRP.
# We create a time evaluable to get the initial state from the SC_center already propagated
x0 = geometry.Vector6(fra, 'Earth', 'SC_center', 'ICRF')
# Define error tolerance
tol = np.full(6, 1e-12)
# Setup inputs for propagator
inp = [prop.PointInput(propPoint, earth, x0, dyn.acc, tol, tol)]
# create a propagator
pro = prop.Propagator('Prop', epoch0, None, inp)
# Propagate without partials
pro.compute(epoch1, False)
# Output final state and spacecraft X axis direction
plot(['SC_center','Spacecraft'])
This shows the effect of SRP on the propagation. Suppose we want to correct the state of the spacecraft to account for the SRP, using the trajectory reference as a guide. We do this by adding events to the propagator, correcting the Spacecraft state to match the SC_center values.
We define our correction action:
class CorrectionAction(interface.EventAction):
def __init__(self, propPoint):
# We have to call C++ trampoline class constructor explicitly
interface.EventAction.__init__(self)
# Assign member variables
self.__propPoint = propPoint
def eval(self, e, crossDir, integDir):
"""Correct the state back to the SC_center value
Parameters
----------
e : tempo.Epoch
The epoch
crossDir : [type]
Cross direction as defined by the event
integDir : [type]
Integration direction
Returns
-------
[type]
[description]
"""
delta = uni.frames.vector6('Spacecraft', 'SC_center', 'ICRF', e)
print(f"adding correction at {e}={delta}")
def actionFunc(): return integ.Action.Continue
deltaState = dict()
deltaState[self.__propPoint] = delta
return interface.create_action(e, actionFunc, deltaState)
We add the event with a simple trigger occurring every hour to our existing propagator:
action = CorrectionAction(propPoint)
trigger = common.SimpleFixedEventTrigger(tempo.EpochRange(epoch0,epoch1).createGrid(3600))
event = common.SimpleEvent(trigger, action)
pro.addEvent(event)
Finally we can re-propagate:
# Propagate without partials
pro.compute(epoch1, False)
# Output final state and spacecraft X axis direction
plot(['SC_center','Spacecraft'])
def Propagator(_uni, _name, _dyn, _coi, _x0, _epoch0, _epoch1, _dt, _delta, _tol=np.full(6, 1e-12)) :
fra = _uni.frames
icrf = fra.axesId('ICRF')
coi = fra.pointId(_coi)
# get or create sc point
if fra.containsPoint(_name):
scPoint = fra.pointId(_name)
else:
scPoint = fra.addPoint(_name, tempo.TimeScale.TDB)
# set up dynamic model
if not _uni.dynamics.contains(_dyn):
raise ValueError(f'Dynamics model {_dyn} not present in universe')
dyn = _uni.dynamics.get(_dyn)
fra.setAlias(dyn.point, scPoint)
fra.setAlias(dyn.coi, earth)
fra.setAlias(dyn.axes, icrf)
for key, alias in dyn.inputs.items():
alias.set(_uni.parameters.get("{}_{}".format(_name, key)))
# create the propagator
propPoint = prop.PropagatorPoint(fra, icrf, scPoint)
inp = [prop.PointInput(propPoint, coi, _x0, dyn.acc, _tol, _tol)]
pro = prop.Propagator('Prop', _epoch0, None, inp)
class Action(interface.EventAction):
def __init__(self, propPoint):
interface.EventAction.__init__(self)
self.__propPoint = propPoint
def eval(self, e, crossDir, integDir):
delta = _delta(e)
def actionFunc(): return integ.Action.Continue
deltaState = dict()
deltaState[self.__propPoint] = delta
return interface.create_action(e, actionFunc, deltaState)
action = Action(propPoint)
epochs = tempo.EpochRange(_epoch0,_epoch1).createGrid(_dt)
trigger = common.SimpleFixedEventTrigger(epochs)
event = common.SimpleEvent(trigger, action)
pro.addEvent(event)
pro.compute(_epoch1, False)
x0 = geometry.Vector6(fra, 'Earth', 'SC_center', 'ICRF')
if not fra.containsPoint('mysc'):
scPoint = fra.addPoint('mysc', tempo.TimeScale.TDB)
# define some action function
def func(e): return fra.vector6('mysc', 'SC_center', 'ICRF', e) /2
Propagator(uni, 'mysc', 'GravityOnly', 'Earth', x0, epoch0, epoch1, 3600,func)
plot(['SC_center','mysc'])
uni.frames.containsPoint('SC')