How to use galpy - 10 common examples

def compare(ic, inttime, numsteps):
    times = numpy.linspace(0, inttime, numsteps)

    ### integrate the orbit in galpy using MWPotential2014 from galpy
    g_orb_obj = galpy.orbit.Orbit([ic[0],ic[3],ic[5],ic[1],ic[4],ic[2]])
    dt = time.time()
    g_orb_obj.integrate(times, g_pot)
    g_orb = g_orb_obj.getOrbit()
    print 'Time to integrate orbit in galpy: %.4g s' % (time.time()-dt)

    ### integrate the orbit with the galpy routine, but using Agama potential instead
    ### (much slower because of repeated transfer of control between C++ and Python
    dt = time.time()
    g_orb_obj.integrate(times[:numsteps//10], a_pot)
    a_orb = g_orb_obj.getOrbit()
    print 'Time to integrate 1/10th of the orbit in galpy using Agama potential: %.4g s' % (time.time()-dt)

    ### integrate the same orbit (note different calling conventions - cartesian coordinates as input)
    ### using both the orbit integration routine and the potential from Agama - much faster
    dt = time.time()
    times_c, c_orb_car = agama.orbit(ic=[ic[0],0,ic[1],ic[3],ic[5],ic[4]], \

### integrate the same orbit (note different calling conventions - cartesian coordinates as input)
    ### using both the orbit integration routine and the potential from Agama - much faster
    dt = time.time()
    times_c, c_orb_car = agama.orbit(ic=[ic[0],0,ic[1],ic[3],ic[5],ic[4]], \
        potential=c_pot, time=inttime, trajsize=numsteps)
    print 'Time to integrate orbit in Agama: %.4g s' % (time.time()-dt)

    ### make it compatible with galpy's convention (native output is in cartesian coordinates)
    c_orb = c_orb_car*1.0
    c_orb[:,0] = (c_orb_car[:,0]**2+c_orb_car[:,1]**2)**0.5
    c_orb[:,3] =  c_orb_car[:,2]

    ### in galpy, this is the only tool that can estimate focal distance,
    ### but it requires the orbit to be computed first
    delta = galpy.actionAngle.estimateDeltaStaeckel(g_pot, g_orb[:,0], g_orb[:,3])
    print "Focal distance estimated from the entire trajectory: Delta=%.4g" % delta

    ### plot the orbit(s) in R,z plane, along with the prolate spheroidal coordinate grid
    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

def __init__(self,*args,**kwargs):
            """
            Initialize a potential from parameters provided in an INI file
            or as named arguments to the constructor.
            Arguments are the same as for regular agama.Potential (see below);
            an extra keyword "normalize=..." has the same meaning as in Galpy:
            if True, normalize such that vc(1.,0.)=1., or,
            if given as a number, such that the force is this fraction of the force
            necessary to make vc(1.,0.)=1.

            """
            galpy.potential.Potential.__init__(self,amp=1.)
            normalize=False
            for key, value in kwargs.items():
                if key=="normalize":
                    normalize=value
                    del kwargs[key]
            self._pot = Potential(*args,**kwargs)  # regular Agama potential
            if normalize or \
                    (isinstance(normalize,(int,float)) \
                        and not isinstance(normalize,bool)):
                self.normalize(normalize)
            self.hasC= False
            self.hasC_dxdv=False
        __init__.__doc__ += Potential.__doc__

def __init__(self,*args,**kwargs):
        """
        NAME:
           __init__
        PURPOSE:
           initialize a potential from parameters provided in an INI file
           or as named arguments to the constructor (see below).
        INPUT:
           normalize - if True, normalize such that vc(1.,0.)=1., or,
           if given as a number, such that the force is this fraction of the force
           necessary to make vc(1.,0.)=1.
        HISTORY:
           2014-12-05 EV
        """
        Potential.__init__(self,amp=1.)
        normalize=False
        for key, value in kwargs.items():
            if key=="normalize":
                normalize=value
                del kwargs[key]
        self._pot = agama.Potential(*args,**kwargs)
        if normalize or \
                (isinstance(normalize,(int,float)) \
                    and not isinstance(normalize,bool)):
            self.normalize(normalize)
        self.hasC= False
        self.hasC_dxdv=False
    __init__.__doc__ += agama.Potential.__doc__

#!/usr/bin/python
### This module allows to use potentials from the AGAMA C++ library as regular galpy potentials

import math
import agama, galpy
from galpy.potential import Potential
class CPotential(galpy.potential.Potential):
    """Class that implements an interface to C++ potentials
    """
    def __init__(self,*args,**kwargs):
        """
        NAME:
           __init__
        PURPOSE:
           initialize a potential from parameters provided in an INI file
           or as named arguments to the constructor (see below).
        INPUT:
           normalize - if True, normalize such that vc(1.,0.)=1., or,
           if given as a number, such that the force is this fraction of the force
           necessary to make vc(1.,0.)=1.
        HISTORY:
           2014-12-05 EV
        """

r=[0.01, 0.0, 0.4, 0.5, 0.3, 0.0, 0.7, 1.0, 1.0, 1.0, 0.9]
    g=[0.01, 0.0, 0.85,1.0, 1.0, 0.9, 1.0, 1.0, 0.85,0.0, 0.9]
    b=[0.01, 1.0, 1.0, 1.0, 0.7, 0.0, 0.0, 0.0, 0.0, 0.0, 0.9]
    matplotlib.pyplot.register_cmap(cmap=matplotlib.colors.LinearSegmentedColormap('sauron', \
        { 'red': zip(f,r,r), 'green': zip(f,g,g), 'blue': zip(f,b,b) } ))
    matplotlib.pyplot.register_cmap(cmap=matplotlib.colors.LinearSegmentedColormap('sauron_r', \
        { 'red': zip(f,r[::-1],r[::-1]), 'green': zip(f,g[::-1],g[::-1]), 'blue': zip(f,b[::-1],b[::-1]) } ))
    del f; del r; del g; del b;  # remove the temporary variables from the module namespace

except ImportError: pass   # no matplotlib - no problem


try:
    # This wrapper class allows to use Agama potentials as regular galpy potentials
    import galpy.potential
    class AgamaPotential(galpy.potential.Potential):
        """
        Class that implements a Galpy interface to Agama potentials
        """
        def __init__(self,*args,**kwargs):
            """
            Initialize a potential from parameters provided in an INI file
            or as named arguments to the constructor.
            Arguments are the same as for regular agama.Potential (see below);
            an extra keyword "normalize=..." has the same meaning as in Galpy:
            if True, normalize such that vc(1.,0.)=1., or,
            if given as a number, such that the force is this fraction of the force
            necessary to make vc(1.,0.)=1.

            """
            galpy.potential.Potential.__init__(self,amp=1.)
            normalize=False

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]))

