import astropy.coordinates as coord
import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
import gala.potential as gp
import gala.dynamics as gd
from gala.units import galactic

pot = gp.LogarithmicPotential(v_c=150*u.km/u.s, q1=1., q2=1., q3=0.9, r_h=0,
                              units=galactic)
w0 = gd.PhaseSpacePosition(pos=[8, 0, 0.]*u.kpc,
                           vel=[75, 150, 50.]*u.km/u.s)

w = gp.Hamiltonian(pot).integrate_orbit(w0, dt=0.5, n_steps=10000)
cyl = w.represent_as('cylindrical')
cyl.plot(['rho', 'z'], linestyle='-')