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

pot = gp.NFWPotential(m=1e11 * u.Msun, r_s=20.0 * u.kpc, units=galactic)
r = np.logspace(np.log10(20.0 / 100.0), np.log10(20 * 100.0), 100) * u.kpc
m_profile = pot.mass_enclosed(r=r)

plt.figure()
plt.loglog(r, m_profile, marker="")  # doctest: +SKIP
plt.xlabel("$r$ [{}]".format(r.unit.to_string(format="latex")))
plt.ylabel("$M(<r)$ [{}]".format(m_profile.unit.to_string(format="latex")))
plt.tight_layout()