GalacticDynamics-Oxford / Agama

Action-based galaxy modeling framework
Other
72 stars 35 forks source link

central density of King profile #32

Open syrte opened 1 year ago

syrte commented 1 year ago

I found a strange behavior with King profile in agama:

import agama
den = agama.Density(type='King', W0=10)
den.density([0., 0, 0]), den.density([1e-100, 0, 0])
# output
# (inf, 0.007955322088520804)

I am sorry if it has already been resolved in recent updates. My agama version is '1.0 compiled on Jan 12 2022'. Anyway, I would be fine with den.density([1e-100, 0, 0]) as a proxy for now :)

BTW, I wanted the density in the center because I am calculating the amplitude of the distribution function (Gieles & Zocchi 2015, eq. 1 and 30).

eugvas commented 1 year ago

yes, that is a consequence of the implementation: the King density or potential (these are not identical) are internally represented respectively as DensitySphericalHarmonic or Multipole potential, and the values are stored on a logarithmic grid in radius (thus excluding zero). As the interpolation spline constructed on this log-grid has a tiny but nonzero slope at the innermost point, it continues to rise to infinity, even though the value at 0+epsilon is finite and could be used as a proxy for the central value. I'll see if this can be rectified by some reasonably simple modification of the interpolation scheme. Btw, the DF of the King model could in principle be constructed using the QuasiSpherical distribution function from the corresponding density and potential, but the result is noisy and often unusable, again as a consequence of interpolation and numerical inaccuracies (the DF is a 1.5-order derivative of the density!). So I would indeed recommend to take the central density and use the analytic expression for the DF.

syrte commented 1 year ago

Got it. Many thanks for your super quick reply! BTW, I noticed that Agama is ~60x faster than limepy when initializing a King profile 👍 This is critical in some applications.