Possible issues with LogarithmicHalo density and Hessian

[nstarman]

See GalacticDynamics/galax#296 (comment).
I'm seeing differences in the galax results (from auto-diff) to the gala numbers (coded here).
This is either a problem on the galax side resulting from machine precision in the auto-diffed function that is then multiplied by a large scale factor from a unit conversion, or there's a mistake in gala! Either way, we should fix this.

[adrn]

Did you ever look into this more? I generated the code you linked from sympy, and just checked again that Gala agrees with sympy. So the problem might be on the Galax side!

import gala.potential as gp
import numpy as np
import sympy as sy

vc = 1.0
rh = 0.1
q1 = 1.0
q2 = 0.95
q3 = 0.92
pot = gp.LogarithmicPotential(vc, rh, q1, q2, q3)

x, y, z = sy.symbols("x y z")

pot_expr = 0.5 * vc**2 * sy.log(rh**2 + x**2 / q1**2 + y**2 / q2**2 + z**2 / q3**2)
sy_pot = sy.lambdify((x, y, z), pot_expr)
sy_dens = sy.lambdify(
    (x, y, z), sy.trace(sy.hessian(pot_expr, (x, y, z))) / (4 * np.pi)
)

rng = np.random.default_rng(seed=42)
trial_x = rng.uniform(-1, 1, (3, 128))

assert np.allclose(sy_pot(*trial_x), pot.energy(trial_x))
assert np.allclose(sy_dens(*trial_x), pot.density(trial_x))

[adrn]

@nstarman did we ever resolve this?

[nstarman]

I think this is resolved by GalacticDynamics/galax@2f11f65, but worth double checking.