potential = scf.SCFPotential(Snlm=S, Tnlm=T, m=M, r_s=a) # M=a=1

# we need an array of positions with shape (3,n_samples) for SCFPotential
xyz = np.vstack((x.ravel(),y.ravel(),z.ravel()))
scf_dens = potential.density(xyz)

# log-spaced contour levels
true_dens = flattened_hernquist_density(x, y, z, M, a, q)
levels = np.logspace(np.log10(true_dens.min()), np.log10(true_dens.max()), 16)

plt.figure(figsize=(6,6))

plt.contourf(x, z, true_dens, cmap='magma',
             levels=levels, locator=ticker.LogLocator())
plt.contour(x, z, scf_dens.reshape(x.shape), colors='w',
            levels=levels, locator=ticker.LogLocator())

plt.title("Isodensity")
plt.xlabel("$x$", fontsize=22)
plt.ylabel("$z$", fontsize=22)
plt.tight_layout()