Universe dynamics model evaluation example¶
In this tutorial we will set up different dynamics model and evaluate them on a spacecraft.
First, load all the godot modules we will need:
import numpy as np
from godot.core import tempo, util, autodif as ad
from godot.model import common
from godot import cosmos
Load the universe configuration and create the universe object.
# avoid verbose output
util.suppressLogger()
uni_config = cosmos.util.load_yaml("universe.yml")
uni = cosmos.Universe(uni_config)
We need to add a model that will store the spacecraft (or test point) state vector w.r.t. an arbitrary point in the frame system. We can change the value of the model as it is constant, and do evaluations.
translation_model = common.ConstantVectorTimeEvaluable(np.zeros(6))
print(translation_model.eval(tempo.Epoch()))
translation_model.set([1, 2, 3, 4, 5, 6])
print(translation_model.eval(tempo.Epoch()))
We can use this model to add a new point in the frame system that is translated from Earth along the ICRF axes.
uni.frames.addPoint("Spacecraft", tempo.TimeScale.TDB)
uni.frames.addTimeEvaluableTranslation("Spacecraft", "Earth", "ICRF", translation_model)
Similarly to before, we can use the frame system to evaluate the state vector w.r.t. any frame system point and along an arbitrary axes.
print(uni.frames.vector6("Earth", "Spacecraft", "ICRF", tempo.Epoch()))
translation_model.set([6, 5, 4, 3, 2, 1])
print(uni.frames.vector6("Earth", "Spacecraft", "ICRF", tempo.Epoch()))
We can retrieve the dynamics models from the universe, but evaluation is not psosible at the moment because the dynamics model does not know about the test point the evaluation should use.
dyn = uni.dynamics.get("EarthPointMass")
try:
dyn.acc.eval(tempo.Epoch())
except RuntimeError as e:
print(e)
We need to alias the spacecraft point such that the dynamics model will use that as test point (with zero mass).
uni.frames.setAlias(dyn.point, uni.frames.pointId("Spacecraft"))
print(dyn.acc.eval(tempo.Epoch()))
translation_model.set([7000.0, 0, 0, 0, 0, 0])
print(dyn.acc.eval(tempo.Epoch()))
We can similarly interate over the available dynamics models and evaluate different ones after making sure the frame point aliasing is done before each model evaluation. Here we also need to alias the evaluation axes for the combined dynamics model.
model_names = ["EarthPointMass", "EarthJ2", "EMS", "EMS_and_SRP"]
for name in model_names:
dyn = uni.dynamics.get(name)
uni.frames.setAlias(dyn.point, uni.frames.pointId("Spacecraft"))
uni.frames.setAlias(dyn.axes, uni.frames.axesId("ICRF"))
print(name, dyn.acc.eval(tempo.Epoch()))
In case the user needs to compute the partials derivatives w.r.t the state of the spacecraft, it can be done by using the eval method with an XEpoch.
print(translation_model.eval(tempo.XEpoch()))
To actually get partial derivatives, the model needs to be set to use autodif variable.
state = [7000.0, 0, 0, 0, 0, 0]
x = ad.Vector(state, "x")
translation_model.set(x)
print(translation_model.eval(tempo.XEpoch()))
Then, we can use the dynamics model the same way as before, with XEpoch and get the partials.
dyn = uni.dynamics.get("EarthPointMass")
uni.frames.setAlias(dyn.point, uni.frames.pointId("Spacecraft"))
a = dyn.acc.eval(tempo.XEpoch())
print(a)
Gradient matrix can be retrieved to get a numpy array.
np.set_printoptions(linewidth=100)
print(a.gradient()["x"])