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)

pl.figure()
pl.loglog(orbit.t[11::10], lyap, marker='')
pl.xlabel("Time [{}]".format(orbit.t.unit))
pl.ylabel(r"$\lambda_{{\rm max}}$ [{}]".format(lyap.unit))
pl.tight_layout()