Add module with Milky-Way-like potentials

HISTORY.txt

@@ -14,6 +14,11 @@ v1.5 (????-??-??)
   Orbit.from_name also supports tab completion for known objects in
   this list in IPython/Jupyter (PR #392).
 
+- Added the galpy.potentials.mwpotentials module with various
+  Milky-Way-like potentials. Currently included are MWPotential2014,
+  McMillan17 for the potential from McMillan (2017), and the three
+  models from Irrgang et al. (2013). See PR #395.
+
 - Added potential.toVerticalPotential to convert any 3D potential to a
   1D potential at a given (R,phi) [generalizes RZToverticalPotential
   to non-axi potentials; RZToverticalPotential retained for backwards

doc/source/reference/potential.rst

@@ -218,13 +218,19 @@ Dissipative forces
 
 .. _potential-mw:
 
-In addition to these classes, a simple Milky-Way-like potential fit to
-data on the Milky Way is included as
-``galpy.potential.MWPotential2014`` (see the ``galpy`` paper for
-details). Note that this potential assumes a circular velocity of 220
-km/s at the solar radius at 8 kpc; see `arXiv/1412.3451
+Milky-Way-like potentials
+-------------------------
+
+``galpy`` contains various simple models for the Milky Way's
+gravitational potential. The recommended model, described in `Bovy
+(2015) <http://arxiv.org/abs/1412.3451>`_, is included as
+``galpy.potential.MWPotential2014``. This potential was fit to a large
+variety of data on the Milky Way and thus serves as both a simple and
+accurate model for the Milky Way's potential (see `Bovy 2015
 <http://arxiv.org/abs/1412.3451>`_ for full information on how this
-potential was fit. This potential is defined as
+potential was fit). Note that this potential assumes a circular
+velocity of 220 km/s at the solar radius at 8 kpc. This potential is
+defined as
 
 >>> bp= PowerSphericalPotentialwCutoff(alpha=1.8,rc=1.9/8.,normalize=0.05)
 >>> mp= MiyamotoNagaiPotential(a=3./8.,b=0.28/8.,normalize=.6)
@@ -272,16 +278,95 @@ Dehnen's gyrfalcON code using ``accname=PowSphwCut+MiyamotoNagai+NFW``
 and
 ``accpars=0,1001.79126907,1.8,1.9#0,306770.418682,3.0,0.28#0,16.0,162.958241887``.
 
-An older version ``galpy.potential.MWPotential`` of a similar
-potential that was *not* fit to data on the Milky Way is defined as
+``galpy`` also contains other models for the Milky Way's potential
+from the literature in the ``galpy.potential.mwpotentials`` module
+(which also contains ``MWPotential2014``). Currently, these are:
+
+* ``McMillan17``: the potential model from `McMillan (2017) <https://ui.adsabs.harvard.edu/abs/2017MNRAS.465...76M>`_
+* ``Irrgang13I``: model I from `Irrgang et al. (2013) <https://ui.adsabs.harvard.edu/abs/2013A%26A...549A.137I>`_, which is an updated version of the classic `Allen & Santillan (1991) <https://ui.adsabs.harvard.edu/abs/1991RMxAA..22..255A>`_
+* ``Irrgang13II`` and ``Irrgang13III``: model II and III from `Irrgang et al. (2013) <https://ui.adsabs.harvard.edu/abs/2013A%26A...549A.137I>`_
+
+Unlike ``MWPotential2014``, these potentials have physical units
+turned on, using as the unit scaling parameters ``ro`` and ``vo`` the
+distance to the Galactic center and the circular velocity at the Sun's
+radius of each potential. These can be obtained using the
+``galpy.util.bovy_conversion.get_physical`` function, e.g.,
+
+>>> from galpy.potential.mwpotentials import McMillan17
+>>> from galpy.util.bovy_conversion import get_physical
+>>> get_physical(McMillan17)
+# {'ro': 8.21, 'vo': 233.1}
+
+This function returns the unit-conversion parameters as a dictionary,
+so they can be easily passed to other functions. For example, when
+integrating an orbit in these potentials and either initializing the
+orbit using observed coordinates or converting the integrated orbit to
+observed coordinates, it is important to use the same unit-conversion
+parameters (otherwise an error will be raised). For example, to obtain the orbit of the Sun in the ``McMillan17`` potential, we do
+
+>>> from galpy.orbit import Orbit
+>>> o= Orbit(**get_physical(McMillan17))
+
+As an example, we integrate the Sun's orbit for 10 Gyr in
+``MWPotential2014``, ``McMillan17`` and ``Irrgang13I``
+
+>>> from galpy.potential.mwpotentials import MWPotential2014, McMillan17, Irrgang13I
+>>> from galpy.orbit import Orbit
+>>> from galpy.util.bovy_conversion import get_physical
+>>> from astropy import units
+>>> times= numpy.linspace(0.,10.,3001)*units.Gyr
+>>> o_mwp14= Orbit(ro=8.,vo=220.) # Need to set these by hand
+>>> o_mcm17= Orbit(**get_physical(McMillan17))
+>>> o_irrI= Orbit(**get_physical(Irrgang13I))
+>>> o_mwp14.integrate(times,MWPotential2014)
+>>> o_mcm17.integrate(times,McMillan17)
+>>> o_irrI.integrate(times,Irrgang13I)
+>>> o_mwp14.plot(lw=0.6)
+>>> o_mcm17.plot(overplot=True,lw=0.6)
+>>> o_irrI.plot(overplot=True,lw=0.6)
+
+which gives
+
+.. image:: ../images/orbit-sun-mwpotentials.png
+   :scale: 40 %
+
+Much of the difference between these orbits is due to the different
+present Galactocentric radius of the Sun, if we simply plot the
+difference with respect to the present Galactocentric radius, they
+agree better
+
+>>> o_mwp14.plot(d1='R-8.',d2='z',lw=0.6,xlabel=r'$R-R_0\,(\mathrm{kpc})$')
+>>> o_mcm17.plot(d1='R-{}'.format(get_physical(McMillan17)['ro']),d2='z',overplot=True,lw=0.6)
+>>> o_irrI.plot(d1='R-{}'.format(get_physical(Irrgang13I)['ro']),d2='z',overplot=True,lw=0.6)
+
+.. image:: ../images/orbit-sun-mwpotentials-vsRsun.png
+   :scale: 40 %
+
+We can also compare the rotation curves of these different models
+
+>>> from galpy.potential import plotRotcurve
+>>> plotRotcurve(MWPotential2014,label=r'$\mathrm{MWPotential2014}$',ro=8.,vo=220.) # need to set ro and vo explicitly, because MWPotential2014 has units turned off
+>>> plotRotcurve(McMillan17,overplot=True,label=r'$\mathrm{McMillan\, (2017)}$')
+>>> plotRotcurve(Irrgang13I,overplot=True,label=r'$\mathrm{Irrgang\ et\ al.\, (2017), model\ I}$')
+>>> legend()
+
+.. image:: ../images/mwpotentials-vcirc.png
+   :scale: 40 %
+
+
+
+
+An older version ``galpy.potential.MWPotential`` of
+``MWPotential2014`` that was *not* fit to data on the Milky Way is
+defined as
 
 >>> 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)
 >>> MWPotential= mp+np+hp
 
