Introduction
The most basic features of galpy are its ability to display rotation
curves and perform orbit integration for arbitrary combinations of
potentials. This section introduce the most basic features of
galpy.potential and galpy.orbit.
Rotation curves
The following code example shows how to initialize a Miyamoto-Nagai disk potential and plot its rotation curve
>>> from galpy.potential import MiyamotoNagaiPotential
>>> mp= MiyamotoNagaiPotential(a=0.5,b=0.0375,normalize=1.)
>>> mp.plotRotcurve(Rrange=[0.01,10.],grid=1001)
The normalize=1. option normalizes the potential such that the
radial force is a fraction normalize=1. of the radial force
necessary to make the circular velocity 1 at R=1.
Similarly we can initialize other potentials and plot the combined
rotation curve
>>> from galpy.potential import NFWPotential, HernquistPotential
>>> mp= MiyamotoNagaiPotential(a=0.5,b=0.0375,normalize=.6)
>>> np= NFWPotential(a=4.5,normalize=.35)
>>> hp= HernquistPotential(a=0.6/8,normalize=0.05)
>>> from galpy.potential import plotRotcurve
>>> plotRotcurve([hp,mp,np],Rrange=[0.01,10.],grid=1001,yrange=[0.,1.2])
Note that the normalize values add up to 1. such that the circular
velocity will be 1 at R=1. The resulting rotation curve is
approximately flat. To show the rotation curves of the three
components do
>>> mp.plotRotcurve(Rrange=[0.01,10.],grid=1001,overplot=True)
>>> hp.plotRotcurve(Rrange=[0.01,10.],grid=1001,overplot=True)
>>> np.plotRotcurve(Rrange=[0.01,10.],grid=1001,overplot=True)
You’ll see the following
As a shortcut the [hp,mp,np] Milky-Way-like potential is defined as
>>> from galpy.potential import MWPotential
This is not the recommended Milky-Way-like potential in
galpy. The (currently) recommended Milky-Way-like potential is
MWPotential2014:
>>> from galpy.potential import MWPotential2014
MWPotential2014 has a more realistic bulge model and is actually
fit to various dynamical constraints on the Milky Way (see
here and the galpy paper).
Units in galpy
Above we normalized the potentials such that they give a circular
velocity of 1 at R=1. These are the standard galpy units (sometimes
referred to as natural units in the documentation). galpy will work
most robustly when using these natural units. When using galpy to
model a real galaxy with, say, a circular velocity of 220 km/s at R=8
kpc, all of the velocities should be scaled as v= V/[220 km/s] and all
of the positions should be scaled as x = X/[8 kpc] when using galpy’s
natural units.
For convenience, a utility module bovy_conversion is included in
galpy that helps in converting between physical units and natural
units for various quantities. For example, in natural units the
orbital time of a circular orbit at R=1 is \(2\pi\); in
physical units this corresponds to
>>> from galpy.util import bovy_conversion
>>> print 2.*numpy.pi*bovy_conversion.time_in_Gyr(220.,8.)
0.223405444283
or about 223 Myr. We can also express forces in various physical
units. For example, for the Milky-Way-like potential defined in galpy,
we have that the vertical force at 1.1 kpc is
>>> from galpy.potential import MWPotential2014, evaluatezforces
>>> -evaluatezforces(1.,1.1/8.,MWPotential2014)*bovy_conversion.force_in_pcMyr2(220.,8.)
2.0259181908629933
which we can also express as an equivalent surface-density by dividing
by \(2\pi G\)
>>> -evaluatezforces(1.,1.1/8.,MWPotential2014)*bovy_conversion.force_in_2piGmsolpc2(220.,8.)
71.658016957792356
Because the vertical force at the solar circle in the Milky Way at 1.1
kpc above the plane is approximately \(70\,(2\pi G\,
M_\odot\,\mathrm{pc}^{-2})\) (e.g., 2013arXiv1309.0809B), this shows
that our Milky-Way-like potential has a realistic disk (at least in
this respect).
bovy_conversion further has functions to convert densities,
masses, surface densities, and frequencies to physical units (actions
are considered to be too obvious to be included); see here for a full list. As a final example, the local dark
matter density in the Milky-Way-like potential is given by
>>> MWPotential2014[2].dens(1.,0.)*bovy_conversion.dens_in_msolpc3(220.,8.)
0.0075419566970079373
or
>>> MWPotential2014[2].dens(1.,0.)*bovy_conversion.dens_in_gevcc(220.,8.)
0.28643101789044584
or about \(0.0075\,M_\odot\,\mathrm{pc}^{-3} \approx
0.3\,\mathrm{GeV\,cm}^{-3}\), in line with current measurements (e.g.,
2012ApJ...756...89B).
When galpy Orbits are initialized using a distance scale ro=
and a velocity scale vo= output quantities returned and plotted in
physical coordinates. Specifically, positions are are returned in
units of kpc, velocities in km/s, energies and the Jacobi integral in
(km/s)^2, the angular momentum o.L() and actions in km/s kpc,
frequencies in 1/Gyr, and times and periods in Gyr.
Orbit integration
We can also integrate orbits in all galpy potentials. Going back to a
simple Miyamoto-Nagai potential, we initialize an orbit as follows
>>> from galpy.orbit import Orbit
>>> mp= MiyamotoNagaiPotential(a=0.5,b=0.0375,amp=1.,normalize=1.)
>>> o= Orbit(vxvv=[1.,0.1,1.1,0.,0.1])
Since we gave Orbit() a five-dimensional initial condition
[R,vR,vT,z,vz], we assume we are dealing with a three-dimensional
axisymmetric potential in which we do not wish to track the
azimuth. We then integrate the orbit for a set of times ts
>>> import numpy
>>> ts= numpy.linspace(0,100,10000)
>>> o.integrate(ts,mp,nethod='odeint')
Now we plot the resulting orbit as
Which gives
The integrator used is not symplectic, so the energy error grows with time, but is small nonetheless
When we use a symplectic leapfrog integrator, we see that the energy
error remains constant
>>> o.integrate(ts,mp,method='leapfrog')
>>> o.plotE(xlabel=r'$t$',ylabel=r'$E(t)/E(0)$')
Because stars have typically only orbited the center of their galaxy
tens of times, using symplectic integrators is mostly unnecessary
(compared to planetary systems which orbits millions or billions of
times). galpy contains fast integrators written in C, which can be
accessed through the method= keyword (e.g.,
integrate(...,method='dopr54_c') is a fast high-order
Dormand-Prince method).
When we integrate for much longer we see how the orbit fills up a
torus (this could take a minute)
>>> ts= numpy.linspace(0,1000,10000)
>>> o.integrate(ts,mp,method='odeint')
>>> o.plot()
As before, we can also integrate orbits in combinations of potentials. Assuming mp, np, and hp were defined as above, we can
>>> ts= numpy.linspace(0,100,10000)
>>> o.integrate(ts,[mp,hp,np])
>>> o.plot()
Energy is again approximately conserved
>>> o.plotE(xlabel=r'$t$',ylabel=r'$E(t)/E(0)$')
Escape velocity curves
Just like we can plot the rotation curve for a potential or a
combination of potentials, we can plot the escape velocity curve. For
example, the escape velocity curve for the Miyamoto-Nagai disk defined
above
>>> mp.plotEscapecurve(Rrange=[0.01,10.],grid=1001)
or of the combination of potentials defined above
>>> from galpy.potential import plotEscapecurve
>>> plotEscapecurve([mp,hp,np],Rrange=[0.01,10.],grid=1001)
For the Milky-Way-like potential MWPotential2014, the
escape-velocity curve is
>>> plotEscapecurve(MWPotential2014,Rrange=[0.01,10.],grid=1001)
At the solar radius, the escape velocity is
>>> from galpy.potential import vesc
>>> vesc(MWPotential2014,1.)
2.3316389848832784
Or, for a local circular velocity of 220 km/s
>>> vesc(MWPotential2014,1.)*220.
512.96057667432126
similar to direct measurements of this (e.g., 2007MNRAS.379..755S and
2014A%26A...562A..91P).
Potentials in galpy
galpy contains a large variety of potentials in galpy.potential
that can be used for orbit integration, the calculation of
action-angle coordinates, as part of steady-state distribution
functions, and to study the properties of gravitational
potentials. This section introduces some of these features.
Potentials and forces
Various 3D and 2D potentials are contained in galpy, list in the
API page. Another way to list the latest overview
of potentials included with galpy is to run
>>> import galpy.potential
>>> print [p for p in dir(galpy.potential) if 'Potential' in p]
['CosmphiDiskPotential',
'DehnenBarPotential',
'DoubleExponentialDiskPotential',
'EllipticalDiskPotential',
'FlattenedPowerPotential',
'HernquistPotential',
....]
(list cut here for brevity). Section Rotation curves explains how to initialize potentials and how to display
the rotation curve of single Potential instances or of combinations of
such instances. Similarly, we can evaluate a Potential instance
>>> from galpy.potential import MiyamotoNagaiPotential
>>> mp= MiyamotoNagaiPotential(a=0.5,b=0.0375,normalize=1.)
>>> mp(1.,0.)
-1.2889062500000001
Most member functions of Potential instances have corresponding
functions in the galpy.potential module that allow them to be
evaluated for lists of multiple Potential
instances. galpy.potential.MWPotential2014 is such a list of three
Potential instances
>>> from galpy.potential import MWPotential2014
>>> print MWPotential2014
[<galpy.potential_src.PowerSphericalPotentialwCutoff.PowerSphericalPotentialwCutoff instance at 0x1089b23b0>, <galpy.potential_src.MiyamotoNagaiPotential.MiyamotoNagaiPotential instance at 0x1089b2320>, <galpy.potential_src.TwoPowerSphericalPotential.NFWPotential instance at 0x1089b2248>]
and we can evaluate the potential by using the evaluatePotentials
function
>>> from galpy.potential import evaluatePotentials
>>> evaluatePotentials(1.,0.,MWPotential2014)
-1.3733506513947895
We can plot the potential of axisymmetric potentials (or of
non-axisymmetric potentials at phi=0) using the plot member
function
which produces the following plot
Similarly, we can plot combinations of Potentials using
plotPotentials, e.g.,
>>> from galpy.potential import plotPotentials
>>> plotPotentials(MWPotential2014,rmin=0.01)
These functions have arguments that can provide custom R and z
ranges for the plot, the number of grid points, the number of
contours, and many other parameters determining the appearance of
these figures.
galpy also allows the forces corresponding to a gravitational
potential to be calculated. Again for the Miyamoto-Nagai Potential
instance from above
>>> mp.Rforce(1.,0.)
-1.0
This value of -1.0 is due to the normalization of the potential such
that the circular velocity is 1. at R=1. Similarly, the vertical force
is zero in the mid-plane
>>> mp.zforce(1.,0.)
-0.0
but not further from the mid-plane
>>> mp.zforce(1.,0.125)
-0.53488743705310848
As explained in Units in galpy, these forces are in
standard galpy units, and we can convert them to physical units using
methods in the galpy.util.bovy_conversion module. For example,
assuming a physical circular velocity of 220 km/s at R=8 kpc
>>> from galpy.util import bovy_conversion
>>> mp.zforce(1.,0.125)*bovy_conversion.force_in_kmsMyr(220.,8.)
-3.3095671288657584 #km/s/Myr
>>> mp.zforce(1.,0.125)*bovy_conversion.force_in_2piGmsolpc2(220.,8.)
-119.72021771473301 #2 \pi G Msol / pc^2
Again, there are functions in galpy.potential that allow for the
evaluation of the forces for lists of Potential instances, such that
>>> from galpy.potential import evaluateRforces
>>> evaluateRforces(1.,0.,MWPotential2014)
-1.0
>>> from galpy.potential import evaluatezforces
>>> evaluatezforces(1.,0.125,MWPotential2014)*bovy_conversion.force_in_2piGmsolpc2(220.,8.)
>>> -69.680720137571114 #2 \pi G Msol / pc^2
We can evaluate the flattening of the potential as
\(\sqrt{|z\,F_R/R\,F_Z|}\) for a Potential instance as well as for
a list of such instances
>>> mp.flattening(1.,0.125)
0.4549542914935209
>>> from galpy.potential import flattening
>>> flattening(MWPotential2014,1.,0.125)
0.61231675305658628
Densities
galpy can also calculate the densities corresponding to gravitational
potentials. For many potentials, the densities are explicitly
implemented, but if they are not, the density is calculated using the
Poisson equation (second derivatives of the potential have to be
implemented for this). For example, for the Miyamoto-Nagai potential,
the density is explicitly implemented
>>> mp.dens(1.,0.)
1.1145444383277576
and we can also calculate this using the Poisson equation
>>> mp.dens(1.,0.,forcepoisson=True)
1.1145444383277574
which are the same to machine precision
>>> mp.dens(1.,0.,forcepoisson=True)-mp.dens(1.,0.)
-2.2204460492503131e-16
Similarly, all of the potentials in galpy.potential.MWPotential2014
have explicitly-implemented densities, so we can do
>>> from galpy.potential import evaluateDensities
>>> evaluateDensities(1.,0.,MWPotential2014)
0.57508603122264867
In physical coordinates, this becomes
>>> evaluateDensities(1.,0.,MWPotential2014)*bovy_conversion.dens_in_msolpc3(220.,8.)
0.1010945632524705 #Msol / pc^3
We can also plot densities
>>> from galpy.potential import plotDensities
>>> plotDensities(MWPotential2014,rmin=0.1,zmax=0.25,zmin=-0.25,nrs=101,nzs=101)
which gives
Another example of this is for an exponential disk potential
>>> from galpy.potential import DoubleExponentialDiskPotential
>>> dp= DoubleExponentialDiskPotential(hr=1./4.,hz=1./20.,normalize=1.)
The density computed using the Poisson equation now requires multiple
numerical integrations, so the agreement between the analytical
density and that computed using the Poisson equation is slightly less good, but still better than a percent
>>> (dp.dens(1.,0.,forcepoisson=True)-dp.dens(1.,0.))/dp.dens(1.,0.)
0.0032522956769123019
The density is
>>> dp.plotDensity(rmin=0.1,zmax=0.25,zmin=-0.25,nrs=101,nzs=101)
and the potential is
>>> dp.plot(rmin=0.1,zmin=-0.25,zmax=0.25)
Clearly, the potential is much less flattened than the density.
Close-to-circular orbits and orbital frequencies
We can also compute the properties of close-to-circular orbits. First
of all, we can calculate the circular velocity and its derivative
>>> mp.vcirc(1.)
1.0
>>> mp.dvcircdR(1.)
-0.163777427566978
or, for lists of Potential instances
>>> from galpy.potential import vcirc
>>> vcirc(MWPotential2014,1.)
1.0
>>> from galpy.potential import dvcircdR
>>> dvcircdR(MWPotential2014,1.)
-0.10091361254334696
We can also calculate the various frequencies for close-to-circular
orbits. For example, the rotational frequency
>>> mp.omegac(0.8)
1.2784598203204887
>>> from galpy.potential import omegac
>>> omegac(MWPotential2014,0.8)
1.2733514576122869
and the epicycle frequency
>>> mp.epifreq(0.8)
1.7774973530267848
>>> from galpy.potential import epifreq
>>> epifreq(MWPotential2014,0.8)
1.7452189766287691
as well as the vertical frequency
>>> mp.verticalfreq(1.0)
3.7859388972001828
>>> from galpy.potential import verticalfreq
>>> verticalfreq(MWPotential2014,1.)
2.7255405754769875
For close-to-circular orbits, we can also compute the radii of the
Lindblad resonances. For example, for a frequency similar to that of
the Milky Way’s bar
>>> mp.lindbladR(5./3.,m='corotation') #args are pattern speed and m of pattern
0.6027911166042229 #~ 5kpc
>>> print mp.lindbladR(5./3.,m=2)
None
>>> mp.lindbladR(5./3.,m=-2)
0.9906190683480501
The None here means that there is no inner Lindblad resonance, the
m=-2 resonance is in the Solar neighborhood (see the section on
the Hercules stream in this documentation).
Using interpolations of potentials
galpy contains a general Potential class interpRZPotential
that can be used to generate interpolations of potentials that can be
used in their stead to speed up calculations when the calculation of
the original potential is computationally expensive (for example, for
the DoubleExponentialDiskPotential). Full details on how to set
this up are given here. Interpolated potentials can
be used anywhere that general three-dimensional galpy potentials can
be used. Some care must be taken with outside-the-interpolation-grid
evaluations for functions that use C to speed up computations.
Adding potentials to the galpy framework
Potentials in galpy can be used in many places such as orbit
integration, distribution functions, or the calculation of
action-angle variables, and in most cases any instance of a potential
class that inherits from the general Potential class (or a list of
such instances) can be given. For example, all orbit integration
routines work with any list of instances of the general Potential
class. Adding new potentials to galpy therefore allows them to be used
everywhere in galpy where general Potential instances can be
used. Adding a new class of potentials to galpy consists of the
following series of steps (some of these are also given in the file
README.dev in the galpy distribution):
- Implement the new potential in a class that inherits from galpy.potential.Potential. The new class should have an __init__ method that sets up the necessary parameters for the class. An amplitude parameter amp= should be taken as an argument for this class and before performing any other setup, the galpy.potential.Potential.__init__(self,amp=amp) method should be called to setup the amplitude. To add support for normalizing the potential to standard galpy units, one can call the galpy.potential.Potential.normalize function at the end of the __init__ function.
The new potential class should implement some of the following
functions:
- _evaluate(self,R,z,phi=0,t=0) which evaluates the
potential itself (without the amp factor, which is added in the
__call__ method of the general Potential class).
- _Rforce(self,R,z,phi=0.,t=0.) which evaluates the radial force
in cylindrical coordinates (-d potential / d R).
- _zforce(self,R,z,phi=0.,t=0.) which evaluates the vertical force
in cylindrical coordinates (-d potential / d z).
- _R2deriv(self,R,z,phi=0.,t=0.) which evaluates the second
(cylindrical) radial derivative of the potential (d^2 potential /
d R^2).
- _z2deriv(self,R,z,phi=0.,t=0.) which evaluates the second
(cylindrical) vertical derivative of the potential (d^2 potential /
d z^2).
- _Rzderiv(self,R,z,phi=0.,t=0.) which evaluates the mixed
(cylindrical) radial and vertical derivative of the potential (d^2
potential / d R d z).
- _dens(self,R,z,phi=0.,t=0.) which evaluates the density. If
not given, the density is computed using the Poisson equation from
the first and second derivatives of the potential (if all are
implemented).
- _mass(self,R,z=0.,t=0.) which evaluates the mass. For
spherical potentials this should give the mass enclosed within the
spherical radius; for axisymmetric potentials this should return
the mass up to R and between -Z and Z. If not given,
the mass is computed by integrating the density (if it is
implemented or can be calculated from the Poisson equation).
- _phiforce(self,R,z,phi=0.,t=0.): the azimuthal force in
cylindrical coordinates (assumed zero if not implemented).
- _phi2deriv(self,R,z,phi=0.,t=0.): the second azimuthal
derivative of the potential in cylindrical coordinates (d^2
potential / d phi^2; assumed zero if not given).
- _Rphideriv(self,R,z,phi=0.,t=0.): the mixed radial and
azimuthal derivative of the potential in cylindrical coordinates
(d^2 potential / d R d phi; assumed zero if not given).
If you want to be able to calculate the concentration for a
potential, you also have to set self._scale to a scale parameter for
your potential.
The code for galpy.potential.MiyamotoNagaiPotential gives a good
template to follow for 3D axisymmetric potentials. Similarly, the
code for galpy.potential.CosmphiDiskPotential provides a good
template for 2D, non-axisymmetric potentials.
After this step, the new potential will work in any part of galpy
that uses pure python potentials. To get the potential to work with
the C implementations of orbit integration or action-angle
calculations, the potential also has to be implemented in C and the
potential has to be passed from python to C.
The __init__ method should be written in such a way that a
relevant object can be initialized using Classname() (i.e.,
there have to be reasonable defaults given for all parameters,
including the amplitude); doing this allows the nose tests for
potentials to automatically check that your Potential’s potential
function, force functions, second derivatives, and density (through
the Poisson equation) are correctly implemented (if they are
implemented). The continuous-integration platform that builds the
galpy codebase upon code pushes will then automatically test all of
this, streamlining push requests of new potentials.
- To add a C implementation of the potential, implement it in a .c file under potential_src/potential_c_ext. Look at potential_src/potential_c_ext/LogarithmicHaloPotential.c for the right format for 3D, axisymmetric potentials, or at potential_src/potential_c_ext/LopsidedDiskPotential.c for 2D, non-axisymmetric potentials.
For orbit integration, the functions such as:
- double LogarithmicHaloPotentialRforce(double R,double Z, double phi,double t,struct potentialArg * potentialArgs)
- double LogarithmicHaloPotentialzforce(double R,double Z, double phi,double t,struct potentialArg * potentialArgs)
are most important. For some of the action-angle calculations
- double LogarithmicHaloPotentialEval(double R,double Z, double phi,double t,struct potentialArg * potentialArgs)
is most important (i.e., for those algorithms that evaluate the potential). The arguments of the potential are passed in a potentialArgs structure that contains args, which are the arguments that should be unpacked. Again, looking at some example code will make this clear. The potentialArgs structure is defined in potential_src/potential_c_ext/galpy_potentials.h.
3. Add the potential’s function declarations to
potential_src/potential_c_ext/galpy_potentials.h
4. (4. and 5. for planar orbit integration) Edit the code under
orbit_src/orbit_c_ext/integratePlanarOrbit.c to set up your new
potential (in the parse_leapFuncArgs function).
5. Edit the code in orbit_src/integratePlanarOrbit.py to set up your
new potential (in the _parse_pot function).
6. Edit the code under orbit_src/orbit_c_ext/integrateFullOrbit.c to
set up your new potential (in the parse_leapFuncArgs_Full function).
7. Edit the code in orbit_src/integrateFullOrbit.py to set up your
new potential (in the _parse_pot function).
8. (for using the actionAngleStaeckel methods in C) Edit the code in
actionAngle_src/actionAngle_c_ext/actionAngle.c to parse the new
potential (in the parse_actionAngleArgs function).
9. Finally, add self.hasC= True to the initialization of the
potential in question (after the initialization of the super class, or
otherwise it will be undone). If you have implemented the necessary
second derivatives for integrating phase-space volumes, also add
self.hasC_dxdv=True.
After following the relevant steps, the new potential class can be
used in any galpy context in which C is used to speed up computations.
Two-dimensional disk distribution functions
galpy contains various disk distribution functions, both in two and
three dimensions. This section introduces the two-dimensional
distribution functions, useful for studying the dynamics of stars that
stay relatively close to the mid-plane of a galaxy. The vertical
motions of these stars may be approximated as being entirely decoupled
from the motion in the plane.
Types of disk distribution functions
galpy contains the following distribution functions for razor-thin
disks: galpy.df.dehnendf and galpy.df.shudf. These are the
distribution functions of Dehnen (1999AJ....118.1201D) and Shu
(1969ApJ...158..505S). Everything
shown below for dehnendf can also be done for shudf.
These disk distribution functions are functions of the energy and the
angular momentum alone. They can be evaluated for orbits, or for a
given energy and angular momentum. At this point, only power-law
rotation curves are supported. A dehnendf instance is initialized
as follows
>>> from galpy.df import dehnendf
>>> dfc= dehnendf(beta=0.)
This initializes a dehnendf instance based on an exponential
surface-mass profile with scale-length 1/3 and an exponential
radial-velocity-dispersion profile with scale-length 1 and a value of
0.2 at R=1. Different parameters for these profiles can be provided as
an initialization keyword. For example,
>>> dfc= dehnendf(beta=0.,profileParams=(1./4.,1.,0.2))
initializes the distribution function with a radial scale length of
1/4 instead.
We can show that these distribution functions have an asymmetric drift
built-in by evaluating the DF at R=1. We first create a set of
orbit-instances and then evaluate the DF at them
>>> from galpy.orbit import Orbit
>>> os= [Orbit([1.,0.,1.+-0.9+1.8/1000*ii]) for ii in range(1001)]
>>> dfro= [dfc(o) for o in os]
>>> plot([1.+-0.9+1.8/1000*ii for ii in range(1001)],dfro)
Or we can plot the two-dimensional density at R=1.
>>> dfro= [[dfc(Orbit([1.,-0.7+1.4/200*jj,1.-0.6+1.2/200*ii])) for jj in range(201)]for ii in range(201)]
>>> dfro= numpy.array(dfro)
>>> from galpy.util.bovy_plot import bovy_dens2d
>>> bovy_dens2d(dfro,origin='lower',cmap='gist_yarg',contours=True,xrange=[-0.7,0.7],yrange=[0.4,1.6],xlabel=r'$v_R$',ylabel=r'$v_T$')
Evaluating moments of the DF
galpy can evaluate various moments of the disk distribution
functions. For example, we can calculate the mean velocities (for the
DF with a scale length of 1/3 above)
>>> dfc.meanvT(1.)
0.91715276979447324
>>> dfc.meanvR(1.)
0.0
and the velocity dispersions
>>> numpy.sqrt(dfc.sigmaR2(1.))
0.19321086259083936
>>> numpy.sqrt(dfc.sigmaT2(1.))
0.15084122011271159
and their ratio
>>> dfc.sigmaR2(1.)/dfc.sigmaT2(1.)
1.6406766813028968
In the limit of zero velocity dispersion (the epicycle approximation)
this ratio should be equal to 2, which we can check as follows
>>> dfccold= dehnendf(beta=0.,profileParams=(1./3.,1.,0.02))
>>> dfccold.sigmaR2(1.)/dfccold.sigmaT2(1.)
1.9947493895454664
We can also calculate higher order moments
>>> dfc.skewvT(1.)
-0.48617143862047763
>>> dfc.kurtosisvT(1.)
0.13338978593181494
>>> dfc.kurtosisvR(1.)
-0.12159407676394096
We already saw above that the velocity dispersion at R=1 is not
exactly equal to the input velocity dispersion (0.19321086259083936
vs. 0.2). Similarly, the whole surface-density and velocity-dispersion
profiles are not quite equal to the exponential input profiles. We can
calculate the resulting surface-mass density profile using
surfacemass, sigmaR2, and sigma2surfacemass. The latter
calculates the product of the velocity dispersion squared and the
surface-mass density. E.g.,
>>> dfc.surfacemass(1.)
0.050820867101511534
We can plot the surface-mass density as follows
>>> Rs= numpy.linspace(0.01,5.,151)
>>> out= [dfc.surfacemass(r) for r in Rs]
>>> plot(Rs, out)
or
>>> plot(Rs,numpy.log(out))
which shows the exponential behavior expected for an exponential
disk. We can compare this to the input surface-mass density
>>> input_out= [dfc.targetSurfacemass(r) for r in Rs]
>>> plot(Rs,numpy.log(input_out)-numpy.log(out))
which shows that there are significant differences between the desired
surface-mass density and the actual surface-mass density. We can do
the same for the velocity-dispersion profile
>>> out= [dfc.sigmaR2(r) for r in Rs]
>>> input_out= [dfc.targetSigma2(r) for r in Rs]
>>> plot(Rs,numpy.log(input_out)-numpy.log(out))
That the input surface-density and velocity-dispersion profiles are
not the same as the output profiles, means that estimates of DF
properties based on these profiles will not be quite
correct. Obviously this is the case for the surface-density and
velocity-dispersion profiles themselves, which have to be explicitly
calculated by integration over the DF rather than by evaluating the
input profiles. This also means that estimates of the asymmetric drift
based on the input profiles will be wrong. We can calculate the
asymmetric drift at R=1 using the asymmetric drift equation derived
from the Jeans equation (eq. 4.228 in Binney & Tremaine 2008), using
the input surface-density and velocity dispersion profiles
>>> dfc.asymmetricdrift(1.)
0.090000000000000024
which should be equal to the circular velocity minus the mean rotational
velocity
>>> 1.-dfc.meanvT(1.)
0.082847230205526756
These are not the same in part because of the difference between the
input and output surface-density and velocity-dispersion profiles (and
because the asymmetricdrift method assumes that the ratio of the
velocity dispersions squared is two using the epicycle approximation;
see above).
Using corrected disk distribution functions
As shown above, for a given surface-mass density and velocity
dispersion profile, the two-dimensional disk distribution functions
only do a poor job of reproducing the desired profiles. We can correct
this by calculating a set of corrections to the input profiles such
that the output profiles more closely resemble the desired profiles
(see 1999AJ....118.1201D). galpy supports
the calculation of these corrections, and comes with some
pre-calculated corrections (these can be found here). For example,
the following initializes a dehnendf with corrections up to 20th
order (the default)
>>> dfc= dehnendf(beta=0.,correct=True)
The following figure shows the difference between the actual
surface-mass density profile and the desired profile for 1, 2, 3, 4,
5, 10, 15, and 20 iterations
and the same for the velocity-dispersion profile
galpy will automatically save any new corrections that you calculate.
All of the methods for an uncorrected disk DF can be used for the
corrected DFs as well. For example, the velocity dispersion is now
>>> numpy.sqrt(dfc.sigmaR2(1.))
0.19999985069451526
and the mean rotation velocity is
>>> dfc.meanvT(1.)
0.90355161181498711
and (correct) asymmetric drift
>>> 1.-dfc.meanvT(1.)
0.09644838818501289
That this still does not agree with the simple dfc.asymmetricdrift
estimate is because of the latter’s using the epicycle approximation
for the ratio of the velocity dispersions.
Oort constants and functions
galpy also contains methods to calculate the Oort functions for
two-dimensional disk distribution functions. These are known as the
Oort constants when measured in the solar neighborhood. They are
combinations of the mean velocities and derivatives thereof. galpy
calculates these by direct integration over the DF and derivatives of
the DF. Thus, we can calculate
>>> dfc= dehnendf(beta=0.)
>>> dfc.oortA(1.)
0.43190780889218749
>>> dfc.oortB(1.)
-0.48524496090228575
The K and C Oort constants are zero for axisymmetric DFs
>>> dfc.oortC(1.)
0.0
>>> dfc.oortK(1.)
0.0
In the epicycle approximation, for a flat rotation curve A =- B =
0.5. The explicit calculates of A and B for warm DFs quantify how
good (or bad) this approximation is
>>> dfc.oortA(1.)+dfc.oortB(1.)
-0.053337152010098254
For the cold DF from above the approximation is much better
>>> dfccold= dehnendf(beta=0.,profileParams=(1./3.,1.,0.02))
>>> dfccold.oortA(1.), dfccold.oortB(1.)
(0.49917556666144003, -0.49992824742490816)
Sampling data from the DF
We can sample from the disk distribution functions using
sample. sample can return either an energy–angular-momentum
pair, or a full orbit initialization. We can sample 4000 orbits for
example as (could take two minutes)
>>> o= dfc.sample(n=4000,returnOrbit=True,nphi=1)
We can then plot the histogram of the sampled radii and compare it to the input surface-mass density profile
>>> Rs= [e.R() for e in o]
>>> hists, bins, edges= hist(Rs,range=[0,2],normed=True,bins=30)
>>> xs= numpy.array([(bins[ii+1]+bins[ii])/2. for ii in range(len(bins)-1)])
>>> plot(xs, xs*exp(-xs*3.)*9.,'r-')
E.g.,
We can also plot the spatial distribution of the sampled disk
>>> xs= [e.x() for e in o]
>>> ys= [e.y() for e in o]
>>> figure()
>>> plot(xs,ys,',')
E.g.,
We can also sample points in a specific radial range (might take a few
minutes)
>>> o= dfc.sample(n=1000,returnOrbit=True,nphi=1,rrange=[0.8,1.2])
and we can plot the distribution of tangential velocities
>>> vTs= [e.vxvv[2] for e in o]
>>> hists, bins, edges= hist(vTs,range=[.5,1.5],normed=True,bins=30)
>>> xs= numpy.array([(bins[ii+1]+bins[ii])/2. for ii in range(len(bins)-1)])
>>> dfro= [dfc(Orbit([1.,0.,x]))/9./numpy.exp(-3.) for x in xs]
>>> plot(xs,dfro,'r-')
The agreement between the sampled distribution and the theoretical
curve is not as good because the sampled distribution has a finite
radial range. If we sample 10,000 points in rrange=[0.95,1.05] the
agreement is better (this takes a long time):
We can also directly sample velocities at a given radius rather than
in a range of radii. Doing this for a correct DF gives
>>> dfc= dehnendf(beta=0.,correct=True)
>>> vrvt= dfc.sampleVRVT(1.,n=10000)
>>> hists, bins, edges= hist(vrvt[:,1],range=[.5,1.5],normed=True,bins=101)
>>> xs= numpy.array([(bins[ii+1]+bins[ii])/2. for ii in range(len(bins)-1)])
>>> dfro= [dfc(Orbit([1.,0.,x])) for x in xs]
>>> plot(xs,dfro/numpy.sum(dfro)/(xs[1]-xs[0]),'r-')
galpy further has support for sampling along a given line of sight in
the disk, which is useful for interpreting surveys consisting of a
finite number of pointings. For example, we can sampled distances
along a given line of sight
>>> ds= dfc.sampledSurfacemassLOS(30./180.*numpy.pi,n=10000)
which is very fast. We can histogram these
>>> hists, bins, edges= hist(ds,range=[0.,3.5],normed=True,bins=101)
and compare it to the predicted distribution, which we can calculate as
>>> xs= numpy.array([(bins[ii+1]+bins[ii])/2. for ii in range(len(bins)-1)])
>>> fd= numpy.array([dfc.surfacemassLOS(d,30.) for d in xs])
>>> plot(xs,fd/numpy.sum(fd)/(xs[1]-xs[0]),'r-')
which shows very good agreement with the sampled distances
galpy can further sample full 4D phase–space coordinates along a
given line of sight through dfc.sampleLOS.
Non-axisymmetric, time-dependent disk distribution functions
galpy also supports the evaluation of non-axisymmetric,
time-dependent two-dimensional DFs. These specific DFs are constructed
by assuming an initial axisymmetric steady state, described by a DF of
the family discussed above, that is then acted upon by a
non-axisymmetric, time-dependent perturbation. The DF at a given time
and phase-space position is evaluated by integrating the orbit
backwards in time in the non-axisymmetric potential until the time of
the initial DF is reached. From Liouville’s theorem, which states that
phase-space volume is conserved along the orbit, we then know that we
can evaluate the non-axisymmetric DF today as the initial DF at the
initial point on the orbit. This procedure was first used by Dehnen
(2000).
This is implemented in galpy as galpy.df.evolveddiskdf. Such a
DF is setup by specifying the initial DF, the non-axisymmetric
potential, and the time of the initial state. For example, we can look
at the effect of an elliptical perturbation to the potential like that
described by Kuijken & Tremaine. To do this, we
set up an elliptical perturbation to a logarithmic potential that is
grown slowly to minimize non-adiabatic effects
>>> from galpy.potential import LogarithmicHaloPotential, EllipticalDiskPotential
>>> lp= LogarithmicHaloPotential(normalize=1.)
>>> ep= EllipticalDiskPotential(twophio=0.05,phib=0.,p=0.,tform=-150.,tsteady=125.)
This perturbation starts to be grown at tform=-150 over a time
period of tsteady=125 time units. We will consider the effect of
this perturbation on a very cold disk (velocity dispersion
\(\sigma_R = 0.0125\,v_c\)) and a warm disk (\(\sigma_R =
0.15\,v_c\)). We set up these two initial DFs
>>> idfcold= dehnendf(beta=0.,profileParams=(1./3.,1.,0.0125))
>>> idfwarm= dehnendf(beta=0.,profileParams=(1./3.,1.,0.15))
and then set up the evolveddiskdf
>>> from galpy.df import evolveddiskdf
>>> edfcold= evolveddiskdf(idfcold,[lp,ep],to=-150.)
>>> edfwarm= evolveddiskdf(idfwarm,[lp,ep],to=-150.)
where we specify that the initial state is at to=-150.
We can now use these evolveddiskdf instances in much the same way
as diskdf instances. One difference is that there is much more
support for evaluating the DF on a grid (to help speed up the rather
slow computations involved). Thus, we can evaluate the mean radial
velocity at R=0.9, phi=22.5 degree, and t=0 by using a grid
>>> mvrcold, gridcold= edfcold.meanvR(0.9,phi=22.5,deg=True,t=0.,grid=True,returnGrid=True,gridpoints=51,nsigma=6.)
>>> mvrwarm, gridwarm= edfcold.meanvR(0.9,phi=22.5,deg=True,t=0.,grid=True,returnGrid=True,gridpoints=51)
>>> print mvrcold, mvrwarm
-0.0358753028951 -0.0294763627935
The cold response agrees well with the analytical calculation, which
predicts that this is \(-0.05/\sqrt{2}\):
>>> print mvrcold+0.05/sqrt(2.)
-0.000519963835811
The warm response is slightly smaller in amplitude
>>> print mvrwarm/mvrcold
0.821633837619
although the numerical uncertainty in mvrwarm is large, because
the grid is not sufficiently fine.
We can then re-use this grid in calculations of other moments of
the DF, e.g.,
>>> print edfcold.meanvT(0.9,phi=22.5,deg=True,t=0.,grid=gridcold)
0.965058551359
>>> print edfwarm.meanvT(0.9,phi=22.5,deg=True,t=0.,grid=gridwarm)
0.915397094614
which returns the mean rotational velocity, and
>>> print edfcold.vertexdev(0.9,phi=22.5,deg=True,t=0.,grid=gridcold)
3.21160878582
>>> print edfwarm.vertexdev(0.9,phi=22.5,deg=True,t=0.,grid=gridwarm)
4.23510254333
which gives the vertex deviation. The reason we have to calculate the
grid out to 6nsigma for the cold response is that the response is
much bigger than the velocity dispersion of the population. This
velocity dispersion is used to automatically to set the grid edges,
but sometimes has to be adjusted to contain the full DF.
evolveddiskdf can also calculate the Oort functions, by directly
calculating the spatial derivatives of the DF. These can also be calculated on a grid, such that we can do
>>> oortacold, gridcold, gridrcold, gridphicold= edfcold.oortA(0.9,phi=22.5,deg=True,t=0.,returnGrids=True,gridpoints=51,derivGridpoints=51,grid=True,derivphiGrid=True,derivRGrid=True,nsigma=6.)
>>> oortawarm, gridwarm, gridrwarm, gridphiwarm= edfwarm.oortA(0.9,phi=22.5,deg=True,t=0.,returnGrids=True,gridpoints=51,derivGridpoints=51,grid=True,derivphiGrid=True,derivRGrid=True)
>>> print oortacold, oortawarm
0.575494559999 0.526389833249
It is clear that these are quite different. The cold calculation is
again close to the analytical prediction, which says that \(A =
A_{\mathrm{axi}}+0.05/(2\sqrt{2})\) where \(A_{\mathrm{axi}} =
1/(2\times0.9)\) in this case:
>>> print oortacold-(0.5/0.9+0.05/2./sqrt(2.))
0.0022613349141670236
These grids can then be re-used for the other Oort functions, for
example,
>>> print edfcold.oortB(0.9,phi=22.5,deg=True,t=0.,grid=gridcold,derivphiGrid=gridphicold,derivRGrid=gridrcold)
-0.574674310521
>>> print edfwarm.oortB(0.9,phi=22.5,deg=True,t=0.,grid=gridwarm,derivphiGrid=gridphiwarm,derivRGrid=gridrwarm)
-0.555546911144
and similar for oortC and oortK. These warm results should
again be considered for illustration only, as the grid is not
sufficiently fine to have a small numerical error.
The grids that have been calculated can also be plotted to show the
full velocity DF. For example,
gives
which demonstrates that the DF is basically the initial DF that has been displaced (by a significant amount compared to the velocity dispersion). The warm velocityd distribution is given by
which returns
The shift of the smooth DF here is much smaller than the velocity
dispersion.
Example: The Hercules stream in the Solar neighborhood as a result of the Galactic bar
We can combine the orbit integration capabilities of galpy with the
provided distribution functions and see the effect of the Galactic bar
on stellar velocities. By backward integrating orbits starting at the
Solar position in a potential that includes the Galactic bar we can
evaluate what the velocity distribution is that we should see today if
the Galactic bar stirred up a steady-state disk. For this we
initialize a flat rotation curve potential and Dehnen’s bar potential
>>> from galpy.potential import LogarithmicHaloPotential, DehnenBarPotential
>>> lp= LogarithmicHaloPotential(normalize=1.)
>>> dp= DehnenBarPotential()
The Dehnen bar potential is initialized to start bar formation four bar
periods before the present day and to have completely formed the bar two
bar periods ago. We can integrate back to the time before
bar-formation:
>>> ts= numpy.linspace(0,dp.tform(),1000)
where dp.tform() is the time of bar-formation (in the usual
time-coordinates).
We initialize orbits on a grid in velocity space and integrate them
>>> ins=[[Orbit([1.,-0.7+1.4/100*jj,1.-0.6+1.2/100*ii,0.]) for jj in range(101)] for ii in range(101)]
>>> int=[[o.integrate(ts,[lp,dp]) for o in j] for j in ins]
We can then evaluate the weight of these orbits by assuming that the
disk was in a steady-state before bar-formation with a Dehnen
distribution function. We evaluate the Dehnen distribution function at
dp.tform() for each of the orbits
>>> dfc= dehnendf(beta=0.,correct=True)
>>> out= [[dfc(o(dp.tform())) for o in j] for j in ins]
>>> out= numpy.array(out)
This gives
>>> from galpy.util.bovy_plot import bovy_dens2d
>>> bovy_dens2d(out,origin='lower',cmap='gist_yarg',contours=True,xrange=[-0.7,0.7],yrange=[0.4,1.6],xlabel=r'$v_R$',ylabel=r'$v_T$')
Now that galpy contains the evolveddiskdf described above,
this whole calculation is encapsulated in this module and can be done
much more easily as
>>> edf= evolveddiskdf(dfc,[lp,dp],to=dp.tform())
>>> mvr, grid= edf.meanvR(1.,grid=True,gridpoints=101,returnGrid=True)
The gridded DF can be accessed as grid.df, which we can plot as before
>>> bovy_dens2d(grid.df.T,origin='lower',cmap='gist_yarg',contours=True,xrange=[grid.vRgrid[0],grid.vRgrid[-1]],yrange=[grid.vTgrid[0],grid.vTgrid[-1]],xlabel=r'$v_R$',ylabel=r'$v_T$')
For more information see 2000AJ....119..800D and
2010ApJ...725.1676B. Note that the
x-axis in the Figure above is defined as minus the x-axis in these
papers.
A closer look at orbit integration
Orbit initialization
Standard initialization
Orbits can be initialized in various coordinate frames. The simplest
initialization gives the initial conditions directly in the
Galactocentric cylindrical coordinate frame (or in the rectangular
coordinate frame in one dimension). Orbit() automatically figures
out the dimensionality of the space from the initial conditions in
this case. In three dimensions initial conditions are given either as
vxvv=[R,vR,vT,z,vz,phi] or one can choose not to specify the
azimuth of the orbit and initialize with
vxvv=[R,vR,vT,z,vz]. Since potentials in galpy are easily
initialized to have a circular velocity of one at a radius equal to
one, initial coordinates are best given as a fraction of the radius at
which one specifies the circular velocity, and initial velocities are
best expressed as fractions of this circular velocity. For example,
>>> o= Orbit(vxvv=[1.,0.1,1.1,0.,0.1,0.])
initializes a fully three-dimensional orbit, while
>>> o= Orbit(vxvv=[1.,0.1,1.1,0.,0.1])
initializes an orbit in which the azimuth is not tracked, as might be
useful for axisymmetric potentials.
In two dimensions, we can similarly specify fully two-dimensional
orbits o=Orbit(vxvv=[R,vR,vT,phi]) or choose not to track the
azimuth and initialize with o= Orbit(vxvv=[R,vR,vT]).
In one dimension we simply initialize with o= Orbit(vxvv=[x,vx]).
Initialization with physical scales
Orbits are normally used in galpy’s natural coordinates. When Orbits
are initialized using a distance scale ro= and a velocity scale
vo=, then many Orbit methods return quantities in physical
coordinates. Specifically, physical distance and velocity scales are
specified as
>>> op= Orbit(vxvv=[1.,0.1,1.1,0.,0.1,0.],ro=8.,vo=220.)
All output quantities will then be automatically be specified in
physical units: kpc for positions, km/s for velocities, (km/s)^2 for
energies and the Jacobi integral, km/s kpc for the angular momentum
o.L() and actions, 1/Gyr for frequencies, and Gyr for times and
periods. See below for examples of this.
Physical units are only used for outputs: internally natural units are
still used and inputs have to also be specified in natural units (for
example, integration times or the time at which an output is requested
must be specified in natural units). If for any output you do not
want the output in physical units, you can specify this by supplying
the keyword argument use_physical=False.
Initialization from observed coordinates
For orbit integration and characterization of observed stars or
clusters, initial conditions can also be specified directly as
observed quantities when radec=True is set. In this case a full
three-dimensional orbit is initialized as o=
Orbit(vxvv=[RA,Dec,distance,pmRA,pmDec,Vlos],radec=True) where RA
and Dec are expressed in degrees, the distance is expressed in kpc,
proper motions are expressed in mas/yr (pmra = pmra’ * cos[Dec] ), and
the line-of-sight velocity is given in km/s. The observed epoch is
currently assumed to be J2000.00. These observed coordinates are
translated to the Galactocentric cylindrical coordinate frame by
assuming a Solar motion that can be specified as either
solarmotion=hogg (default; 2005ApJ...629..268H),
solarmotion=dehnen (1998MNRAS.298..387D) or
solarmotion=shoenrich (2010MNRAS.403.1829S). A circular
velocity can be specified as vo=220 in km/s and a value for the
distance between the Galactic center and the Sun can be given as
ro=8.0 in kpc (e.g., 2012ApJ...759..131B). While the
inputs are given in physical units, the orbit is initialized assuming
a circular velocity of one at the distance of the Sun (that is, the
orbit’s position and velocity is scaled to galpy’s natural units
after converting to the Galactocentric coordinate frame, using the
specified ro= and vo=). The parameters of the coordinate
transformations are stored internally, such that they are
automatically used for relevant outputs (for example, when the RA of
an orbit is requested). An example of all of this is:
>>> o= Orbit(vxvv=[20.,30.,2.,-10.,20.,50.],radec=True,ro=8.,vo=220.)
However, the internally stored position/velocity vector is
>>> print o.vxvv
[1.1476649101960512, 0.20128601278731811, 1.8303776114906387, -0.13107066602923434, 0.58171049004255293, 0.14071341020496472]
and is therefore in natural units.
Similarly, one can also initialize orbits from Galactic coordinates
using o= Orbit(vxvv=[glon,glat,distance,pmll,pmbb,Vlos],lb=True),
where glon and glat are Galactic longitude and latitude expressed in
degrees, and the proper motions are again given in mas/yr ((pmll =
pmll’ * cos[glat] ):
>>> o= Orbit(vxvv=[20.,30.,2.,-10.,20.,50.],lb=True,ro=8.,vo=220.)
>>> print o.vxvv
[0.79998509943955398, 0.075939950035477488, 0.52838231795389867, 0.12812499999999999, 0.89052135379600328, 0.092696334097541536]
When radec=True or lb=True is set, velocities can also be specified in
Galactic coordinates if UVW=True is set. The input is then
vxvv=[RA,Dec,distance,U,V,W], where the velocities are expressed
in km/s. U is, as usual, defined as -vR (minus vR).
When orbits are initialized using radec=True or lb=True,
physical scales ro= and vo= are automatically specified
(because they have defaults of ro=8 and vo=220). Therefore,
all output quantities will be specified in physical units (see
above). If you do want to get outputs in galpy’s natural coordinates,
you can turn this behavior off by doing
>>> o.turn_physical_off()
All outputs will then be specified in galpy’s natural coordinates.
Orbit integration
After an orbit is initialized, we can integrate it for a set of times
ts, given as a numpy array. For example, in a simple logarithmic
potential we can do the following
>>> from galpy.potential import LogarithmicHaloPotential
>>> lp= LogarithmicHaloPotential(normalize=1.)
>>> o= Orbit(vxvv=[1.,0.1,1.1,0.,0.1,0.])
>>> import numpy
>>> ts= numpy.linspace(0,100,10000)
>>> o.integrate(ts,lp)
to integrate the orbit from t=0 to t=100, saving the orbit at
10000 instances.
If we initialize the Orbit using a distance scale ro= and a
velocity scale vo=, then Orbit plots and outputs will use physical
coordinates (currently, times, positions, and velocities)
>>> op= Orbit(vxvv=[1.,0.1,1.1,0.,0.1,0.],ro=8.,vo=220.) #Use Vc=220 km/s at R= 8 kpc as the normalization
>>> op.integrate(ts,lp) #times are still specified in natural coordinates
Displaying the orbit
After integrating the orbit, it can be displayed by using the
plot() function. The quantities that are plotted when plot()
is called depend on the dimensionality of the orbit: in 3D the (R,z)
projection of the orbit is shown; in 2D either (X,Y) is plotted if the
azimuth is tracked and (R,vR) is shown otherwise; in 1D (x,vx) is
shown. E.g., for the example given above,
gives
If we do the same for the Orbit that has physical distance and
velocity scales associated with it, we get the following
If we call op.plot(use_physical=False), the quantities will be
displayed in natural galpy coordinates.
Other projections of the orbit can be displayed by specifying the
quantities to plot. E.g.,
>>> o.plot(d1='x',d2='y')
gives the projection onto the plane of the orbit:
while
>>> o.plot(d1='R',d2='vR')
gives the projection onto (R,vR):
We can also plot the orbit in other coordinate systems such as
Galactic longitude and latitude
>>> o.plot('k.',d1='ll',d2='bb')
which shows
or RA and Dec
>>> o.plot('k.',d1='ra',d2='dec')
See the documentation of the o.plot function and the o.ra(), o.ll(),
etc. functions on how to provide the necessary parameters for the
coordinate transformations.
Orbit characterization
The properties of the orbit can also be found using galpy. For
example, we can calculate the peri- and apocenter radii of an orbit,
its eccentricity, and the maximal height above the plane of the orbit
>>> o.rap(), o.rperi(), o.e(), o.zmax()
(1.2581455175173673,0.97981663263371377,0.12436710999105324,0.11388132751079502)
We can also calculate the energy of the orbit, either in the potential
that the orbit was integrated in, or in another potential:
>>> o.E(), o.E(pot=mp)
(0.6150000000000001, -0.67390625000000015)
where mp is the Miyamoto-Nagai potential of Introduction:
Rotation curves.
For the Orbit op that was initialized above with a distance scale
ro= and a velocity scale vo=, these outputs are all in
physical units
>>> op.rap(), op.rperi(), op.e(), op.zmax()
(10.065158988860341,7.8385312810643057,0.12436696983841462,0.91105035688072711) #kpc
>>> op.E(), op.E(pot=mp)
(29766.000000000004, -32617.062500000007) #(km/s)^2
We can also show the energy as a function of time (to check energy
conservation)
gives
We can specify another quantity to plot the energy against by
specifying d1=. We can also show the vertical energy, for example,
as a function of R
>>> o.plotEz(d1='R',normed=True)
Often, a better approximation to an integral of the motion is given by
Ez/sqrt(density[R]). We refer to this quantity as EzJz and we can plot its
behavior
>>> o.plotEzJz(d1='R',normed=True)
Accessing the raw orbit
The value of R, vR, vT, z, vz, x, vx,
y, vy, phi, and vphi at any time can be obtained by
calling the corresponding function with as argument the time (the same
holds for other coordinates ra, dec, pmra, pmdec,
vra, vdec, ll, bb, pmll, pmbb, vll,
vbb, vlos, dist, helioX, helioY, helioZ,
U, V, and W). If no time is given the initial condition is
returned, and if a time is requested at which the orbit was not saved
spline interpolation is used to return the value. Examples include
>>> o.R(1.)
1.1545076874679474
>>> o.phi(99.)
88.105603035901169
>>> o.ra(2.,obs=[8.,0.,0.],ro=8.)
array([ 285.76403985])
>>> o.helioX(5.)
array([ 1.24888927])
>>> o.pmll(10.,obs=[8.,0.,0.,0.,245.,0.],ro=8.,vo=230.)
array([-6.45263888])
For the Orbit op that was initialized above with a distance scale
ro= and a velocity scale vo=, the first of these would be
>>> op.R(1.)
9.2360614837829225 #kpc
which we can also access in natural coordinates as
>>> op.R(1.,use_physical=False)
1.1545076854728653
We can also specify a different distance or velocity scale on the fly,
e.g.,
>>> op.R(1.,ro=4.) #different velocity scale would be vo=
4.6180307418914612
We can also initialize an Orbit instance using the phase-space
position of another Orbit instance evaulated at time t. For
example,
will initialize a new Orbit instance with as initial condition the phase-space position of orbit o at time=10..
The whole orbit can also be obtained using the function getOrbit
which returns a matrix of phase-space points with dimensions [ntimes,ndim].
Fast orbit integration
The standard orbit integration is done purely in python using standard
scipy integrators. When fast orbit integration is needed for batch
integration of a large number of orbits, a set of orbit integration
routines are written in C that can be accessed for most potentials, as
long as they have C implementations, which can be checked by using the
attribute hasC
>>> mp= MiyamotoNagaiPotential(a=0.5,b=0.0375,amp=1.,normalize=1.)
>>> mp.hasC
True
Fast C integrators can be accessed through the method= keyword of
the orbit.integrate method. Currently available integrators are
which are Runge-Kutta and Dormand-Prince methods. There are also a
number of symplectic integrators available
- leapfrog_c
- symplec4_c
- symplec6_c
The higher order symplectic integrators are described in Yoshida
(1993).
For most applications I recommend dopr54_c. For example, compare
>>> o= Orbit(vxvv=[1.,0.1,1.1,0.,0.1])
>>> timeit(o.integrate(ts,mp))
1 loops, best of 3: 553 ms per loop
>>> timeit(o.integrate(ts,mp,method='dopr54_c'))
galpyWarning: Using C implementation to integrate orbits
10 loops, best of 3: 25.6 ms per loop
As this example shows, galpy will issue a warning that C is being
used. Speed-ups by a factor of 20 are typical.
Integration of the phase-space volume
galpy further supports the integration of the phase-space volume
through the method integrate_dxdv, although this is currently only
implemented for two-dimensional orbits (planarOrbit). As an
example, we can check Liouville’s theorem explicitly. We initialize
the orbit
>>> o= Orbit(vxvv=[1.,0.1,1.1,0.])
and then integrate small deviations in each of the four
phase-space directions
>>> ts= numpy.linspace(0.,28.,1001) #~1 Gyr at the Solar circle
>>> o.integrate_dxdv([1.,0.,0.,0.],ts,mp,method='dopr54_c',rectIn=True,rectOut=True)
>>> dx= o.getOrbit_dxdv()[-1,:] # evolution of dxdv[0] along the orbit
>>> o.integrate_dxdv([0.,1.,0.,0.],ts,mp,method='dopr54_c',rectIn=True,rectOut=True)
>>> dy= o.getOrbit_dxdv()[-1,:]
>>> o.integrate_dxdv([0.,0.,1.,0.],ts,mp,method='dopr54_c',rectIn=True,rectOut=True)
>>> dvx= o.getOrbit_dxdv()[-1,:]
>>> o.integrate_dxdv([0.,0.,0.,1.],ts,mp,method='dopr54_c',rectIn=True,rectOut=True)
>>> dvy= o.getOrbit_dxdv()[-1,:]
We can then compute the determinant of the Jacobian of the mapping
defined by the orbit integration from time zero to the final time
>>> tjac= numpy.linalg.det(numpy.array([dx,dy,dvx,dvy]))
This determinant should be equal to one
>>> print tjac
0.999999991189
>>> numpy.fabs(tjac-1.) < 10.**-8.
True
The calls to integrate_dxdv above set the keywords rectIn= and
rectOut= to True, as the default input and output uses phase-space
volumes defined as (dR,dvR,dvT,dphi) in cylindrical coordinates. When
rectIn or rectOut is set, the in- or output is in rectangular
coordinates ([x,y,vx,vy] in two dimensions).
Implementing the phase-space integration for three-dimensional
FullOrbit instances is straightforward and is part of the longer
term development plan for galpy. Let the main developer know if
you would like this functionality, or better yet, implement it
yourself in a fork of the code and send a pull request!
Example: The eccentricity distribution of the Milky Way’s thick disk
A straightforward application of galpy’s orbit initialization and
integration capabilities is to derive the eccentricity distribution of
a set of thick disk stars. We start by downloading the sample of SDSS
SEGUE (2009AJ....137.4377Y) thick disk
stars compiled by Dierickx et al. (2010arXiv1009.1616D) at
http://www.mpia-hd.mpg.de/homes/rix/Data/Dierickx-etal-tab2.txt
After reading in the data (RA,Dec,distance,pmRA,pmDec,vlos; see above)
as a vector vxvv with dimensions [6,ndata] we (a) define the
potential in which we want to integrate the orbits, and (b) integrate
each orbit and save its eccentricity (running this for all 30,000-ish
stars will take about half an hour)
>>> lp= LogarithmicHaloPotential(normalize=1.)
>>> ts= nu.linspace(0.,20.,10000)
>>> mye= nu.zeros(ndata)
>>> for ii in range(len(e)):
... o= Orbit(vxvv[ii,:],radec=True,vo=220.,ro=8.) #Initialize
... o.integrate(ts,lp) #Integrate
... mye[ii]= o.e() #Calculate eccentricity
We then find the following eccentricity distribution
The eccentricity calculated by galpy compare well with those
calculated by Dierickx et al., except for a few objects
The script that calculates and plots everything can be downloaded
here.
Action-angle coordinates
galpy can calculate actions and angles for a large variety of
potentials (any time-independent potential in principle). These are
implemented in a separate module galpy.actionAngle, and the
preferred method for accessing them is through the routines in this
module. There is also some support for accessing the actionAngle
routines as methods of the Orbit class.
Action-angle coordinates can be calculated for the following
potentials/approximations:
- Isochrone potential
- Spherical potentials
- Adiabatic approximation
- Staeckel approximation
- A general orbit-integration-based technique
There are classes corresponding to these different
potentials/approximations and actions, frequencies, and angles can
typically be calculated using these three methods:
- __call__: returns the actions
- actionsFreqs: returns the actions and the frequencies
- actionsFreqsAngles: returns the actions, frequencies, and angles
These are not all implemented for each of the cases above yet.
The adiabatic and Staeckel approximation have also been implemented in
C and using grid-based interpolation, for extremely fast action-angle
calculations (see below).
Action-angle coordinates for the isochrone potential
The isochrone potential is the only potential for which all of the
actions, frequencies, and angles can be calculated analytically. We
can do this in galpy by doing
>>> from galpy.potential import IsochronePotential
>>> from galpy.actionAngle import actionAngleIsochrone
>>> ip= IsochronePotential(b=1.,normalize=1.)
>>> aAI= actionAngleIsochrone(ip=ip)
aAI is now an instance that can be used to calculate action-angle
variables for the specific isochrone potential ip. Calling this
instance returns \((J_R,L_Z,J_Z)\)
>>> aAI(1.,0.1,1.1,0.1,0.) #inputs R,vR,vT,z,vz
(array([ 0.00713759]), array([ 1.1]), array([ 0.00553155]))
or for a more eccentric orbit
>>> aAI(1.,0.5,1.3,0.2,0.1)
(array([ 0.13769498]), array([ 1.3]), array([ 0.02574507]))
Note that we can also specify phi, but this is not necessary
>>> aAI(1.,0.5,1.3,0.2,0.1,0.)
(array([ 0.13769498]), array([ 1.3]), array([ 0.02574507]))
We can likewise calculate the frequencies as well
>>> aAI.actionsFreqs(1.,0.5,1.3,0.2,0.1,0.)
(array([ 0.13769498]),
array([ 1.3]),
array([ 0.02574507]),
array([ 1.29136096]),
array([ 0.79093738]),
array([ 0.79093738]))
The output is \((J_R,L_Z,J_Z,\Omega_R,\Omega_\phi,\Omega_Z)\). For
any spherical potential, \(\Omega_\phi =
\mathrm{sgn}(L_Z)\Omega_Z\), such that the last two frequencies are the
same.
We obtain the angles as well by calling
>>> aAI.actionsFreqsAngles(1.,0.5,1.3,0.2,0.1,0.)
(array([ 0.13769498]),
array([ 1.3]),
array([ 0.02574507]),
array([ 1.29136096]),
array([ 0.79093738]),
array([ 0.79093738]),
array([ 0.57101518]),
array([ 5.96238847]),
array([ 1.24999949]))
The output here is
\((J_R,L_Z,J_Z,\Omega_R,\Omega_\phi,\Omega_Z,\theta_R,\theta_\phi,\theta_Z)\).
To check that these are good action-angle variables, we can calculate
them along an orbit
>>> from galpy.orbit import Orbit
>>> o= Orbit([1.,0.5,1.3,0.2,0.1,0.])
>>> ts= numpy.linspace(0.,100.,1001)
>>> o.integrate(ts,ip)
>>> jfa= aAI.actionsFreqsAngles(o.R(ts),o.vR(ts),o.vT(ts),o.z(ts),o.vz(ts),o.phi(ts))
which works because we can provide arrays for the R etc. inputs.
We can then check that the actions are constant over the orbit
>>> plot(ts,numpy.log10(numpy.fabs((jfa[0]-numpy.mean(jfa[0])))))
>>> plot(ts,numpy.log10(numpy.fabs((jfa[1]-numpy.mean(jfa[1])))))
>>> plot(ts,numpy.log10(numpy.fabs((jfa[2]-numpy.mean(jfa[2])))))
which gives
The actions are all conserved. The angles increase linearly with time
>>> plot(ts,jfa[6],'b.')
>>> plot(ts,jfa[7],'g.')
>>> plot(ts,jfa[8],'r.')
Action-angle coordinates for spherical potentials
Action-angle coordinates for any spherical potential can be calculated
using a few orbit integrations. These are implemented in galpy in the
actionAngleSpherical module. For example, we can do
>>> from galpy.potential import LogarithmicHaloPotential
>>> lp= LogarithmicHaloPotential(normalize=1.)
>>> from galpy.actionAngle import actionAngleSpherical
>>> aAS= actionAngleSpherical(pot=lp)
For the same eccentric orbit as above we find
>>> aAS(1.,0.5,1.3,0.2,0.1,0.)
(array([ 0.22022112]), array([ 1.3]), array([ 0.02574507]))
>>> aAS.actionsFreqs(1.,0.5,1.3,0.2,0.1,0.)
(array([ 0.22022112]),
array([ 1.3]),
array([ 0.02574507]),
array([ 0.87630459]),
array([ 0.60872881]),
array([ 0.60872881]))
>>> aAS.actionsFreqsAngles(1.,0.5,1.3,0.2,0.1,0.)
(array([ 0.22022112]),
array([ 1.3]),
array([ 0.02574507]),
array([ 0.87630459]),
array([ 0.60872881]),
array([ 0.60872881]),
array([ 0.40443857]),
array([ 5.85965048]),
array([ 1.1472615]))
We can again check that the actions are conserved along the orbit and
that the angles increase linearly with time:
>>> o.integrate(ts,lp)
>>> jfa= aAS.actionsFreqsAngles(o.R(ts),o.vR(ts),o.vT(ts),o.z(ts),o.vz(ts),o.phi(ts),fixed_quad=True)
where we use fixed_quad=True for a faster evaluation of the
required one-dimensional integrals using Gaussian quadrature. We then
plot the action fluctuations
>>> plot(ts,numpy.log10(numpy.fabs((jfa[0]-numpy.mean(jfa[0])))))
>>> plot(ts,numpy.log10(numpy.fabs((jfa[1]-numpy.mean(jfa[1])))))
>>> plot(ts,numpy.log10(numpy.fabs((jfa[2]-numpy.mean(jfa[2])))))
which gives
showing that the actions are all conserved. The angles again increase
linearly with time
>>> plot(ts,jfa[6],'b.')
>>> plot(ts,jfa[7],'g.')
>>> plot(ts,jfa[8],'r.')
We can check the spherical action-angle calculations against the
analytical calculations for the isochrone potential. Starting again
from the isochrone potential used in the previous section
>>> ip= IsochronePotential(b=1.,normalize=1.)
>>> aAI= actionAngleIsochrone(ip=ip)
>>> aAS= actionAngleSpherical(pot=ip)
we can compare the actions, frequencies, and angles computed using
both
>>> aAI.actionsFreqsAngles(1.,0.5,1.3,0.2,0.1,0.)
(array([ 0.13769498]),
array([ 1.3]),
array([ 0.02574507]),
array([ 1.29136096]),
array([ 0.79093738]),
array([ 0.79093738]),
array([ 0.57101518]),
array([ 5.96238847]),
array([ 1.24999949]))
>>> aAS.actionsFreqsAngles(1.,0.5,1.3,0.2,0.1,0.)
(array([ 0.13769498]),
array([ 1.3]),
array([ 0.02574507]),
array([ 1.29136096]),
array([ 0.79093738]),
array([ 0.79093738]),
array([ 0.57101518]),
array([ 5.96238838]),
array([ 1.2499994]))
or more explicitly comparing the two
>>> [r-s for r,s in zip(aAI.actionsFreqsAngles(1.,0.5,1.3,0.2,0.1,0.),aAS.actionsFreqsAngles(1.,0.5,1.3,0.2,0.1,0.))]
[array([ 6.66133815e-16]),
array([ 0.]),
array([ 0.]),
array([ -4.53851845e-10]),
array([ 4.74775219e-10]),
array([ 4.74775219e-10]),
array([ -1.65965242e-10]),
array([ 9.04759645e-08]),
array([ 9.04759649e-08])]
Action-angle coordinates using the adiabatic approximation
For non-spherical, axisymmetric potentials galpy contains multiple
methods for calculating approximate action–angle coordinates. The
simplest of those is the adiabatic approximation, which works well for
disk orbits that do not go too far from the plane, as it assumes that
the vertical motion is decoupled from that in the plane (e.g.,
2010MNRAS.401.2318B).
Setup is similar as for other actionAngle objects
>>> from galpy.potential import MWPotential2014
>>> from galpy.actionAngle import actionAngleAdiabatic
>>> aAA= actionAngleAdiabatic(pot=MWPotential2014)
and evaluation then proceeds similarly as before
>>> aAA(1.,0.1,1.1,0.,0.05)
(0.01351896260559274, 1.1, 0.0004690133479435352)
We can again check that the actions are conserved along the orbit
>>> from galpy.orbit import Orbit
>>> ts=numpy.linspace(0.,100.,1001)
>>> o= Orbit([1.,0.1,1.1,0.,0.05])
>>> o.integrate(ts,MWPotential2014)
>>> js= aAA(o.R(ts),o.vR(ts),o.vT(ts),o.z(ts),o.vz(ts))
This takes a while. The adiabatic approximation is also implemented in
C, which leads to great speed-ups. Here is how to use it
>>> timeit(aAA(1.,0.1,1.1,0.,0.05))
10 loops, best of 3: 73.7 ms per loop
>>> aAA= actionAngleAdiabatic(pot=MWPotential2014,c=True)
>>> timeit(aAA(1.,0.1,1.1,0.,0.05))
1000 loops, best of 3: 1.3 ms per loop
or about a 50 times speed-up. For arrays the speed-up is even more
impressive
>>> s= numpy.ones(100)
>>> timeit(aAA(1.*s,0.1*s,1.1*s,0.*s,0.05*s))
10 loops, best of 3: 37.8 ms per loop
>>> aAA= actionAngleAdiabatic(pot=MWPotential2014) #back to no C
>>> timeit(aAA(1.*s,0.1*s,1.1*s,0.*s,0.05*s))
1 loops, best of 3: 7.71 s per loop
or a speed-up of 200! Back to the previous example, you can run it
with c=True to speed up the computation
>>> aAA= actionAngleAdiabatic(pot=MWPotential2014,c=True)
>>> js= aAA(o.R(ts),o.vR(ts),o.vT(ts),o.z(ts),o.vz(ts))
We can plot the radial- and vertical-action fluctuation as a function
of time
>>> plot(ts,numpy.log10(numpy.fabs((js[0]-numpy.mean(js[0]))/numpy.mean(js[0]))))
>>> plot(ts,numpy.log10(numpy.fabs((js[2]-numpy.mean(js[2]))/numpy.mean(js[2]))))
which gives
The radial action is conserved to about half a percent, the vertical
action to two percent.
Another way to speed up the calculation of actions using the adiabatic
approximation is to tabulate the actions on a grid in (approximate)
integrals of the motion and evaluating new actions by interpolating on
this grid. How this is done in practice is described in detail in the
galpy paper. To setup this grid-based interpolation method, which is
contained in actionAngleAdiabaticGrid, do
>>> from galpy.actionAngle import actionAngleAdiabaticGrid
>>> aAG= actionAngleAdiabaticGrid(pot=MWPotential2014,nR=31,nEz=31,nEr=51,nLz=51,c=True)
where c=True specifies that we use the C implementation of
actionAngleAdiabatic for speed. We can now evaluate in the same
was as before, for example
>>> aAA(1.,0.1,1.1,0.,0.05), aAG(1.,0.1,1.1,0.,0.05)
((array([ 0.01352523]), array([ 1.1]), array([ 0.00046909])),
(0.013527010324238781, 1.1, 0.00047747359874375148))
which agree very well. To look at the timings, we first switch back to
not using C and then list all of the relevant timings:
>>> aAA= actionAngleAdiabatic(pot=MWPotential2014,c=False)
# Not using C, direct calculation
>>> timeit(aAA(1.*s,0.1*s,1.1*s,0.*s,0.05*s))
1 loops, best of 3: 9.05 s per loop
>>> aAA= actionAngleAdiabatic(pot=MWPotential2014,c=True)
# Using C, direct calculation
>>> timeit(aAA(1.*s,0.1*s,1.1*s,0.*s,0.05*s))
10 loops, best of 3: 39.7 ms per loop
# Grid-based calculation
>>> timeit(aAG(1.*s,0.1*s,1.1*s,0.*s,0.05*s))
1000 loops, best of 3: 1.09 ms per loop
Thus, in this example (and more generally) the grid-based calculation
is significantly faster than even the direct implementation in C. The
overall speed up between the direct Python version and the grid-based
version is larger than 8,000; the speed up between the direct C
version and the grid-based version is 36. For larger arrays of input
phase-space positions, the latter speed up can increase to 150. For
simpler, fully analytical potentials the speed up will be slightly
less, but for MWPotential2014 and other more complicated
potentials (such as those involving a double-exponential disk), the
overhead of setting up the grid is worth it when evaluating more than
a few thousand actions.
The adiabatic approximation works well for orbits that stay close to
the plane. The orbit we have been considering so far only reaches a
height two percent of \(R_0\), or about 150 pc for \(R_0 = 8\)
kpc.
>>> o.zmax()*8.
0.17903686455491979
For orbits that reach distances of a kpc and more from the plane, the
adiabatic approximation does not work as well. For example,
>>> o= Orbit([1.,0.1,1.1,0.,0.25])
>>> o.integrate(ts,MWPotential2014)
>>> o.zmax()*8.
1.3506059038621048
and we can again calculate the actions along the orbit
>>> js= aAA(o.R(ts),o.vR(ts),o.vT(ts),o.z(ts),o.vz(ts))
>>> plot(ts,numpy.log10(numpy.fabs((js[0]-numpy.mean(js[0]))/numpy.mean(js[0]))))
>>> plot(ts,numpy.log10(numpy.fabs((js[2]-numpy.mean(js[2]))/numpy.mean(js[2]))))
which gives
The radial action is now only conserved to about ten percent and the
vertical action to approximately five percent.
Warning
Frequencies and angles using the adiabatic approximation are not implemented at this time.
Action-angle coordinates using the Staeckel approximation
A better approximation than the adiabatic one is to locally
approximate the potential as a Staeckel potential, for which actions,
frequencies, and angles can be calculated through numerical
integration. galpy contains an implementation of the algorithm of
Binney (2012; 2012MNRAS.426.1324B), which
accomplishes the Staeckel approximation for disk-like (i.e., oblate)
potentials without explicitly fitting a Staeckel potential. For all
intents and purposes the adiabatic approximation is made obsolete by
this new method, which is as fast and more precise. The only advantage
of the adiabatic approximation over the Staeckel approximation is that
the Staeckel approximation requires the user to specify a focal
length \(\Delta\) to be used in the Staeckel
approximation. However, this focal length can be easily estimated from
the second derivatives of the potential (see Sanders 2012;
2012MNRAS.426..128S).
Starting from the second orbit example in the adiabatic section above,
we first estimate a good focal length of the MWPotential2014 to
use in the Staeckel approximation. We do this by averaging (through
the median) estimates at positions around the orbit (which we
integrated in the example above)
>>> from galpy.actionAngle import estimateDeltaStaeckel
>>> estimateDeltaStaeckel(o.R(ts),o.z(ts),pot=MWPotential2014)
0.40272708556203662
We will use \(\Delta = 0.4\) in what follows. We set up the
actionAngleStaeckel object
>>> from galpy.actionAngle import actionAngleStaeckel
>>> aAS= actionAngleStaeckel(pot=MWPotential2014,delta=0.4,c=False) #c=True is the default
and calculate the actions
>>> aAS(o.R(),o.vR(),o.vT(),o.z(),o.vz())
(0.019212848866725911, 1.1000000000000001, 0.015274597971510892)
The adiabatic approximation from above gives
>>> aAA(o.R(),o.vR(),o.vT(),o.z(),o.vz())
(array([ 0.01686478]), array([ 1.1]), array([ 0.01590001]))
The actionAngleStaeckel calculations are sped up in two ways. First,
the action integrals can be calculated using Gaussian quadrature by
specifying fixed_quad=True
>>> aAS(o.R(),o.vR(),o.vT(),o.z(),o.vz(),fixed_quad=True)
(0.01922167296633687, 1.1000000000000001, 0.015276825017286706)
which in itself leads to a ten times speed up
>>> timeit(aAS(o.R(),o.vR(),o.vT(),o.z(),o.vz(),fixed_quad=False))
10 loops, best of 3: 129 ms per loop
>>> timeit(aAS(o.R(),o.vR(),o.vT(),o.z(),o.vz(),fixed_quad=True))
100 loops, best of 3: 10.3 ms per loop
Second, the actionAngleStaeckel calculations have also been
implemented in C, which leads to even greater speed-ups, especially
for arrays
>>> aAS= actionAngleStaeckel(pot=MWPotential2014,delta=0.4,c=True)
>>> s= numpy.ones(100)
>>> timeit(aAS(1.*s,0.1*s,1.1*s,0.*s,0.05*s))
10 loops, best of 3: 35.1 ms per loop
>>> aAS= actionAngleStaeckel(pot=MWPotential2014,delta=0.4,c=False) #back to no C
>>> timeit(aAS(1.*s,0.1*s,1.1*s,0.*s,0.05*s,fixed_quad=True))
1 loops, best of 3: 496 ms per loop
or a fifteen times speed up. The speed up is not that large because
the bulge model in MWPotential2014 requires expensive special
functions to be evaluated. Computations could be sped up ten times
more when using a simpler bulge model.
Similar to actionAngleAdiabaticGrid, we can also tabulate the
actions on a grid of (approximate) integrals of the motion and
interpolate over this look-up table when evaluating new actions. The
details of how this look-up table is setup and used are again fully
explained in the galpy paper. To use this grid-based Staeckel
approximation, contained in actionAngleStaeckelGrid, do
>>> from galpy.actionAngle import actionAngleStaeckelGrid
>>> aASG= actionAngleStaeckelGrid(pot=MWPotential2014,delta=0.4,nE=51,npsi=51,nLz=61,c=True)
where c=True makes sure that we use the C implementation of the
Staeckel method to calculate the grid. Because this is a fully
three-dimensional grid, setting up the grid takes longer than it does
for the adiabatic method (which only uses two two-dimensional
grids). We can then evaluate actions as before
>>> aAS(o.R(),o.vR(),o.vT(),o.z(),o.vz()), aASG(o.R(),o.vR(),o.vT(),o.z(),o.vz())
((0.019212848866725911, 1.1000000000000001, 0.015274597971510892),
(0.019221119033345408, 1.1000000000000001, 0.015022528662310393))
These actions agree very well. We can compare the timings of these
methods as above
>>> timeit(aAS(1.*s,0.1*s,1.1*s,0.*s,0.05*s,fixed_quad=True))
1 loops, best of 3: 576 ms per loop # Not using C, direct calculation
>>> aAS= actionAngleStaeckel(pot=MWPotential2014,delta=0.4,c=True)
>>> timeit(aAS(1.*s,0.1*s,1.1*s,0.*s,0.05*s))
100 loops, best of 3: 17.8 ms per loop # Using C, direct calculation
>>> timeit(aASG(1.*s,0.1*s,1.1*s,0.*s,0.05*s))
100 loops, best of 3: 3.45 ms per loop # Grid-based calculation
This demonstrates that the grid-based interpolation again leeds to a
significant speed up, even over the C implementation of the direct
calculation. This speed up becomes more significant for larger array
input, although it saturates at about 25 times (at least for
MWPotential2014).
We can now go back to checking that the actions are conserved along
the orbit (going back to the c=False version of
actionAngleStaeckel)
>>> aAS= actionAngleStaeckel(pot=MWPotential2014,delta=0.4,c=False)
>>> js= aAS(o.R(ts),o.vR(ts),o.vT(ts),o.z(ts),o.vz(ts),fixed_quad=True)
>>> plot(ts,numpy.log10(numpy.fabs((js[0]-numpy.mean(js[0]))/numpy.mean(js[0]))))
>>> plot(ts,numpy.log10(numpy.fabs((js[2]-numpy.mean(js[2]))/numpy.mean(js[2]))))
which gives
The radial action is now conserved to better than a percent and the
vertical action to only a fraction of a percent. Clearly, this is much
better than the five to ten percent errors found for the adiabatic
approximation above.
For the Staeckel approximation we can also calculate frequencies and
angles through the actionsFreqs and actionsFreqsAngles
methods.
Warning
Frequencies and angles using the Staeckel approximation
are only implemented in C. So use c=True in the setup of the
actionAngleStaeckel object.
Warning
Angles using the Staeckel approximation in galpy are such
that (a) the radial angle starts at zero at pericenter and
increases then going toward apocenter; (b) the vertical angle
starts at zero at z=0 and increases toward positive zmax. The
latter is a different convention from that in Binney (2012), but is
consistent with that in actionAngleIsochrone and
actionAngleSpherical.
>>> aAS= actionAngleStaeckel(pot=MWPotential2014,delta=0.4,c=True)
>>> o= Orbit([1.,0.1,1.1,0.,0.25,0.]) #need to specify phi for angles
>>> aAS.actionsFreqsAngles(o.R(),o.vR(),o.vT(),o.z(),o.vz(),o.phi())
(array([ 0.01922167]),
array([ 1.1]),
array([ 0.01527683]),
array([ 1.11317796]),
array([ 0.82538032]),
array([ 1.34126138]),
array([ 0.37758087]),
array([ 6.17833493]),
array([ 6.13368239]))
and we can check that the angles increase linearly along the orbit
>>> o.integrate(ts,MWPotential2014)
>>> jfa= aAS.actionsFreqsAngles(o.R(ts),o.vR(ts),o.vT(ts),o.z(ts),o.vz(ts),o.phi(ts))
>>> plot(ts,jfa[6],'b.')
>>> plot(ts,jfa[7],'g.')
>>> plot(ts,jfa[8],'r.')
or
>>> plot(jfa[6],jfa[8],'b.')
Action-angle coordinates using an orbit-integration-based approximation
The adiabatic and Staeckel approximations used above are good for
stars on close-to-circular orbits, but they break down for more
eccentric orbits (specifically, orbits for which the radial and/or
vertical action is of a similar magnitude as the angular
momentum). This is because the approximations made to the potential in
these methods (that it is separable in R and z for the adiabatic
approximation and that it is close to a Staeckel potential for the
Staeckel approximation) break down for such orbits. Unfortunately,
these methods cannot be refined to provide better approximations for
eccentric orbits.
galpy contains a new method for calculating actions, frequencies, and
angles that is completely general for any static potential. It can
calculate the actions to any desired precision for any orbit in such
potentials. The method works by employing an auxiliary isochrone
potential and calculates action-angle variables by arithmetic
operations on the actions and angles calculated in the auxiliary
potential along an orbit (integrated in the true potential). Full
details can be found in Appendix A of Bovy (2014).
We setup this method for a logarithmic potential as follows
>>> from galpy.actionAngle import actionAngleIsochroneApprox
>>> from galpy.potential import LogarithmicHaloPotential
>>> lp= LogarithmicHaloPotential(normalize=1.,q=0.9)
>>> aAIA= actionAngleIsochroneApprox(pot=lp,b=0.8)
b=0.8 here sets the scale parameter of the auxiliary isochrone
potential; this potential can also be specified as an
IsochronePotential instance through ip=). We can now calculate the
actions for an orbit similar to that of the GD-1 stream
>>> obs= numpy.array([1.56148083,0.35081535,-1.15481504,0.88719443,-0.47713334,0.12019596]) #orbit similar to GD-1
>>> aAIA(*obs)
(array([ 0.16605011]), array([-1.80322155]), array([ 0.50704439]))
An essential requirement of this method is that the angles calculated
in the auxiliary potential go through the full range
\([0,2\pi]\). If this is not the case, galpy will raise a warning
>>> aAIA= actionAngleIsochroneApprox(pot=lp,b=10.8)
>>> aAIA(*obs)
galpyWarning: Full radial angle range not covered for at least one object; actions are likely not reliable
(array([ 0.08985167]), array([-1.80322155]), array([ 0.50849276]))
Therefore, some care should be taken to choosing a good auxiliary
potential. galpy contains a method to estimate a decent scale
parameter for the auxiliary scale parameter, which works similar to
estimateDeltaStaeckel above except that it also gives a minimum
and maximum b if multiple R and z are given
>>> from galpy.actionAngle import estimateBIsochrone
>>> from galpy.orbit import Orbit
>>> o= Orbit(obs)
>>> ts= numpy.linspace(0.,100.,1001)
>>> o.integrate(ts,lp)
>>> estimateBIsochrone(o.R(ts),o.z(ts),pot=lp)
(0.78065062339131952, 1.2265541473461612, 1.4899326335155412) #bmin,bmedian,bmax over the orbit
Experience shows that a scale parameter somewhere in the range
returned by this function makes sure that the angles go through the
full \([0,2\pi]\) range. However, even if the angles go through
the full range, the closer the angles increase to linear, the better
the converenge of the algorithm is (and especially, the more accurate
the calculation of the frequencies and angles is, see below). For
example, for the scale parameter at the upper and of the range
>>> aAIA= actionAngleIsochroneApprox(pot=lp,b=1.5)
>>> aAIA(*obs)
(array([ 0.01120145]), array([-1.80322155]), array([ 0.50788893]))
which does not agree with the previous calculation. We can inspect how
the angles increase and how the actions converge by using the
aAIA.plot function. For example, we can plot the radial versus the
vertical angle in the auxiliary potential
>>> aAIA.plot(*obs,type='araz')
which gives
and this clearly shows that the angles increase very non-linearly,
because the auxiliary isochrone potential used is too far from the
real potential. This causes the actions to converge only very
slowly. For example, for the radial action we can plot the converge as a function of integration time
>>> aAIA.plot(*obs,type='jr')
which gives
This Figure clearly shows that the radial action has not converged
yet. We need to integrate much longer in this auxiliary potential to
obtain convergence and because the angles increase so non-linearly, we also need to integrate the orbit much more finely:
>>> aAIA= actionAngleIsochroneApprox(pot=lp,b=1.5,tintJ=1000,ntintJ=800000)
>>> aAIA(*obs)
(array([ 0.01711635]), array([-1.80322155]), array([ 0.51008058]))
>>> aAIA.plot(*obs,type='jr')
which shows slow convergence
Finding a better auxiliary potential makes convergence much faster
and also allows the frequencies and the angles to be calculated by
removing the small wiggles in the auxiliary angles vs. time (in the
angle plot above, the wiggles are much larger, such that removing them
is hard). The auxiliary potential used above had b=0.8, which
shows very quick converenge and good behavior of the angles
>>> aAIA= actionAngleIsochroneApprox(pot=lp,b=0.8)
>>> aAIA.plot(*obs,type='jr')
gives
and
>>> aAIA.plot(*obs,type='araz')
gives
We can remove the periodic behavior from the angles, which clearly
shows that they increase close-to-linear with time
>>> aAIA.plot(*obs,type='araz',deperiod=True)
We can then calculate the frequencies and the angles for this orbit as
>>> aAIA.actionsFreqsAngles(*obs)
(array([ 0.16392384]),
array([-1.80322155]),
array([ 0.50999882]),
array([ 0.55808933]),
array([-0.38475753]),
array([ 0.42199713]),
array([ 0.18739688]),
array([ 0.3131815]),
array([ 2.18425661]))
This function takes as an argument maxn= the maximum n for which
to remove sinusoidal wiggles. So we can raise this, for example to 4
from 3
>>> aAIA.actionsFreqsAngles(*obs,maxn=4)
(array([ 0.16392384]),
array([-1.80322155]),
array([ 0.50999882]),
array([ 0.55808776]),
array([-0.38475733]),
array([ 0.4219968]),
array([ 0.18732009]),
array([ 0.31318534]),
array([ 2.18421296]))
Clearly, there is very little change, as most of the wiggles are of
low n.
Warning
While the orbit-based actionAngle technique in principle works for triaxial potentials, angles and frequencies for non-axisymmetric potentials are not implemented yet.
This technique also works for triaxial potentials, but using those
requires the code to also use the azimuthal angle variable in the
auxiliary potential (this is unnecessary in axisymmetric potentials as
the z component of the angular momentum is conserved). We can
calculate actions for triaxial potentials by specifying that
nonaxi=True:
>>> aAIA(*obs,nonaxi=True)
(array([ 0.16605011]), array([-1.80322155]), array([ 0.50704439]))
galpy currently does not contain any triaxial potentials, so we cannot
illustrate this here with any real triaxial potentials.
Accessing action-angle coordinates for Orbit instances
While the recommended way to access the actionAngle routines is
through the methods in the galpy.actionAngle modules, action-angle
coordinates can also be cacluated for galpy.orbit.Orbit
instances. This is illustrated here briefly. We initialize an Orbit
instance
>>> from galpy.orbit import Orbit
>>> from galpy.potential import MWPotential2014
>>> o= Orbit([1.,0.1,1.1,0.,0.25,0.])
and we can then calculate the actions (default is to use the adiabatic
approximation)
>>> o.jr(MWPotential2014), o.jp(MWPotential2014), o.jz(MWPotential2014)
(0.01685643005901713, 1.1, 0.015897730620467752)
o.jp here gives the azimuthal action (which is the z component
of the angular momentum for axisymmetric potentials). We can also use
the other methods described above, but note that these require extra
parameters related to the approximation to be specified (see above):
>>> o.jr(MWPotential2014,type='staeckel',delta=0.4), o.jp(MWPotential2014,type='staeckel',delta=0.4), o.jz(MWPotential2014,type='staeckel',delta=0.4)
(array([ 0.01922167]), array([ 1.1]), array([ 0.01527683]))
>>> o.jr(MWPotential2014,type='isochroneApprox',b=0.8), o.jp(MWPotential2014,type='isochroneApprox',b=0.8), o.jz(MWPotential2014,type='isochroneApprox',b=0.8)
(array([ 0.01906609]), array([ 1.1]), array([ 0.01528049]))
These two methods give very precise actions for this orbit (both are
converged to about 1%) and they agree very well
>>> (o.jr(MWPotential2014,type='staeckel',delta=0.4)-o.jr(MWPotential2014,type='isochroneApprox',b=0.8))/o.jr(MWPotential2014,type='isochroneApprox',b=0.8)
array([ 0.00816012])
>>> (o.jz(MWPotential2014,type='staeckel',delta=0.4)-o.jz(MWPotential2014,type='isochroneApprox',b=0.8))/o.jz(MWPotential2014,type='isochroneApprox',b=0.8)
array([-0.00024])
Warning
Once an action, frequency, or angle is calculated for a given type of calculation (e.g., staeckel), the parameters for that type are fixed in the Orbit instance. Call o.resetaA() to reset the action-angle instance used when using different parameters (i.e., different delta= for staeckel or different b= for isochroneApprox.
We can also calculate the frequencies and the angles. This requires
using the Staeckel or Isochrone approximations, because frequencies
and angles are currently not supported for the adiabatic
approximation. For example, the radial frequency
>>> o.Or(MWPotential2014,type='staeckel',delta=0.4)
1.1131779637307115
>>> o.Or(MWPotential2014,type='isochroneApprox',b=0.8)
1.1134635974560649
and the radial angle
>>> o.wr(MWPotential2014,type='staeckel',delta=0.4)
0.37758086786371969
>>> o.wr(MWPotential2014,type='isochroneApprox',b=0.8)
0.38159809018175395
which again agree to 1%. We can also calculate the other frequencies,
angles, as well as periods using the functions o.Op, o.oz,
o.wp, o.wz, o.Tr, o.Tp, o.Tz.
Example: Evidence for a Lindblad resonance in the Solar neighborhood
We can use galpy to calculate action-angle coordinates for a set of
stars in the Solar neighborhood and look for unexplained features. For
this we download the data from the Geneva-Copenhagen Survey
(2009A&A...501..941H; data available
at viZier). Since
the velocities in this catalog are given as U,V, and W, we use the
radec and UVW keywords to initialize the orbits from the raw
data. For each object ii
>>> o= Orbit(vxvv[ii,:],radec=True,uvw=True,vo=220.,ro=8.)
We then calculate the actions and angles for each object in a flat
rotation curve potential
>>> lp= LogarithmicHaloPotential(normalize=1.)
>>> myjr[ii]= o.jr(lp)
etc.
Plotting the radial action versus the angular momentum
>>> plot.bovy_plot(myjp,myjr,'k.',ms=2.,xlabel=r'$J_{\phi}$',ylabel=r'$J_R$',xrange=[0.7,1.3],yrange=[0.,0.05])
shows a feature in the distribution
If instead we use a power-law rotation curve with power-law index 1
>>> pp= PowerSphericalPotential(normalize=1.,alpha=-2.)
>>> myjr[ii]= o.jr(pp)
We find that the distribution is stretched, but the feature remains
Code for this example can be found here (note that this code uses a particular
download of the GCS data set; if you use your own version, you will
need to modify the part of the code that reads the data). For more
information see 2010MNRAS.409..145S.
Three-dimensional disk distribution functions
galpy contains a fully three-dimensional disk distribution:
galpy.df.quasiisothermaldf, which is an approximately isothermal
distribution function expressed in terms of action–angle variables
(see 2010MNRAS.401.2318B and
2011MNRAS.413.1889B). Recent
research shows that this distribution function provides a good model
for the DF of mono-abundance sub-populations (MAPs) of the Milky Way
disk (see 2013MNRAS.434..652T and
2013ApJ...779..115B). This
distribution function family requires action-angle coordinates to
evaluate the DF, so galpy.df.quasiisothermaldf makes heavy use of
the routines in galpy.actionAngle (in particular those in
galpy.actionAngleAdiabatic and
galpy.actionAngle.actionAngleStaeckel).
Setting up the DF and basic properties
The quasi-isothermal DF is defined by a gravitational potential and a
set of parameters describing the radial surface-density profile and
the radial and vertical velocity dispersion as a function of
radius. In addition, we have to provide an instance of a
galpy.actionAngle class to calculate the actions for a given
position and velocity. For example, for a
galpy.potential.MWPotential2014 potential using the adiabatic
approximation for the actions, we import and define the following
>>> from galpy.potential import MWPotential2014
>>> from galpy.actionAngle import actionAngleAdiabatic
>>> from galpy.df import quasiisothermaldf
>>> aA= actionAngleAdiabatic(pot=MWPotential2014,c=True)
and then setup the quasiisothermaldf instance
>>> qdf= quasiisothermaldf(1./3.,0.2,0.1,1.,1.,pot=MWPotential2014,aA=aA,cutcounter=True)
which sets up a DF instance with a radial scale length of
\(R_0/3\), a local radial and vertical velocity dispersion of
\(0.2\,V_c(R_0)\) and \(0.1\,V_c(R_0)\), respectively, and a
radial scale lengths of the velocity dispersions of
\(R_0\). cutcounter=True specifies that counter-rotating stars
are explicitly excluded (normally these are just exponentially
suppressed). As for the two-dimensional disk DFs, these parameters are
merely input (or target) parameters; the true density and velocity
dispersion profiles calculated by evaluating the relevant moments of
the DF (see below) are not exactly exponential and have scale lengths
and local normalizations that deviate slightly from these input
parameters. We can estimate the DF’s actual radial scale length near
\(R_0\) as
>>> qdf.estimate_hr(1.)
0.32908034635647182
which is quite close to the input value of 1/3. Similarly, we can
estimate the scale lengths of the dispersions
>>> qdf.estimate_hsr(1.)
1.1913935820372923
>>> qdf.estimate_hsz(1.)
1.0506918075359255
The vertical profile is fully specified by the velocity dispersions
and radial density / dispersion profiles under the assumption of
dynamical equilibrium. We can estimate the scale height of this DF at
a given radius and height as follows
>>> qdf.estimate_hz(1.,0.125)
0.021389597757156088
Near the mid-plane this vertical scale height becomes very large
because the vertical profile flattens, e.g.,
>>> qdf.estimate_hz(1.,0.125/100.)
1.006386030587223
or even
>>> qdf.estimate_hz(1.,0.)
187649.98447377066
which is basically infinity.
Evaluating moments
We can evaluate various moments of the DF giving the density, mean
velocities, and velocity dispersions. For example, the mean radial
velocity is again everywhere zero because the potential and the DF are
axisymmetric
>>> qdf.meanvR(1.,0.)
0.0
Likewise, the mean vertical velocity is everywhere zero
>>> qdf.meanvz(1.,0.)
0.0
The mean rotational velocity has a more interesting dependence on
position. Near the plane, this is the same as that calculated for a similar two-dimensional disk DF (see Evaluating moments of the DF)
>>> qdf.meanvT(1.,0.)
0.91988346380781227
However, this value decreases as one moves further from the plane. The
quasiisothermaldf allows us to calculate the average rotational
velocity as a function of height above the plane. For example,
>>> zs= numpy.linspace(0.,0.25,21)
>>> mvts= numpy.array([qdf.meanvT(1.,z) for z in zs])
which gives
We can also calculate the second moments of the DF. We can check
whether the radial and velocity dispersions at \(R_0\) are close
to their input values
>>> numpy.sqrt(qdf.sigmaR2(1.,0.))
0.20807112565801389
>>> numpy.sqrt(qdf.sigmaz2(1.,0.))
0.090453510526130904
and they are pretty close. We can also calculate the mixed R and z
moment, for example,
>>> qdf.sigmaRz(1.,0.125)
0.0
or expressed as an angle (the tilt of the velocity ellipsoid)
>>> qdf.tilt(1.,0.125)
0.0
This tilt is zero because we are using the adiabatic approximation. As
this approximation assumes that the motions in the plane are decoupled
from the vertical motions of stars, the mixed moment is zero. However,
this approximation is invalid for stars that go far above the
plane. By using the Staeckel approximation to calculate the actions,
we can model this coupling better. Setting up a quasiisothermaldf
instance with the Staeckel approximation
>>> from galpy.actionAngle import actionAngleStaeckel
>>> aAS= actionAngleStaeckel(pot=MWPotential2014,delta=0.45,c=True)
>>> qdfS= quasiisothermaldf(1./3.,0.2,0.1,1.,1.,pot=MWPotential2014,aA=aAS,cutcounter=True)
we can similarly calculate the tilt
>>> qdfS.tilt(1.,0.125)
5.9096430410862419
or about 5 degrees. As a function of height, we find
>>> tilts= numpy.array([qdfS.tilt(1.,z) for z in zs])
>>> plot(zs,tilts)
which gives
We can also calculate the density and surface density (the zero-th
velocity moments). For example, the vertical density
>>> densz= numpy.array([qdf.density(1.,z) for z in zs])
and
>>> denszS= numpy.array([qdfS.density(1.,z) for z in zs])
We can compare the vertical profiles calculated using the adiabatic
and Staeckel action-angle approximations
>>> semilogy(zs,densz/densz[0])
>>> semilogy(zs,denszS/denszS[0])
which gives
Similarly, we can calculate the radial profile of the surface density
>>> rs= numpy.linspace(0.5,1.5,21)
>>> surfr= numpy.array([qdf.surfacemass_z(r) for r in rs])
>>> surfrS= numpy.array([qdfS.surfacemass_z(r) for r in rs])
and compare them with each other and an exponential with scale length
1/3
>>> semilogy(rs,surfr/surfr[10])
>>> semilogy(rs,surfrS/surfrS[10])
>>> semilogy(rs,numpy.exp(-(rs-1.)/(1./3.)))
which gives
The two radial profiles are almost indistinguishable and are very
close, if somewhat shallower, than the pure exponential profile.
General velocity moments, including all higher order moments, are
implemented in quasiisothermaldf.vmomentdensity.
Evaluating and sampling the full probability distribution function
We can evaluate the distribution itself by calling the object, e.g.,
>>> qdf(1.,0.1,1.1,0.1,0.) #input: R,vR,vT,z,vz
array([ 16.86790643])
or as a function of rotational velocity, for example in the mid-plane
>>> vts= numpy.linspace(0.,1.5,101)
>>> pvt= numpy.array([qdfS(1.,0.,vt,0.,0.) for vt in vts])
>>> plot(vts,pvt/numpy.sum(pvt)/(vts[1]-vts[0]))
which gives
This is, however, not the true distribution of rotational velocities
at R =0 and z =0, because it is conditioned on zero radial and
vertical velocities. We can calculate the distribution of rotational
velocities marginalized over the radial and vertical velocities as
>>> qdfS.pvT(1.,1.,0.) #input vT,R,z
14.677231196899195
or as a function of rotational velocity
>>> pvt= numpy.array([qdfS.pvT(vt,1.,0.) for vt in vts])
overplotting this over the previous distribution gives
>>> plot(vts,pvt/numpy.sum(pvt)/(vts[1]-vts[0]))
which is slightly different from the conditioned
distribution. Similarly, we can calculate marginalized velocity
probabilities `pvR, pvz, pvRvT, pvRvz, and
pvTvz. These are all multiplied with the density, such that
marginalizing these over the remaining velocity component results in
the density.
We can sample velocities at a given location using
quasiisothermaldf.sampleV (there is currently no support for
sampling locations from the density profile, although that is rather
trivial):
>>> vs= qdfS.sampleV(1.,0.,n=10000)
>>> hist(vs[:,1],normed=True,histtype='step',bins=101,range=[0.,1.5])
gives
which shows very good agreement with the green (marginalized over vR
and vz) curve (as it should).
