Adding a dynamic model to a propagation¶

Here we will add an acceleration to an already existing universe object.

The acceleration can be added to a dynamics model of type combined in the universe:

yml
dynamics:
  - name: GravityOnly
    type: SystemGravity
    config:
      model: SunEarthMarsGravity
      switch: true
  - name: SunEarthMarsDynamics
    type: Combined
    config:
      - GravityOnly

In this case the model SunEarthMarsDynamics can be added to, but GravityOnly cannot.

First perform all necessary imports

In [1]:

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

from godot.core import tempo, util, constants
from godot.core import autodif as ad
from godot.core.autodif import bridge as br
from godot.model import common
from godot import cosmos

next load the universe.

In [2]:

# avoid verbose output
util.suppressLogger()
uni_config = cosmos.util.load_yaml('universe.yml')
uni = cosmos.Universe(uni_config)

Note that in the universe configuration the dynamics model named SunEarthMarsDynamics is 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

In [3]:

tra_config = cosmos.util.load_yaml('trajectory.yml')
tra = cosmos.Trajectory(uni, tra_config)

Now we plot the trajectory

In [4]:

def plot():
    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):
        try:
            return uni.frames.vector3('Sun', 'SC_center', 'EMC', e) / constants.AU
        except:
            return [np.nan] * 3

    grid = tra.range().createGrid(1.0 * tempo.SecondsInDay)
    t = [e.mjd() for e in grid]
    x = np.asarray([pos(e) for e in grid])

    plt.figure(figsize=(8, 8))
    
    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.plot(x[:, 0], x[:, 1], '-k', linewidth=2, label="Transfer")
    plt.plot(x[0, 0], x[0, 1], 'ok')
    plt.plot(x[-1, 0], x[-1, 1], 'ok')
    plt.xlim([-2, 2])
    plt.ylim([-2, 2])
    plt.legend()
    plt.grid()

In [5]:

tra.compute(False)
plot()

Now we can add a new dynamic model. We create a common.Vector3TimeEvaluble class and compute a Sun, SC radial acceleration of fixed magnitude.

In [6]:

from typing import Union
class ConstantRadialAcceleration(common.Vector3TimeEvaluable):
    def __init__(self, uni, mag):
        # We have to call C++ trampoline class constructor explicitly
        common.Vector3TimeEvaluable.__init__(self)
        self.mag = mag
        self.uni = uni
    def eval(self, epoch: Union[tempo.Epoch, tempo.XEpoch]) -> Union[np.ndarray, ad.Vector]:
        pos = self.uni.frames.vector3("Sun", "SC_center", "ICRF", epoch)
        if isinstance(self.mag, float):
            m = self.mag
        else:
            m = self.mag.eval(epoch)
        acc = pos / br.norm(pos) * m
        return acc

Next we construct an object, get the combined SunEarthMarsDynamics dynamics model and add our acceleration model.

In [7]:

# Construct a constant radial acceleraiton
acc_radial = ConstantRadialAcceleration(uni, 1e-7)
# Get the dynamics models (must be a combined model)
dyn = uni.dynamics.get("SunEarthMarsDynamics")
# Add the constant acceleration to the combined dynamics 
dyn.acc.add(acc_radial)

Now we can propagate the new trajectory

In [8]:

tra.compute(False)
plot()

Of course we can create and pass a parameter to the acceleration model for use in optimisation and estimation:

In [9]:

mag = ad.Scalar(1e-7, 'mag')
magFunc = common.ConstantScalarTimeEvaluable(mag)

# Construct a constant radial acceleraiton
acc_radial_gradients = ConstantRadialAcceleration(uni, magFunc)
# Get the dynamics models (must be a combined model)
dyn = uni.dynamics.get("SunEarthMarsDynamics")
# Add the constant acceleration to the combined dynamics 
dyn.acc.add(acc_radial_gradients)

Now, if we compute the trajectory we will see that the gradient with respect to the parameter mag is available for use.

In [10]:

tra.compute(True)

pos = uni.frames.vector3("Sun", "SC_center", "ICRF", tempo.XEpoch("2026-12-31 TDB"))
print(common.resolve(pos))

         Value |           mag
  3.084142e+05 |  4.444340e+13
  1.694370e+08 |  3.220936e+13
  7.345317e+07 |  1.396576e+13