-``galpy.potential.MWPotential2014`` supersedes
-``galpy.potential.MWPotential``.
+but ``galpy.potential.MWPotential2014`` supersedes
+``galpy.potential.MWPotential`` and its use is no longer recommended.
 
 2D potentials
 -------------

galpy/potential/DiskSCFPotential.py

@@ -354,7 +354,7 @@ def _phiforce(self,R,z,phi=0.,t=0.):
         """
         return self._scf.phiforce(R,z,phi=phi,use_physical=False)
 
-    def _R2deriv(self,R,z,phi=0.,t=0.): #pragma: no cover
+    def _R2deriv(self,R,z,phi=0.,t=0.):
         """
         NAME:
            _R2deriv
@@ -370,16 +370,14 @@ def _R2deriv(self,R,z,phi=0.,t=0.): #pragma: no cover
         HISTORY:
            2016-12-26 - Written - Bovy (UofT/CCA)
         """
-        raise AttributeError
-        # Implementation above does not work bc SCF.R2deriv is not implemented
         r= numpy.sqrt(R**2.+z**2.)
         out= self._scf.R2deriv(R,z,phi=phi,use_physical=False)
         for a,ds,d2s,H in zip(self._Sigma_amp,self._dSigmadR,self._d2SigmadR2,
                               self._Hz):
             out+= 4.*numpy.pi*a*H(z)/r**2.*(d2s(r)*R**2.+z**2./r*ds(r))
         return out
 
-    def _z2deriv(self,R,z,phi=0.,t=0.): #pragma: no cover
+    def _z2deriv(self,R,z,phi=0.,t=0.):
         """
         NAME:
            _z2deriv
