Ctrl+K

(65803) Didymos propagation test

(65803) Didymos propagation test#

Let’s start by importing the necessary libraries#

[1]:

from grss import prop, utils

[2]:

import numpy as np
import matplotlib.pyplot as plt
np.set_printoptions(precision=40, linewidth=np.inf)

We’ll then define the initial barycentric equatorial Cartesian state of the asteroid plus any nongravitational accelerations acting on it.#

[3]:

t0 = 59956.5
pos = [-4.739057923859238E-01, 1.095212783390737E+00, 5.270678006987460E-01]
vel = [-1.653652752454120E-02, -8.564528317342794E-04, 6.471815241318629E-04]
ng_params = prop.NongravParameters()
ng_params.a1 = 0.0
ng_params.a2 = -1.042309691002E-14
ng_params.a3 = 0.0
ng_params.alpha = 1.0
ng_params.k = 0.0
ng_params.m = 2.0
ng_params.n = 0.0
ng_params.r0_au = 1.0
didymos = prop.IntegBody("(65803) Didymos", t0, 0.0, 0.0, pos, vel, ng_params)

Let’s choose the planetary kernels, final time, and initialize the simulation.#

[4]:

de_kernel = 440
de_kernel_path = utils.default_kernel_path
tf = t0 + 3000
prop_sim = prop.PropSimulation("(65803) Didymos propagation test", t0, de_kernel, de_kernel_path)

We need to define the times at which we want to retrieve the state of the asteroid, plus some choices for the evaluated states:#

  • Whether we want to retrieve the geometric or apparent state of the asteroid.

  • Whether the evaluation times are in the TDB or UTC time scale.

  • Whether we want to use a converged light time solution if the apparent state is chosen.

[5]:

t_temp = t0
t_eval = []
while t_temp <= tf:
    t_eval.append(t_temp)
    t_temp += 1.0

eval_apparent_state = False
t_eval_utc = False
converged_light_time = False
prop_sim.set_integration_parameters(tf, t_eval, t_eval_utc, eval_apparent_state, converged_light_time)

Finally, we add the asteroid to the simulation and propagate it.#

[6]:

prop_sim.add_integ_body(didymos)
prop_sim.integrate()

Plot the results#

[7]:

arr = np.array(prop_sim.xIntegEval)
# plot xy and xz projections of the orbit in a 2x1 subplot
plt.figure(figsize=(6, 6), dpi=150)
axs = plt.gca()
prop.plot_solar_system(axs, xy_plane=True, alpha=0.5)
axs.plot(arr[:,0], arr[:,1], '-', lw=1.0, label="Didymos")
axs.set_xlabel("x [AU]")
axs.set_ylabel("y [AU]")
axs.set_aspect('equal')
plt.legend()
plt.show(block=False)
plt.close()
../../../_images/tests_python_prop_didymos_13_0.png

Make sure the final state is correct#

[8]:

jpl = np.array([ 1.1601957603540887e+11,  9.1763914866120560e+10,  3.4431641402258896e+10, -2.1333356367475884e+04,  2.4445528524826328e+04,  1.2446287491647476e+04])
grss = np.array(prop_sim.xInteg) * prop_sim.consts.du2m
grss[3:6] /= 86400.0
pos_diff = np.linalg.norm((jpl-grss)[:3])
assert pos_diff < 100.0, f"Position difference {pos_diff} m is too large"

[ ]:
Contents