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

mw_gala = gp.MilkyWayPotential()
mw_bovy = gp.BovyMWPotential2014()
w0 = gd.PhaseSpacePosition(pos=[25., 0, 0]*u.kpc,
                           vel=[0, 0, 200.]*u.km/u.s)
orbit_gala = mw_gala.integrate_orbit(w0, dt=1., n_steps=1000)
orbit_bovy = mw_bovy.integrate_orbit(w0, dt=1., n_steps=1000)

fig, ax = plt.subplots(1, 1, figsize=(6, 6))
orbit_gala.plot(['x', 'z'], label='Gala', marker='', axes=[ax])
orbit_bovy.plot(['x', 'z'], label='Bovy2015', marker='', axes=[ax])
plt.legend(loc='best')
plt.tight_layout()