import astropy.units as u
import numpy as np
import gala.dynamics as gd
xy = np.arange(16).reshape(2, 8) * u.kpc
vxy = np.arange(16).reshape(2, 8) / 10. * u.kpc / u.Myr
w = gd.PhaseSpacePosition(pos=xy, vel=vxy)
fig = w.plot()