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()
w0 = gd.PhaseSpacePosition(pos=[-8.1, 0, 0.02] * u.kpc,
                           vel=[13, 245, 8.] * u.km/u.s)
orbit = mw.integrate_orbit(w0, dt=1*u.Myr, t1=0, t2=2*u.Gyr)

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