#!/usr/bin/python

# this demo script shows how to use Agama library in galpy:
# creating a galpy-compatible potential, integrating an orbit and using the action finder

import galpy.potential, galpy.actionAngle, galpy.orbit
import pygama as agama
import numpy, matplotlib, matplotlib.pyplot as plt, time
matplotlib.rcParams['legend.frameon']=False

### set up galpy potential
g_pot = galpy.potential.MWPotential2014

### set up equivalent potential from the Agama library
agama.setUnits( mass=1., length=8., velocity=220)
a_pot = agama.AgamaPotential("../data/MWPotential2014galpy.ini")
c_pot = a_pot._pot   # the instance of raw Agama potential

### initialization of the action finder needs to be done once for the given potential
dt = time.time()
c_actfinder = agama.ActionFinder(c_pot, interp=False)
print 'Time to set up action finder: %.4g s' % (time.time()-dt)
### we have a faster but less accurate "interpolated action finder", which takes a bit longer to initialize
dt = time.time()
i_actfinder = agama.ActionFinder(c_pot, interp=True)
print 'Time to set up interpolated action finder: %.4g s' % (time.time()-dt)

### conversion from prolate spheroidal to cylindrical coords

#First evaluate the dust map
        ds= numpy.linspace(range[0],range[1],nds)
        if distmod:
            adust= self(l,b,10.**(ds/5.-2.))
        else:
            adust= self(l,b,ds)
        #Add labels
        if distmod:
            kwargs['xlabel']= r'$\mathrm{Distance\ modulus}$'
        else:
            kwargs['xlabel']= r'$D\,(\mathrm{kpc})$'
        if not self._filter is None:
            kwargs['ylabel']= r'$A_{%s}\,(\mathrm{mag})$' % (self._filter.split(' ')[-1])
        else:
            kwargs['ylabel']= r'$E(B-V)\,(\mathrm{mag})$'
        return bovy_plot.bovy_plot(ds,adust,*args,**kwargs)

def galpy_nfw_orbit():
    # Setting up the potential
    nfw = NFWPotential(conc=C, mvir=M_200, H=70.0, wrtcrit=True, overdens=200)
    nfw.turn_physical_on()
    vxvv = [
        8.0 * units.kpc,
        0.0 * units.km / units.s,
        240.0 * units.km / units.s,
        0.0 * units.pc,
        5.0 * units.km / units.s,
    ]

    # Calculating the orbit
    ts = np.linspace(0.0, 0.58, 1000) * units.Gyr
    o = Orbit(vxvv=vxvv)
    o.integrate(ts, nfw, method="odeint")

    return o