import astropy.units as u
import matplotlib.pyplot as pl
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, r_h=0.1*u.kpc,
                              q1=1., q2=0.8, q3=0.6, units=galactic)
w0 = gd.PhaseSpacePosition(pos=[5.5,0.,0.]*u.kpc,
                           vel=[0.,140.,25]*u.km/u.s)
lyap,orbit = gd.fast_lyapunov_max(w0, pot, dt=1., n_steps=200000)
fig = orbit[:,0].plot(marker=',', linestyle='none', alpha=0.1)