Code to replicate galpy issue #361 where for certain orbits in a triaxial potential the ActionAngleAdiabatic form returns 0 for the radial action and an anomalously high value for the quasi-isothermal DF.

galpy_aA_adiabatic_triaxial_issue.py

import copy, tqdm, pdb
import numpy as np
from matplotlib import pyplot as plt
from astropy import units as apu
from galpy import orbit
from galpy import potential
from galpy import df
from galpy.actionAngle import actionAngleAdiabatic
from galpy.util import bovy_conversion as gpconv

### Setup

# Times for halo evolution
t_evolve = 10 # Gyr
tform = -9 # Gyr ago
tsteady = 8 # Gyr after tform
times = -np.array([0,t_evolve]) * apu.Gyr

# Get the parameters of MWPotential2014
mwbulge = copy.deepcopy(potential.MWPotential2014[0])
mwdisk = copy.deepcopy(potential.MWPotential2014[1])
mwhalo = copy.deepcopy(potential.MWPotential2014[2])

# Make the triaxial halo
trihalo = potential.TriaxialNFWPotential(amp=mwhalo._amp, a=mwhalo.a, 
    b=2.0, pa=np.linspace(0, np.pi/2, 5)[3], c=1.0)

# Wrap the old halo in a decaying DSW
mwhalo_decay_dsw = potential.DehnenSmoothWrapperPotential(pot=mwhalo, 
    tform=tform*apu.Gyr, tsteady=tsteady*apu.Gyr, decay=True)
    
# Wrap the triaxial halo in a growing DSW:
trihalo_dsw = potential.DehnenSmoothWrapperPotential(pot=trihalo,
    tform=tform*apu.Gyr, tsteady=tsteady*apu.Gyr)

# Make the triaxial potential
tripot = [mwbulge, mwdisk, mwhalo_decay_dsw, trihalo_dsw]
potential.turn_physical_off(tripot)

# Make action angle object
qdf_aA= actionAngleAdiabatic(pot=potential.MWPotential2014, c=True)
qdf = df.quasiisothermaldf( hr= 2*apu.kpc,
                            sr= 46*(apu.km/apu.s),
                            sz= 28*(apu.km/apu.s),
                            hsr= 9.8*(apu.kpc),
                            hsz= 7.6*(apu.kpc),
                            pot= potential.MWPotential2014, 
                            aA= qdf_aA)
    
### Integrate orbits    
                        
dvT = 10.
dvR = 10.
vR_low = -40
vR_hi = 40
vT_low = 180
vT_hi = 260

vR_range = np.arange( vR_low, vR_hi, dvR )
vT_range = np.arange( vT_low, vT_hi, dvT )
dfp = np.zeros((len(vR_range),len(vT_range)))
actp = np.zeros((3,len(vR_range),len(vT_range)))

## Calculate the perturbed DF and actions
for i in tqdm.tqdm( range( len(vR_range) ) ):
    for j in range( len(vT_range) ):
        o = orbit.Orbit(vxvv=[8.*apu.kpc,
                              vR_range[i]*apu.km/apu.s,
                              vT_range[j]*apu.km/apu.s,
                              0.*apu.kpc,
                              0.*apu.km/apu.s,
                              0.*apu.radian
                              ]
                        )
        o.integrate(times, tripot)
        
        # Actions
        actp[0,i,j],actp[1,i,j],actp[2,i,j] = qdf_aA( o(times[-1]) )
        # DF
        dfp[i,j] = qdf( o(times[-1]) )
    ###j
###i

### Plot

# Plot perturbed DF
fig = plt.figure()
ax = fig.add_subplot(111)
img_arr = np.rot90( dfp )
img = ax.imshow(img_arr, interpolation='nearest',
    extent=[vR_low, vR_hi, vT_low, vT_hi],
    cmap='Blues')
cbar = plt.colorbar(img, ax=ax)
cbar.set_label(r'$DF$')
ax.set_ylabel(r'$V_{\phi}$ [km/s]')
ax.set_xlabel(r'$V_{R}$ [km/s]')
plt.savefig('test_DF.pdf')
plt.close('all')

# Plot perturbed radial action
fig = plt.figure()
ax = fig.add_subplot(111)
img_arr = np.rot90( actp[0] )
img = ax.imshow(img_arr, interpolation='nearest',
    extent=[vR_low, vR_hi, vT_low, vT_hi],
    cmap='Blues')
cbar = plt.colorbar(img, ax=ax)
cbar.set_label(r'$J_{R}$')
ax.set_ylabel(r'$V_{\phi}$ [km/s]')
ax.set_xlabel(r'$V_{R}$ [km/s]')
plt.savefig('test_JR.pdf')