import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
import gala.potential as gp
from gala.units import galactic, solarsystem

p1 = gp.KeplerPotential(m=1*u.Msun, origin=[1, 0, 0]*u.au,
                        units=solarsystem)
p2 = gp.KeplerPotential(m=0.5*u.Msun, origin=[-2, 0, 0]*u.au,
                        units=solarsystem)

pot = gp.CCompositePotential(p1=p1, p2=p2)
fig, ax = plt.subplots(1, 1, figsize=(5, 5))
grid = np.linspace(-5, 5, 100)
pot.plot_contours(grid=(grid, grid, 0.), ax=ax) # doctest: +SKIP
ax.set_xlabel("$x$ [kpc]") # doctest: +SKIP
ax.set_ylabel("$y$ [kpc]") # doctest: +SKIP
fig.tight_layout()