@@ -395,8 +393,6 @@ def _z2deriv(self,R,z,phi=0.,t=0.): #pragma: no cover
         HISTORY:
            2016-12-26 - Written - Bovy (UofT/CCA)
         """
-        raise AttributeError
-        # Implementation above does not work bc SCF.z2deriv is not implemented
         r= numpy.sqrt(R**2.+z**2.)
         out= self._scf.z2deriv(R,z,phi=phi,use_physical=False)
         for a,s,ds,d2s,h,H,dH in zip(self._Sigma_amp,
@@ -406,7 +402,7 @@ def _z2deriv(self,R,z,phi=0.,t=0.): #pragma: no cover
                                  +2.*ds(r)*dH(z)*z/r+s(r)*h(z))
         return out
 
-    def _Rzderiv(self,R,z,phi=0.,t=0.): #pragma: no cover
+    def _Rzderiv(self,R,z,phi=0.,t=0.):
         """
         NAME:
            _Rzderiv
@@ -422,17 +418,15 @@ def _Rzderiv(self,R,z,phi=0.,t=0.): #pragma: no cover
         HISTORY:
            2016-12-26 - Written - Bovy (UofT/CCA)
         """
-        raise AttributeError
-        # Implementation above does not work bc SCF.Rzderiv is not implemented
         r= numpy.sqrt(R**2.+z**2.)
         out= self._scf.Rzderiv(R,z,phi=phi,use_physical=False)
-        for a,ds,d2s,H,dH in zip(self._Sigma_amp,self._dsigmadR,
+        for a,ds,d2s,H,dH in zip(self._Sigma_amp,self._dSigmadR,
                                  self._d2SigmadR2,self._Hz,self._dHzdz):
             out+= 4.*numpy.pi*a*(H(z)*R*z/r**2.*(d2s(r)-ds(r)/r)
                                  +ds(r)*dH(z)*R/r)
         return out
 
-    def _phi2deriv(self,R,z,phi=0.,t=0.): #pragma: no cover
+    def _phi2deriv(self,R,z,phi=0.,t=0.):
         """
         NAME:
            _phi2deriv
