Propagating Massive Bodies

This notebook shows how to propagate massive bodies using godot. This will be done by simulating the Didymos system consisting of an asteroid with a moon.

This notebook also shows how the Asteroid system propagation is performed in Python, and then the propagated asteroid orbits and attitudes are used in the Hera spacecraft Trajectory propagation and subsequent optimsiation.

First we import the necessary libraries.

In [1]:
from godot.core import tempo, util
from godot.model import prop, common
from godot import cosmos
import yaml
import numpy as np
from aux import plotting, summary, geometry

util.suppressLogger()

Then we set up the universe. This is where most of the configuration takes place. In the configuration file didymos_universe.yml, we do several things:

  • In constants, we set the GM and the SOI (sphere of influence) of both Didymos and Didymoon.
  • In frames, we define the initial orbits of both asteroids (twice!). For Didymos the orbit around the sun, and for Didymoon the orbit around Didymos. One of these we will overwrite with our propagation in Python.
  • Also in frames, we define the orientation of the asteroids as AxesLocalOrbital/Pos frames.
  • In parameters we define a parameter correction to the constant GM of both bodies.
  • In sphericalHarmonics, we define some spherical harmonics for each twin using tables files.
  • In gravity, we define 3 different systems. One just with Didymos, one just with Didymoon and one with both. Didymos needs the gravity of just Didymoon, Didymoon of just Didymos and the spacecraft the gravity of both.
  • In bodies, we then link the gravity systems with the propagated bodies, defining the GMs as the sum of the constant and the parameter GM.
  • In dynamics, the gravity systems are mapped to dynamics so they can be used by the spacecraft.
In [2]:
# set up universe
with open("didymos_universe.yml", "r") as f:
    uniConfig = yaml.safe_load(f)
uni = cosmos.Universe(uniConfig)

Now we can load the spacecraft trajectory file

In [3]:
# set up trajectory and compute
traConfig = cosmos.util.load_yaml("hera_trajectory_opt.yml")
tra = cosmos.Trajectory(uni, traConfig)

Performing the Didymos System Orbit/Attitude Propagation

We will propagate the Didymos System, that is orbit and attitude of both bodies. The orbit propagation will use the above systems, the attitude will be a torque free propagation of the axes.

In [4]:
# generate a new model for didymoon
didymoon_point = uni.frames.pointId("didymoon")
# generate a new model for didymos
didymos_point = uni.frames.pointId("didymos")
# generate a new model for didymoon
didymoon_axes = uni.frames.axesId("didymoonAxes")
# generate a new model for didymos
didymos_axes = uni.frames.axesId("didymosAxes")

Axes Propagation Setup

First we setup the torque free propagation of the asteroids.

For didymoon, we get the initial quaternion for the bodies from the nominal models from the universe. Since it is probably tidally locked this is a reasonable orientation model.

In [5]:
# get inertial axes
icrf = uni.frames.axesId("ICRF")

# define the initial state
epoch0 = tempo.parseEpoch("2020-09-10T00:00:00 TDB")

# get the initial quaternion for didymoon from the universe nominal models (LocalOrbital/Pos)
didymoon_q0 = uni.frames.quaternionState(icrf, uni.frames.axesId("didymoonLOFAxes"), epoch0)

For didymos, we introduce a completely new orientation model and overwrite the nominal model with this one.

In [6]:
from aux import geometry
alpha = uni.parameters.add('alpha', np.pi/4)
delta = uni.parameters.add('delta', 0.)
omega = uni.parameters.add('omega', 0.)
vAngular = np.array([0., 0.5, 1.])
vAngular = 2 * np.pi* vAngular / (86400 * np.linalg.norm(vAngular))
alphaDot = uni.parameters.add("vAlpha", vAngular[0])
deltaDot = uni.parameters.add("vDelta", vAngular[1])
omegaDot = uni.parameters.add("vOmega", vAngular[2])
didymos_q0 = geometry.Quaternion(alpha, delta, omega, alphaDot, deltaDot, omegaDot)

next we define shape/inertia models for the bodies. We define these as ellipoids:

In [7]:
from aux import shape

# add a density parameter
rho = uni.parameters.add("density", 1e12)

# create 2 instances of ellipsoids for the two bodies
didymos_ellipsoid = shape.Ellipsoid(3, 4, 5, rho)
didymoon_ellipsoid = shape.Ellipsoid(0.5, 0.4, 0.5, rho)

finally we create the propagator inputs for the axes.

In [8]:
# specify the torque on the asteroids to be zero.
torque = common.ConstantVector3TimeEvaluable(np.zeros(3))

# set up propagator axes/input for didymos
didymos_prop_axes = prop.PropagatorAxes(uni.frames, didymos_axes, icrf)
didymos_input_axes = prop.AxesInput(didymos_prop_axes, didymos_q0, torque, didymos_ellipsoid.inertia(), np.full(7, 1e-9), np.full(7, 1e-9))

# set up propagator axes/input for didymoon
didymoon_prop_axes = prop.PropagatorAxes(uni.frames, didymoon_axes, icrf)
didymoon_input_axes = prop.AxesInput(didymoon_prop_axes, didymoon_q0, torque, didymoon_ellipsoid.inertia(), np.full(7, 1e-9), np.full(7, 1e-9))

