Setup¶
Import basic python and additional godot packages.
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
Import basic godot packages regarding automatic differentiation.
from godot.core import num, tempo
from godot.core import autodif as ad
from godot.core.autodif import bridge as br
Import model packages for propagation.
from godot.model import common, prop
# optionally avoid verbose logging messages
import godot.core.util as util
util.suppressLogger()
Definition of differential equation¶
Here we define a class that implement the right-hand-side for integrating the cosine function to yield the sine function in the end. We parametrize the equation with the frequency of the sine wave:
f = sin(w * t)
The class has to implement the common.ScalarTimeEvaluable interface to be able to interact with the propagator.
class Equation(common.ScalarTimeEvaluable):
def __init__(self, epoch0, freq):
# We have to call C++ trampoline class constructor explicitly
common.ScalarTimeEvaluable.__init__(self)
# Assign member variables
self.__epoch0 = epoch0
self.__freq = freq
def eval(self, epoch):
# Computes derivative of d/dx sin(w * t) = w * cos(w * t)
t = epoch - self.__epoch0
w = self.__freq.eval(epoch)
return w * br.cos(w * t)
Script to setup propagator¶
We need to carefully setup the propagator inputs, initial states and differential equations for each thing(tm) we want to propagate. In this example there is only one scalar quantity.
sca = prop.PropagatorScalarTimeEvaluable()
First, we setup the initial conditions, epoch and state. The state is created as a time evaluable object as required by the propagator input.
epoch0 = tempo.Epoch(0.0, tempo.TimeScale.TDB)
x0 = ad.Scalar(0.0, 'x0')
x0Func = common.ConstantScalarTimeEvaluable(x0)
We setup another parameter for the frequency of the sine functions, which is defined in 1/s = Hz. We choose a frequency equivalent to one day.
freq = ad.Scalar(1.0 / 86400.0, 'freq')
freqFunc = common.ConstantScalarTimeEvaluable(freq)
Finally, we create the RHS object.
rhs = Equation(epoch0, freqFunc)
With these we can create the propagator input that defines the quantity to propagate. We also need to define the relative and absolute propagation error tolerances.
rtol = 1e-9
atol = 1e-9
inp = [prop.ScalarInput(sca, x0Func, rhs, rtol, atol)]
At this point the propagator can be created using this input. We also define a parameter for the delta in initial epoch.
dt0 = ad.Scalar(0.0, 'dt0')
dt0Func = common.ConstantScalarTimeEvaluable(dt0)
propagator = prop.Propagator('Sine', epoch0, dt0Func, inp)
Propagate the differential equation¶
To compute the propagation, we define the final epoch for the propagation. We choose two whole periods for the sine function, and carry out the propagation. We also specify with False flag that we do not wish to compute partial derivatives.
T = num.TwoPi / freq.value()
n = 2
epoch = epoch0 + n * T
propagator.compute(epoch, False)
To get the propagated state at a certain epoch, we can use the propagator input.
x = sca.eval(epoch)
print(x)
The expected value is indeed zero, as we propagated two periods of the sine function. The residual is the propagation error and can be controlled via the error tolerances specified above.
Plotting the solution¶
To get a full picture of the evolution of the angle, we can generate data and plot it for the range of the propagation. We will define a function to generate the solution as we will reuse this later.
def generatePlotData(pro):
grid = pro.range().createGrid(60.0)
t = [e.mjd() for e in grid]
x = [sca.eval(e) for e in grid]
return t, x
We plot it using matplotlib.
t, x = generatePlotData(propagator)
plt.plot(t, x)
plt.grid()
plt.xlabel('Time (days)')
plt.ylabel('Sine function')
plt.show()
Propagating backwards or both directions¶
We can easily propagate backwards the same way, or using a range to propagate both directions.
ran = tempo.EpochRange(epoch0 - n * T, epoch0 + n * T)
propagator.compute(ran, False)
Let's look at the complete solution via plotting again.
t, x = generatePlotData(propagator)
plt.plot(t, x)
plt.grid()
plt.xlabel('Time (days)')
plt.ylabel('Sine function')
plt.show()
Propagate including partial derivatives¶
To compute partial derivatives we simplfy have to repropagate with partials flag set to True. The same final state is computed.
propagator.compute(ran, True)
x = sca.eval(epoch)
print(x)
To get the partial derivatives of the final state, we have to create and XEpoch and use that for the evaluation.
xepoch = tempo.XEpoch(epoch)
y = sca.eval(xepoch)
print(y)
This will yield the final state with partial derivative of itself. This can be mapped, which is out of the scope of this tutorial. We can resolve the partial derivatives to the final parameter with resolve.
z = common.resolve(y)
print(z)
Visualize partial derivatives evolution¶
We can visualize partial derivatives by computing their evolution.
grid = propagator.range().createGrid(60.0)
xgrid = map(lambda e: tempo.XEpoch(e), grid)
t = [e.mjd() for e in grid]
x = [common.resolve(sca.eval(e)) for e in xgrid]
We create separate vectors for the value and each partial derivative.
value = [entry.value() for entry in x]
ddt0 = [entry.at(dt0.leaf()) for entry in x]
dx0 = [entry.at(x0.leaf()) for entry in x]
dfreq = [entry.at(freq.leaf()) for entry in x]
Finally, we plot the evolution.
plt.subplot(4, 1, 1)
plt.plot(t, value)
plt.grid()
plt.ylabel('Value')
plt.subplot(4, 1, 2)
plt.plot(t, ddt0)
plt.grid()
plt.ylabel('dt0 partial')
plt.subplot(4, 1, 3)
plt.plot(t, dx0)
plt.grid()
plt.ylabel('x0 partial')
plt.subplot(4, 1, 4)
plt.plot(t, dfreq)
plt.grid()
plt.xlabel('Time (days)')
plt.ylabel('freq partial')
plt.show()
