import astropy.units as u
import gala.potential as gp
import gala.dynamics as gd
from gala.units import galactic
pot = gp.MiyamotoNagaiPotential(m=2.5e11, a=6.5, b=0.26, units=galactic)
w0 = gd.PhaseSpacePosition(
    pos=[11.0, 0.0, 0.2] * u.kpc, vel=[0.0, 200, 100] * u.km / u.s
)
orbit = gp.Hamiltonian(pot).integrate_orbit(w0, dt=1.0, n_steps=1000)
_ = orbit.represent_as("cylindrical").plot(["rho", "z"])