import numpy as np
import matplotlib.pyplot as plt
from astropy.time import Time
from mskpy import KeplerState
import cometsuite as cs
date = Time(2450643.5417, format="jd")
rc = [4.36955224e7, -1.64073549e8, -2.73990869e7]
vc = [26.40398803, -21.02309023, -1.63325062]
re = [5.61856527e7, -1.41307139e8, -1.23261993e3]
ve = [2.71884498e1, 1.09043893e1, -6.22859821e-04]
comet = KeplerState(rc, vc, date, name="Encke")
earth = KeplerState(re, ve, date, name="Earth")
beta = [0.001, 0.002, 0.004, 0.006, 0.008, 0.01, 0.1]
ndays = 200
steps = 101
sim = cs.quick_syndynes(comet,
                        date,
                        beta=beta,
                        ndays=ndays,
                        steps=steps,
                        observer=earth)
ax = plt.gca()
plt.setp(ax,
         xlabel='Position angle',
         ylabel=r'$\rho$ (arcsec)',
         rmax=(2 * 3600))  # 4 degree FOV