import astropy.units as u
import numpy as np
import gala.dynamics as gd

np.random.seed(42)
x = np.random.uniform(-10,10,size=(3,128))
v = np.random.uniform(-200,200,size=(3,128))
w = gd.PhaseSpacePosition(pos=x*u.kpc,
                          vel=v*u.km/u.s)
fig = w.plot()