Angle-Action Estimation by Locally Fitting Stackel Potentials

written in C++ by Jason Sanders (jason.sanders@physics.ox.ac.uk).

Function

Given a list of 6D phase-space points this code estimates the corresponding actions and angles in a given general axisymmetric potential. For each phase-space point the code

1. Integrates the orbit in the true potential to find the region the orbit probes,
2. fits a Stäckel potential to this region,
3. estimates the angles and actions in the true potential as those calculated in the fit potential.

For more information see the associated [ArXiv link (paper)] or the pdf included in the download (doc/AnglesActionsInAxisymmetricPotential-JasonSanders.pdf).

Download

All the relevant files are included in aalfs.tar.gz.

Compiling

The code should simply compile with the makefile in the download. However, you may have to change the compiler used in the makefile.

Use

The code compiles to one executable xv2aa. This executable takes two command-line arguments - the first is a list of input 6D phase-space points $(x,y,z,v_{x},v_{y},v_{z})$ and the second specifies an output file. The list of actions and angles $(J_{R},L_{z},J_{z},\theta _{R},\theta _{\phi },\theta _{z})$ is given in the first six columns of the output file. For details of the other columns see src/xv2aa.cpp.

Testing

There is a test file (MNJR0.001Jz0.001Lz1.0.txt) from which the file TEST_OUT.txt should be reproduced by the file OUTPUT.txt with the command:

./xv2aa MNJR0.001Jz0.001Lz1.0.txt OUTPUT.txt

which takes approx. 3.8s on 2.5GHz Intel Core i5 and results in errors in the actions of 0.1%.

Referencing

If you use this code please reference the paper for which it was developed: [ArXivReference] or [ADSReference].

Other

Potentials

This distributed version of the code uses a simple potential (Miyamoto-Nagai) as a test potential (defined in PotMN.cpp). The user will probably want to specify their own potential which they can do by altering the functions in this file. The user must specify the potential and its gradients as a function of the polar coordinates $(R,z)$. The code is designed to work with a realistic Galactic potential such as those defined by Dehnen & Binney 1998 or McMillan 2011. To use these potentials you need to use Walter Dehnen's GalPot code which is available from Walter.

Double Integration

The code involves one double integration. The implementation here is a simple double Gauss Legendre quadrature (DoubleGaussLegendreQuad) but for a better routine one can use the cubature library. The relevant code to alter is in fitpot.cpp.