import numpy as np
import gala.dynamics as gd
import gala.potential as gp
from gala.units import galactic

disk = gp.MiyamotoNagaiPotential(m=1E11, a=6.5, b=0.27, units=galactic)
bulge = gp.HernquistPotential(m=3E10, c=0.7, units=galactic)
pot = gp.CompositePotential(disk=disk, bulge=bulge)

grid = np.linspace(-3.,3.,100)
fig = pot.plot_contours(grid=(grid,0,grid))