@@ -448,8 +442,6 @@ def _phi2deriv(self,R,z,phi=0.,t=0.): #pragma: no cover
         HISTORY:
            2016-12-26 - Written - Bovy (UofT/CCA)
         """
-        raise AttributeError
-        # Implementation above does not work bc SCF.phi2deriv is not implemented
         return self._scf.phi2deriv(R,z,phi=phi,use_physical=False)
 
     def _dens(self,R,z,phi=0.,t=0.):

galpy/potential/Irrgang13.py

@@ -0,0 +1,77 @@
+# Milky-Way mass models from Irrgang et al. (2013)
+import numpy
+from ..potential import PlummerPotential
+from ..potential import MiyamotoNagaiPotential
+from ..potential import SCFPotential
+from ..potential import NFWPotential
+from ..potential import scf_compute_coeffs_spherical
+from ..util import bovy_conversion
+# Their mass unit
+mgal_in_msun= 1e5/bovy_conversion._G
+# Model I: updated version of Allen & Santillan
+# Unit normalizations
+ro, vo= 8.4, 242.
+Irrgang13I_bulge= PlummerPotential(\
+    amp=409.*mgal_in_msun/bovy_conversion.mass_in_msol(vo,ro),
+    b=0.23/ro,ro=ro,vo=vo)
+Irrgang13I_disk= MiyamotoNagaiPotential(\
+    amp=2856.*mgal_in_msun/bovy_conversion.mass_in_msol(vo,ro),
+    a=4.22/ro,b=0.292/ro,ro=ro,vo=vo)
+# The halo is a little more difficult, because the Irrgang13I halo model is 
+# not in galpy, so we use SCF to represent it (because we're lazy...). The 
+# sharp cut-off in the Irrgang13I halo model makes SCF difficult, so we 
+# replace it with a smooth cut-off; this only affects the very outer halo
+def Irrgang13I_halo_dens(\
+    r,amp=1018*mgal_in_msun/bovy_conversion.mass_in_msol(vo,ro),
+    ah=2.562/ro,gamma=2.,Lambda=200./ro):
+    r_over_ah_gamma= (r/ah)**(gamma-1.)
+    return amp/4./numpy.pi/ah*r_over_ah_gamma*(r_over_ah_gamma+gamma)/r**2\
+        /(1.+r_over_ah_gamma)**2.\
+        *((1.-numpy.tanh((r-Lambda)/(Lambda/20.)))/2.)
+a_for_scf= 20.
+# scf_compute_coeffs_spherical currently seems to require a function of 3 parameters...
+Acos= scf_compute_coeffs_spherical(\
+    lambda r,z,p: Irrgang13I_halo_dens(r),40,a=a_for_scf)[0]
+Irrgang13I_halo= SCFPotential(Acos=Acos,a=a_for_scf,ro=ro,vo=vo)
+# Final model I
+Irrgang13I= Irrgang13I_bulge+Irrgang13I_disk+Irrgang13I_halo
+
+# Model II
+# Unit normalizations
+ro, vo= 8.35, 240.4
+Irrgang13II_bulge= PlummerPotential(\
+    amp=175.*mgal_in_msun/bovy_conversion.mass_in_msol(vo,ro),
+    b=0.184/ro,ro=ro,vo=vo)
+Irrgang13II_disk= MiyamotoNagaiPotential(\
+    amp=2829.*mgal_in_msun/bovy_conversion.mass_in_msol(vo,ro),
+    a=4.85/ro,b=0.305/ro,ro=ro,vo=vo)
+# Again use SCF because the Irrgang13II halo model is not in galpy; because 
+# the halo model is quite different from Hernquist both in the inner and outer
+# part, need quite a few basis functions...
+def Irrgang13II_halo_dens(\
+    r,amp=69725*mgal_in_msun/bovy_conversion.mass_in_msol(vo,ro),
+    ah=200./ro):
+    return amp/4./numpy.pi*ah**2./r**2./(r**2.+ah**2.)**1.5
+a_for_scf= 0.15
+# scf_compute_coeffs_spherical currently seems to require a function of 3 parameters...
+Acos= scf_compute_coeffs_spherical(\
+    lambda r,z,p: Irrgang13II_halo_dens(r),75,a=a_for_scf)[0]
+Irrgang13II_halo= SCFPotential(Acos=Acos,a=a_for_scf,ro=ro,vo=vo)
+# Final model II
+Irrgang13II= Irrgang13II_bulge+Irrgang13II_disk+Irrgang13II_halo
+
+# Model III
+# Unit normalizations
+ro, vo= 8.33, 239.7
+Irrgang13III_bulge= PlummerPotential(\
+    amp=439.*mgal_in_msun/bovy_conversion.mass_in_msol(vo,ro),
+    b=0.236/ro,ro=ro,vo=vo)
+Irrgang13III_disk= MiyamotoNagaiPotential(\
+    amp=3096.*mgal_in_msun/bovy_conversion.mass_in_msol(vo,ro),
+    a=3.262/ro,b=0.289/ro,ro=ro,vo=vo)
+Irrgang13III_halo= NFWPotential(\
+    amp=142200.*mgal_in_msun/bovy_conversion.mass_in_msol(vo,ro),
+    a=45.02/ro,ro=ro,vo=vo)
+# Final model III
+Irrgang13III= Irrgang13III_bulge+Irrgang13III_disk+Irrgang13III_halo
+