(1566) Icarus propagation test

import grss
prop = grss.prop
fit = grss.fit

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

body_id = '1566'
init_sol, init_cov, nongrav_info = fit.get_sbdb_info(body_id)

e = init_sol['e']
q = init_sol['q']
tp = init_sol['tp']
om = init_sol['om']
w = init_sol['w']
i = init_sol['i']

cometary_elements = [e, q, tp, om, w, i]
ng_params = prop.NongravParameters()
ng_params.a1 = init_sol.get('a1', 0.0)
ng_params.a2 = init_sol.get('a2', 0.0)
ng_params.a3 = init_sol.get('a3', 0.0)

de_kernel = 440
de_kernel_path = grss.utils.default_kernel_path
t0 = init_sol['t']
daysInYear = 365.25
numYears = 10
numDays = numYears * daysInYear
tf = t0 + numDays
prop_sim = prop.PropSimulation("(1566) Icarus propagation test", t0, de_kernel, de_kernel_path)

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)

icarus = prop.IntegBody("(1566) Icarus", t0, 0.0, 0.0, cometary_elements, ng_params)
prop_sim.add_integ_body(icarus)

prop_sim.integrate()

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=0.2)
axs.set_xlabel("x [AU]")
axs.set_ylabel("y [AU]")
axs.set_aspect('equal')
plt.show(block=False)
plt.close()