Add dynamical friction

.coveragerc

@@ -13,8 +13,7 @@ exclude_lines =
     if __name__ == .__main__.:
 
 omit = 
-    galpy/snapshot_src*
-    galpy/snapshot.py
+    galpy/snapshot/*
     galpy/util/bovy_plot.py
     galpy/util/bovy_ars.py
     galpy/util/multi.py

.coveragerc_travis

@@ -13,8 +13,7 @@ exclude_lines =
     if __name__ == .__main__.:
 
 omit = 
-    galpy/snapshot_src*
-    galpy/snapshot.py
+    galpy/snapshot/*
     galpy/util/bovy_plot.py
     galpy/util/bovy_ars.py
     galpy/util/multi.py

.gitignore

@@ -41,7 +41,7 @@ coverage_html_report
 
 # galpy #
 #########
-galpy/df_src/data/dfcorrection*.sav
+galpy/df/data/dfcorrection*.sav
 
 
 # Build directories #

.travis.yml

@@ -10,10 +10,10 @@ env: #split tests
     - REQUIRES_ASTROQUERY=false
     - PYTHON_COVREPORTS_VERSION=3.6 # Version for which reports are uploaded
   matrix:
-    - TEST_FILES='tests/ --ignore=tests/test_qdf.py --ignore=tests/test_pv2qdf.py --ignore=tests/test_diskdf.py --ignore=tests/test_orbit.py --ignore=tests/test_streamdf.py --ignore=tests/test_streamgapdf.py --ignore=tests/test_evolveddiskdf.py --ignore=tests/test_quantity.py --ignore=tests/test_nemo.py --ignore=tests/test_coords.py' REQUIRES_PYNBODY=true
+    - TEST_FILES='tests/ --ignore=tests/test_qdf.py --ignore=tests/test_pv2qdf.py --ignore=tests/test_diskdf.py --ignore=tests/test_orbit.py --ignore=tests/test_streamdf.py --ignore=tests/test_streamgapdf.py --ignore=tests/test_evolveddiskdf.py --ignore=tests/test_quantity.py --ignore=tests/test_nemo.py --ignore=tests/test_coords.py --ignore=tests/test_jeans.py' REQUIRES_PYNBODY=true
     - TEST_FILES='tests/test_quantity.py tests/test_coords.py' REQUIRES_ASTROPY=true # needs to be separate for different config
     - TEST_FILES='tests/test_orbit.py' REQUIRES_PYNBODY=true REQUIRES_ASTROPY=true REQUIRES_ASTROQUERY=true
-    - TEST_FILES='tests/test_evolveddiskdf.py'
+    - TEST_FILES='tests/test_evolveddiskdf.py tests/test_jeans.py'
     - TEST_FILES='tests/test_diskdf.py'
     - TEST_FILES='tests/test_qdf.py tests/test_pv2qdf.py tests/test_streamgapdf.py'
     - TEST_FILES='tests/test_streamdf.py'
@@ -60,7 +60,7 @@ before_install:
  - sh -e /etc/init.d/xvfb start $For plotting tests
 #Download corrections for some tests
  - curl -O https://github.s3.amazonaws.com/downloads/jobovy/galpy/galpy-dfcorrections.tar.gz
- - tar xvzf galpy-dfcorrections.tar.gz -C ./galpy/df_src/data/
+ - tar xvzf galpy-dfcorrections.tar.gz -C ./galpy/df/data/
 # command to install dependencies
 install:
  - pip install 'coverage==4.1' --force-reinstall # necessary bc of backward incompatible change in 4.2 about combining reports
@@ -71,8 +71,8 @@ install:
  - if [[ $TRAVIS_PYTHON_VERSION == $PYTHON_COVREPORTS_VERSION ]]; then easy_install --upgrade coveralls; fi
  - if $REQUIRES_PYNBODY; then pip install git+git://github.com/pynbody/pynbody.git; fi
  # clone my version of the torus code, don't do this for one test, to make sure the code installs without the torus code
- - if [[ $TEST_FILES != 'tests/test_evolveddiskdf.py' ]]; then git clone https://github.com/jobovy/Torus.git galpy/actionAngle_src/actionAngleTorus_c_ext/torus; fi
- - if [[ $TEST_FILES != 'tests/test_evolveddiskdf.py' ]]; then cd galpy/actionAngle_src/actionAngleTorus_c_ext/torus; fi
+ - if [[ $TEST_FILES != 'tests/test_evolveddiskdf.py' ]]; then git clone https://github.com/jobovy/Torus.git galpy/actionAngle/actionAngleTorus_c_ext/torus; fi
+ - if [[ $TEST_FILES != 'tests/test_evolveddiskdf.py' ]]; then cd galpy/actionAngle/actionAngleTorus_c_ext/torus; fi
  - if [[ $TEST_FILES != 'tests/test_evolveddiskdf.py' ]]; then git checkout galpy; fi
  - if [[ $TEST_FILES != 'tests/test_evolveddiskdf.py' ]]; then cd -; fi
  - if $REQUIRES_ASTROPY; then conda install astropy; fi
@@ -88,7 +88,7 @@ script:
 after_success:
  # Generate lcov output 
  - if [[ $TRAVIS_PYTHON_VERSION == $PYTHON_COVREPORTS_VERSION ]]; then lcov --capture --base-directory . --directory build/temp.linux-x86_64-3.6/galpy/ --no-external --output-file coverage_full.info; fi
- - if [[ $TRAVIS_PYTHON_VERSION == $PYTHON_COVREPORTS_VERSION ]]; then lcov --remove coverage_full.info 'galpy/actionAngle_src/actionAngleTorus_c_ext/torus/*' -o coverage.info; fi
+ - if [[ $TRAVIS_PYTHON_VERSION == $PYTHON_COVREPORTS_VERSION ]]; then lcov --remove coverage_full.info 'galpy/actionAngle/actionAngleTorus_c_ext/torus/*' -o coverage.info; fi
  # Codecov
  - if [[ $TRAVIS_PYTHON_VERSION == $PYTHON_COVREPORTS_VERSION ]]; then bash <(curl -s https://codecov.io/bash) -X gcov; fi
  # coveralls: combine, generate json, and upload

HISTORY.txt

@@ -1,6 +1,17 @@
 v1.4 (??????)
 ==============
 
+- Added ChandrasekharDynamicalFrictionForce, an implementation of
+  dynamical friction based on the classical Chandrasekhar formula
+  (with recent tweaks from the literature to better represent the
+  results from N-body simulations).
+
+- Added galpy.df.jeans with tools for Jeans modeling. Currently only
+  contains the function sigmar to calculate the velocity dispersion
+  using the spherical Jeans equation in a given potential, density
+  profile, and anisotropy profile (anisotropy can be radially
+  varying).
+
 - Added CorotatingRotationWrapperPotential to galpy.potential: a
   wrapper to make a pattern wind up over time such that it is always
   corotating (for use in simulating corotating spiral structure like
@@ -16,6 +27,9 @@ v1.4 (??????)
   components (e.g., a bar) to previously defined potentials (which may
   be lists themselves): new_pot= [pot,bar_pot].
 
+- Add from_name class method of Orbit, which allows Orbit instances to
+  be initialized using the coordinates of a named object found in SIMBAD.
+
 - Defined Orbit initialization without any arguments (or, more
   generally, without specifying the vxvv initial phase-space input) to
   return an Orbit instances representing the Sun. Can therefore setup
@@ -48,12 +62,22 @@ v1.4 (??????)
 - Tweaked coordinate-transformations to Galactocentric coordinates to
   be consistent with astropy's.
 
+- Introduced general Force class of which Potential and
+  DissipativeForce inherit, the former for forces that derive from a
+  potential, the latter for those that do not.
+
+- Introduced DissipativeForce, a superclass for all dissipative
+  forces. ChandrasekharDynamicalFrictionForce is currently the only
+  class that inherits from DissipativeForce.
+
+- Re-arranged the package structure to better comply with the standard
+  layout. All subpackages (e.g., galpy.potential) are now contained in
+  subdirectories of the same name (e.g., galpy/potential/ rather than
+  the old galpy/potential_src/).
+
 - Made the code fully compilable on Windows with MSVC and test Windows
   builds automatically on appveyor (#333).
 
-- Add from_name class method of Orbit, which allows Orbit instances to
-  be initialized using the coordinates of a named object found in SIMBAD.
-
 v1.3 (2018-02-06)
 ==================
 

MANIFEST.in

@@ -1,9 +1,9 @@
 include README.rst README.dev README.nemo LICENSE HISTORY.txt AUTHORS.txt
 include gsl-config.bat
-include galpy/df_src/data/*.sav
-include galpy/actionAngle_src/actionAngle_c_ext/*.h
-include galpy/actionAngle_src/actionAngleTorus_c_ext/*.h
-include galpy/orbit_src/orbit_c_ext/*.h
-include galpy/potential_src/potential_c_ext/*.h
+include galpy/df/data/*.sav
+include galpy/actionAngle/actionAngle_c_ext/*.h
+include galpy/actionAngle/actionAngleTorus_c_ext/*.h
+include galpy/orbit/orbit_c_ext/*.h
+include galpy/potential/potential_c_ext/*.h
 include galpy/util/*.h
 include galpy/util/interp_2d/*.h

README.dev

@@ -5,26 +5,26 @@ Adding a potential to the C integrator
 --------------------------------------
 
 1) Implement the potential in a .c file under
-potential_src/potential_c_ext. Look at
-potential_src/potential_c_ext/LogarithmicHaloPotential.c for the right
+potential/potential_c_ext. Look at
+potential/potential_c_ext/LogarithmicHaloPotential.c for the right
 format
 
 2) Add your new potential to
-potential_src/potential_c_ext/galpy_potentials.h
+potential/potential_c_ext/galpy_potentials.h
 
-3) Edit the code under orbit_src/orbit_c_ext/integratePlanarOrbit.c to
+3) Edit the code under orbit/orbit_c_ext/integratePlanarOrbit.c to
 set up your new potential (in the 'parse_leapFuncArgs' function)
 
-4) Edit the code in orbit_src/integratePlanarOrbit.py to set up your
+4) Edit the code in orbit/integratePlanarOrbit.py to set up your
 new potential
 
-5) Edit the code under orbit_src/orbit_c_ext/integrateFullOrbit.c to
+5) Edit the code under orbit/orbit_c_ext/integrateFullOrbit.c to
 set up your new potential (in the 'parse_leapFuncArgs_Full' function)
 
-6) Edit the code in orbit_src/integrateFullOrbit.py to set up your
+6) Edit the code in orbit/integrateFullOrbit.py to set up your
 new potential
 
-7) Edit the code in actionAngle_src/actionAngle_c_ext/actionAngle.c to
+7) Edit the code in actionAngle/actionAngle_c_ext/actionAngle.c to
 parse the new potential
 
 8) Finally, add 'self.hasC= True' to the initialization of the

README.rst

@@ -143,7 +143,7 @@ profiles more closely (see `1999AJ....118.1201D
 these corrections is expensive, and a large set of precalculated
 corrections can be found `here
 <http://github.com/downloads/jobovy/galpy/galpy-dfcorrections.tar.gz>`__
-\[tar.gz archive\]. Install these by downloading them and unpacking them into the galpy/df_src/data directory before running the setup.py installation. E.g.::
+\[tar.gz archive\]. Install these by downloading them and unpacking them into the galpy/df/data directory before running the setup.py installation. E.g.::
 
    curl -O https://github.s3.amazonaws.com/downloads/jobovy/galpy/galpy-dfcorrections.tar.gz
-   tar xvzf galpy-dfcorrections.tar.gz -C ./galpy/df_src/data/
+   tar xvzf galpy-dfcorrections.tar.gz -C ./galpy/df/data/

doc/source/installation.rst

@@ -100,8 +100,8 @@ To install the TorusMapper code, *before* running the installation of
 galpy, navigate to the top-level galpy directory (which contains the
 ``setup.py`` file) and do::
 
-	     git clone https://github.com/jobovy/Torus.git galpy/actionAngle_src/actionAngleTorus_c_ext/torus
-	     cd galpy/actionAngle_src/actionAngleTorus_c_ext/torus
+	     git clone https://github.com/jobovy/Torus.git galpy/actionAngle/actionAngleTorus_c_ext/torus
+	     cd galpy/actionAngle/actionAngleTorus_c_ext/torus
 	     git checkout galpy
 	     cd -
 

doc/source/orbit.rst

@@ -192,6 +192,8 @@ natural coordinates, you can turn this behavior off by doing
 
 All outputs will then be specified in galpy's natural coordinates.
 
+.. _orbfromname:
+
 Initialization from an object's name
 ****************************************
 
@@ -735,3 +737,138 @@ the estimations in one batch using the ``actionAngle`` interface, which takes co
 >>> par = aAS.EccZmaxRperiRap(Rphiz[:,0], vRvTvz[:,0], vRvTvz[:,1], Rphiz[:,2], vRvTvz[:,2], Rphiz[:,1], delta=deltas)
 
 The above code calculates the parameters in roughly 100ms on a single core.
+
+**NEW in v1.4** Example: The orbit of the Large Magellanic Cloud in the presence of dynamical friction
+--------------------------------------------------------------------------------------------------------
+
+As a further example of what you can do with galpy, we investigate the
+Large Magellanic Cloud's (LMC) past and future orbit. Because the LMC
+is a massive satellite of the Milky Way, its orbit is affected by
+dynamical friction, a frictional force of gravitational origin that
+occurs when a massive object travels through a sea of low-mass objects
+(halo stars and dark matter in this case). First we import all the
+necessary packages:
+
+>>> from astropy import units
+>>> from galpy.potential import MWPotential2014, ChandrasekharDynamicalFrictionForce
+>>> from galpy.orbit import Orbit
+
+(also do ``%pylab inline`` if running this in a jupyter notebook or
+turn on the ``pylab`` option in ipython for plotting). We can load the
+current phase-space coordinates for the LMC using the
+``Orbit.from_name`` function described :ref:`above <orbfromname>`:
+
+>>> o= Orbit.from_name('LMC')
+
+We will use ``MWPotential2014`` as our Milky-Way potential
+model. Because the LMC is in fact unbound in ``MWPotential2014``, we
+increase the halo mass by 50% to make it bound (this corresponds to a
+Milky-Way halo mass of :math:`\approx 1.2\,\times 10^{12}\,M_\odot`, a
+not unreasonable value). We can hack this together as
+
+>>> MWPotential2014[2]._amp*= 1.5
+
+(Note that this is *not* a generally recommended route for changing
+the mass of an object, since it relies on editing a private
+attribute). Let us now integrate the orbit backwards in time for 10
+Gyr and plot it:
+
+>>> ts= numpy.linspace(0.,-10.,1001)*units.Gyr
+>>> o.integrate(ts,MWPotential2014)
+>>> o.plot(d1='t',d2='r')
+
+.. image:: images/lmc-mwp14.png
+        :scale: 50 %
+
+We see that the LMC is indeed bound, with an apocenter just over 250
+kpc. Now let's add dynamical friction for the LMC, assuming that its
+mass if :math:`5\times 10^{10}\,M_\odot`. We setup the
+dynamical-friction object:
+
+>>> cdf= ChandrasekharDynamicalFrictionForce(GMs=5.*10.**10.*units.Msun,rhm=5.*units.kpc,
+					     dens=MWPotential2014)
+
+Dynamical friction depends on the velocity distribution of the halo,
+which is assumed to be an isotropic Gaussian distribution with a
+radially-dependent velocity dispersion. If the velocity dispersion is
+not given (like in the example above), it is computed from the
+spherical Jeans equation. We have set the half-mass radius to 5 kpc
+for definiteness. We now make a copy of the orbit instance above and
+integrate it in the potential that includes dynamical friction:
+
+>>> odf= o()
+>>> odf.integrate(ts,[MWPotential2014,cdf])
+
+(Note that specifying the forces as the list ``[MWPotential2014,cdf]``
+works even though ``MWPotential2014`` is itself a list of potentials,
+because we can use nested lists of potentials or forces wherever a
+list is allowed in ``galpy``). Overlaying the orbits, we can see the
+difference in the evolution:
+
+>>> o.plot(d1='t',d2='r',label=r'$\mathrm{No\ DF}$')
+>>> odf.plot(d1='t',d2='r',overplot=True,label=r'$\mathrm{DF}, M=5\times10^{10}\,M_\odot$')
+>>> ylim(0.,400.)
+>>> legend()
+
+.. image:: images/lmc-mwp14-plusdynfric-51010msun.png
+        :scale: 50 %
+
+We see that dynamical friction removes energy from the LMC's orbit,
+such that its past apocenter is now around 400 kpc rather than 250
+kpc! The period of the orbit is therefore also much longer. Clearly,
+dynamical friction has a big impact on the orbit of the LMC.
+
+Recent observations have suggested that the LMC may be even more
+massive than what we have assumed so far, with masses over
+:math:`10^{11}\,M_\odot` seeming in good agreement with various
+observations. Let's see how a mass of :math:`10^{11}\,M_\odot` changes
+the past orbit of the LMC. We can change the mass of the LMC used in
+the dynamical-friction calculation as
+
+>>> cdf.GMs= 10.**11.*units.Msun
+
+This way of changing the mass is preferred over re-initializing the
+``ChandrasekharDynamicalFrictionForce`` object, because it avoids
+having to solve the Jeans equation again to obtain the velocity
+dispersion. Then we integrate the orbit and overplot it on the
+previous results:
+
+>>> odf2= o()
+>>> odf2.integrate(ts,[MWPotential2014,cdf])
+
+and
+
+>>> o.plot(d1='t',d2='r',label=r'$\mathrm{No\ DF}$')
+>>> odf.plot(d1='t',d2='r',overplot=True,label=r'$\mathrm{DF}, M=5\times10^{10}\,M_\odot$')
+>>> odf2.plot(d1='t',d2='r',overplot=True,label=r'$\mathrm{DF}, M=1\times10^{11}\,M_\odot$')
+>>> ylim(0.,740.)
+>>> legend()
+
+which gives
+
+.. image:: images/lmc-mwp14-plusdynfric-1011msun.png
+        :scale: 50 %
+
+Now the apocenter increases to about 600 kpc and the LMC doesn't
+perform a full orbit over the last 10 Gyr.
+
+Finally, let's see what will happen in the future if the LMC is as
+massive as :math:`10^{11}\,M_\odot`. We simply flip the sign of the
+integration times to get the future trajectory:
+
+>>> odf2.integrate(-ts[-ts < 9*units.Gyr],[MWPotential2014,cdf])
+>>> odf2.plot(d1='t',d2='r')
+
+.. image:: images/lmc-mwp14-plusdynfric-1011msun-future.png
+   :scale: 50 %
+
+Because of the large effect of dynamical friction, the LMC will merge
+with the Milky-Way in about 4 Gyr after a few more pericenter
+passages. Note that we have not taken any mass-loss into
+account. Because mass-loss would lead to a smaller dynamical-friction
+force, this would somewhat increase the merging timescale, but
+dynamical friction will inevitably lead to the merger of the LMC with
+the Milky Way.
+
+.. WARNING::
+   When using dynamical friction, if the radius gets very small, the integration sometimes becomes very erroneous, which can lead to a big, unphysical kick (even though we turn off friction at very small radii); this is the reason why we have limited the future integration to 9 Gyr in the example above. When using dynamical friction, inspect the full orbit to make sure to catch whether a merger has happened.