import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
import gala.dynamics as gd
import gala.potential as gp

mw = gp.MilkyWayPotential(version="latest")
w0 = gd.PhaseSpacePosition(
    pos=[-8.1, 0, 0.02] * u.kpc,
    vel=[13, 245, 8.0] * u.km / u.s,
)
orbit = mw.integrate_orbit(w0, dt=1 * u.Myr, t1=0, t2=2 * u.Gyr)

orbit.plot(["x", "y"])