Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
plt.figure(figsize=(12,8))
plt.axes([0.04, 0.54, 0.45, 0.45])
plotCoords(delta, 1.5)
plt.plot(g_orb[:,0],g_orb[:,3], 'b', label='galpy') # R,z
plt.plot(c_orb[:,0],c_orb[:,3], 'g', label='Agama')
plt.plot(a_orb[:,0],a_orb[:,3], 'r', label='galpy using Agama potential')
plt.xlabel("R/8kpc")
plt.ylabel("z/8kpc")
plt.xlim(0, 1.2)
plt.ylim(-1,1)
plt.legend(loc='lower left')
### create galpy action/angle finder for the given value of Delta
### note: using c=False in the routine below is much slower but apparently more accurate,
### comparable to the Agama for the same value of delta
g_actfinder = galpy.actionAngle.actionAngleStaeckel(pot=g_pot, delta=delta, c=True)
### find the actions for each point of the orbit using galpy action finder
dt = time.time()
g_act = g_actfinder(g_orb[:,0],g_orb[:,1],g_orb[:,2],g_orb[:,3],g_orb[:,4],fixed_quad=True)
print 'Time to compute actions in galpy: %.4g s' % (time.time()-dt)
print 'Jr = %.6g +- %.4g, Jz = %.6g +- %.4g' % \
(numpy.mean(g_act[0]), numpy.std(g_act[0]), numpy.mean(g_act[2]), numpy.std(g_act[2]))
### use the Agama action routine for the same value of Delta as in galpy -
### the result is almost identical but computed much faster
dt = time.time()
c_act = agama.actions(point=c_orb_car, potential=c_pot, fd=delta) # explicitly specify focal distance
print 'Time to compute actions in Agama using Galpy-estimated focal distance: %.4g s' % (time.time()-dt)
print 'Jr = %.6g +- %.4g, Jz = %.6g +- %.4g' % \
(numpy.mean(c_act[:,0]), numpy.std(c_act[:,0]), numpy.mean(c_act[:,1]), numpy.std(c_act[:,1]))