Error in ActionFinder initialization: Interpolator: energy or angular momentum of a circular orbit are not monotonic

[yunpeng-jia]

Dear Dr. Vasiliev, I got an error when using ActionFinder(python) to calculate actions. The potential contains two components, a halo(analytic) and a disk(from Gadget snapshots files). However, the disk contains 1000 large mass particles that distributed non-axisymmetric in the disk(the disk contains 5 million disk particles). Using type='Multipole' or 'CylSpline' gives similar errors for some of snapshots files. By the way, agama.ActionFinder works well for most of the snapshots files. Really appreciate your time here. Regards, Yunpeng Jia

Here are the details: pot_disk_nbody = agama.Potential(type='Multipole', particles=(disk_cor, disk_mass), symmetry='Axisymmetric') pot_halo_Hernquist = agama.Potential(type='Dehnen', mass=1.0e12, gamma=1, axisRatioY=1, axisRatioZ=1, scaleRadius=47) pot_nbody = agama.Potential(pot_disk_nbody, pot_halo_Hernquist) af_fudge = agama.ActionFinder(pot_nbody,interp=True) `

_RuntimeError: Error in ActionFinder initialization: Interpolator: energy or angular momentum of a circular orbit are not monotonic with radius at R=0.0341387 Stack trace (18 functions, starting from the most recent one): /home/jia/anaconda3/lib/python3.11/site-packages/agama/agama.so(+0x122a0c) [0x7f55803e6a0c] /home/jia/anaconda3/lib/python3.11/site-packages/agama/agama.so(actions::ActionFinderAxisymFudge::ActionFinderAxisymFudge(std::shared_ptr const&, bool)+0xff) [0x7f55805433ff] /home/jia/anaconda3/lib/python3.11/site-packages/agama/agama.so(actions::createActionFinder(std::shared_ptr const&, bool)+0x2cc) [0x7f5580530ffc] /home/jia/anaconda3/lib/python3.11/site-packages/agama/agama.so(pygama::ActionFinder_init(pygama::ActionFinderObject, _object, _object*)+0xbc) [0x7f558040c75c] python(_PyObject_MakeTpCall+0x22b) [0x502d2b] python(_PyEval_EvalFrameDefault+0x755) [0x50f025] python() [0x5c82ce] python(PyEval_EvalCode+0x9f) [0x5c79cf] python() [0x5e8807] python() [0x5e4e40] python() [0x5f9132] python(_PyRun_SimpleFileObject+0x19f) [0x5f871f] python(_PyRun_AnyFileObject+0x43) [0x5f8473] python(Py_RunMain+0x2ee) [0x5f2fee] python(Py_BytesMain+0x39) [0x5b6e19] /lib/x86_64-linux-gnu/libc.so.6(__libc_startmain+0xf3) [0x7f55926dd083] python() [0x5b6c6f]

[eugvas]

In the latest commit I added a fix to a long-standing but rare case when this problem may occur with Multipole potentials created from analytic density profiles, but I don't think this would resolve your issue (but please check anyway). If it doesn't help, you may need to manually adjust the grid radii (rmin/rmax) when constructing the potential from an N-body snapshot. Normally the automatic grid selection procedure works reasonably well, but there may be exceptions. If nothing else helps, can you send me the complete example (including the snapshot files shared somewhere) so that I can investigate it further.

[yunpeng-jia]

Thank you for your prompt reply. However, the errors persist even after I manually set the grid radii (rmin=, rmax=). I have sent you an email with an example that includes a snapshot file and the Python code.

[eugvas]

Thanks for sending the snapshot; I've replied by email, but not sure if you've received it, so repeating here. The problem is that the snapshot is not properly centered, see the attached image - the offset is less than 1 kpc (1 grid cell), but sufficient to derail the method. It may be because of substructures that the center of mass does not coincide with the center of density; to alleviate the problem, I usually do something like this (take a median rather than a mean, and use progressively reducing cutoff in radius):

def recenter(pos, radius=np.inf):
    select = np.sum(pos**2, axis=1) < radius**2
    x = np.median(pos[select], axis=0)
    pos -= x

and then

recenter(pos)       # entire snapshot
recenter(pos, 20.0) # only particles within 20 kpc
recenter(pos, 5.0)  # gradually reducing the max radius, may be repeat a couple of times

Note that for CylSpline the plane of the disk should be aligned with x-y (and also for Multipole if you use axisymmetric geometry, which is needed for the action finder to work anyway; but multipole is not recommended for highly flattened systems, unless you go to high order in lmax >~12-20). This requirement seems to be satisfied in your case.