# How to use galpy - 10 common examples

## To help you get started, we’ve selected a few galpy examples, based on popular ways it is used in public projects.

Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.

GalacticDynamics-Oxford / Agama / pytests / example_galpy.py View on Github
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]], \
GalacticDynamics-Oxford / Agama / pytests / example_galpy.py View on Github
### 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
GalacticDynamics-Oxford / Agama / pytests / pygama.py View on Github
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__
GalacticDynamics-Oxford / Agama / pytests / galpy_agama.py View on Github
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__
GalacticDynamics-Oxford / Agama / pytests / galpy_agama.py View on Github
#!/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
"""
GalacticDynamics-Oxford / Agama / pytests / pygama.py View on Github
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
GalacticDynamics-Oxford / Agama / pytests / example_galpy.py View on Github
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()
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]))
GalacticDynamics-Oxford / Agama / pytests / example_galpy.py View on Github
#!/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
jobovy / mwdust / mwdust / DustMap3D.py View on Github
#First evaluate the dust map
ds= numpy.linspace(range[0],range[1],nds)
if distmod:
else:
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)
SWIFTSIM / swiftsim / examples / GravityTests / NFW_Halo / makePlots.py View on Github
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

## galpy

Galactic Dynamics in python

BSD-3-Clause