Point Propagator Setup

Now we repeat the above, this time to setup the propagator points.

We get the initial states from the nominal models in the universe.

In [9]:
sun = uni.frames.pointId("Sun")

# get the initial states
didymos_s0 = uni.frames.vector6(sun, uni.frames.pointId("didymosKEP"), icrf, epoch0)
didymoon_s0 = uni.frames.vector6(sun, uni.frames.pointId("didymoonKEP"), icrf, epoch0)

These propagations use dynamic models from the universe, which need to have some aliases set to make them work.

In [10]:
didymos_dynamics = uni.dynamics.get("didymos")
didymoon_dynamics = uni.dynamics.get("didymoon")
for dyn in [didymos_dynamics, didymoon_dynamics]:
    uni.frames.setAlias(dyn.coi, sun)
    uni.frames.setAlias(dyn.axes, icrf)

Now we can set up the point propagator inputs

In [11]:
# create the didymos propagation input
didymos_prop_point = prop.PropagatorPoint(uni.frames, icrf, didymos_point)
didymos_input_point = prop.PointInput(didymos_prop_point, sun, didymos_s0, didymos_dynamics.acc, np.full(6, 1e-9), np.full(6, 1e-9))

# create the didymoon propagation input
didymoon_prop_point = prop.PropagatorPoint(uni.frames, icrf, didymoon_point)
didymoon_input_point = prop.PointInput(didymoon_prop_point, sun, didymoon_s0, didymoon_dynamics.acc, np.full(6, 1e-9), np.full(6, 1e-9))

We retrieve the initial state of the twins and use it to propagate them.

In [12]:
# collect the propagator inputs
inputs = [didymos_input_point, didymoon_input_point, didymos_input_axes, didymoon_input_axes]

dt = common.ConstantScalarTimeEvaluable(0.)

# create the system propagator
system_prop = prop.Propagator("", epoch0, dt, inputs)

# propagate
epoch1 = tempo.parseEpoch("2020-10-31T00:00:00 TDB")
system_prop.compute(epoch1, True)

Now we can plot the orientations of the two bodies.

In [13]:
plotting.plot_axes(uni, uni.frames.axesId("didymosAxes"),uni.frames.axesId("ICRF"),tempo.EpochRange(epoch0,epoch1))
plotting.plot_axes(uni, uni.frames.axesId("didymoonAxes"),uni.frames.axesId("ICRF"),tempo.EpochRange(epoch0,epoch1))

Spacecraft Propagation

Now we also want to add a spacecraft that flies by the two bodies. This is done using the trajectory class.

In [14]:
uni.parameters.track(uni.parameters.contains())
tra.compute(False)

We print out the initiali state from the trajectory in various frames and representations

In [15]:
start = tempo.Epoch("2020-09-27T00:00:00 TDB")
end = tempo.Epoch("2020-10-10T00:00:00 TDB")
box=50
plotting.plot_points(uni,uni.frames.pointId("didymoon"),uni.frames.pointId("HERA_center"),uni.frames.pointId("didymos"), uni.frames.axesId("ICRF"),tempo.EpochRange(start,end), box=box )
print("Sun-Hera ICRF state @",start)
s = uni.frames.vector6(uni.frames.pointId("Sun") , uni.frames.pointId("HERA_center"), uni.frames.axesId("ICRF"), start)
print("pos_x: {:.12e} km".format(s[0]))
print("pos_y: {:.12e} km".format(s[1]))
print("pos_z: {:.12e} km".format(s[2]))
print("vel_x: {:.15e} km/s".format(s[3]))
print("vel_y: {:.15e} km/s".format(s[4]))
print("vel_z: {:.15e} km/s".format(s[5]))
summary.trajectory(uni,tra)

Sun-Hera ICRF state @ 2020-09-27T00:00:00.000000 TDB
pos_x: 1.073895129565e+08 km
pos_y: 1.083911171453e+08 km
pos_z: -4.260835673821e+06 km
vel_x: -2.298947555279856e+01 km/s
vel_y: 2.584223723156880e+01 km/s
vel_z: 1.754179042507382e+00 km/s
launch 	 2020-09-27T00:00:00.000000 TDB 	 0.0
man1_end 	 2020-09-29T00:00:00.000000 TDB 	 3.548058012515837e-06
man1 	 2020-09-29T00:00:00.000000 TDB 	 3.548058012515837e-06
man1_start 	 2020-09-29T00:00:00.000000 TDB 	 3.548058012515837e-06
match1_left 	 2020-09-30T00:00:00.000000 TDB 	 0.00014931378167770988
match1_right 	 2020-09-30T00:00:00.000000 TDB 	 0.00014931378167770988
man2_end 	 2020-10-01T00:00:00.000000 TDB 	 7.958084385626637e-13
man2 	 2020-10-01T00:00:00.000000 TDB 	 7.958084385626637e-13
man2_start 	 2020-10-01T00:00:00.000000 TDB 	 7.958084385626637e-13
arrival 	 2020-10-03T00:00:00.000000 TDB 	 0.0
end 	 2020-10-10T00:00:00.000000 TDB 	 4.950215345172641e-13