Merge c8a9194 into 2160d5e

.travis.yml

@@ -11,12 +11,12 @@ env: #split tests
     - REQUIRES_ASTROQUERY=false
     - PYTHON_COVREPORTS_VERSION=3.8 # 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_amuse.py --ignore=tests/test_coords.py --ignore=tests/test_jeans.py --ignore=tests/test_orbits.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_amuse.py --ignore=tests/test_coords.py --ignore=tests/test_jeans.py --ignore=tests/test_orbits.py --ignore=tests/test_dynamfric.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 tests/test_orbits.py' REQUIRES_PYNBODY=true REQUIRES_ASTROPY=true REQUIRES_ASTROQUERY=true
-    - TEST_FILES='tests/test_evolveddiskdf.py tests/test_jeans.py'
-    - TEST_FILES='tests/test_diskdf.py'
+    - TEST_FILES='tests/test_evolveddiskdf.py tests/test_jeans.py tests/test_dynamfric.py'
     - TEST_FILES='tests/test_qdf.py tests/test_pv2qdf.py tests/test_streamgapdf.py'
+    - TEST_FILES='tests/test_diskdf.py'
     - TEST_FILES='tests/test_streamdf.py'
 matrix: # only run crucial tests for python 2.7, 3.6, 3.7
   include:

HISTORY.txt

@@ -1,4 +1,4 @@
-v2.0 (????-??-??)
+v1.6 (????-??-??)
 =================
 
 - Added HomogeneousSpherePotential, the potential of a constant
@@ -8,9 +8,18 @@ v2.0 (????-??-??)
   avoid use of math, don't use scipy's numpy interface, internal
   imports for all galpy imports.
 
+- Implemented ChandrasekharDynamicalFrictionForce in C, introducing a
+  general scheme for easily implementing velocity-dependent forces in
+  C without requiring all forces to take velocity arguments (#420).
+
 - Fixed AMUSE compatibility with Potentials instantiated with physical
   inputs (#405).
 
+- Fix how DiskSCFPotential instances with a hole in the surface
+  density get passed to C (was wrong when type == 'exp', but 'Rhole'
+  was in the list of parameters, hole was not passed; meant that
+  McMillan17 was wrong in C in v1.5).
+
 - Add time dependence to all relevant Potential functions and method
   (#404).
 

doc/source/potential.rst

@@ -970,10 +970,6 @@ implementation of the classic Chandrasekhar dynamical-friction
 formula, with recent tweaks to better represent the results from
 *N*-body simulations.
 
-Note that there is currently no support for implementing dissipative
-forces in C. Thus, only Python-based integration methods are available
-for any dissipative forces.
-
 .. WARNING::
    Dissipative forces can currently only be used for 3D orbits in ``galpy``. The code should throw an error when they are used for 2D orbits.
 
@@ -992,7 +988,7 @@ used. Adding a new class of potentials to galpy consists of the
 following series of steps (for steps to add a new wrapper potential,
 also see :ref:`the next section <addwrappot>`):
 
-1. 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=`` and two units parameters ``ro=`` and ``vo=`` should be taken as an argument for this class and before performing any other setup, the   ``galpy.potential.Potential.__init__(self,amp=amp,ro=ro,vo=vo,amp_units=)`` method should   be called to setup the amplitude and the system of units; the ``amp_units=`` keyword specifies the physical units of the amplitude parameter (e.g., ``amp_units='velocity2'`` when the units of the amplitude are velocity-squared) 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. 
+1. Implement the new potential in a class that inherits from ``galpy.potential.Potential`` (velocity-dependent forces should inherit from ``galpy.potential.DissipativeForce`` instead; see below for a brief discussion on differences in implementing such forces). The new class should have an ``__init__`` method that sets up the necessary parameters for the class. An amplitude parameter ``amp=`` and two units parameters ``ro=`` and ``vo=`` should be taken as an argument for this class and before performing any other setup, the   ``galpy.potential.Potential.__init__(self,amp=amp,ro=ro,vo=vo,amp_units=)`` method should   be called to setup the amplitude and the system of units; the ``amp_units=`` keyword specifies the physical units of the amplitude parameter (e.g., ``amp_units='velocity2'`` when the units of the amplitude are velocity-squared) 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. 
 
 .. _addpypot:
 
@@ -1048,8 +1044,8 @@ also see :ref:`the next section <addwrappot>`):
     pattern speed (used to compute the Jacobi integral for orbits).
 
   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.
+  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
@@ -1069,7 +1065,7 @@ also see :ref:`the next section <addwrappot>`):
   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.
+  potential has to be passed from python to C (see below).
 
   The ``__init__`` method should be written in such a way that a
   relevant object can be initialized using ``Classname()`` (i.e.,
@@ -1086,7 +1082,10 @@ also see :ref:`the next section <addwrappot>`):
   ``hasC=True`` for potentials for which the forces and potential are
   implemented in C (see below); ``self.hasC_dxdv=True`` for potentials
   for which the (planar) second derivatives are implemented in C;
-  ``self.isNonAxi=True`` for non-axisymmetric potentials.
+  ``self.hasC_dens=True`` for potentials for which the density is
+  implemented in C as well (necessary for them to work with dynamical
+  friction in C); ``self.isNonAxi=True`` for non-axisymmetric
+  potentials.
 
 2. To add a C implementation of the potential, implement it in a .c file under ``potential/potential_c_ext``. Look at ``potential/potential_c_ext/LogarithmicHaloPotential.c`` for the right format for 3D, axisymmetric potentials, or at ``potential/potential_c_ext/LopsidedDiskPotential.c`` for 2D, non-axisymmetric potentials. 
 
@@ -1098,7 +1097,10 @@ also see :ref:`the next section <addwrappot>`):
  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/potential_c_ext/galpy_potentials.h``.
+ is most important (i.e., for those algorithms that evaluate the potential). If you want your potential to be able to be used as the density for the :ref:`ChandrasekharDynamicalFrictionForce <dynamfric_potential>` implementation in C, you need to implement the density in C as well
+
+ * double LogarithmicHaloPotentialDens(double R,double Z, double phi,double t,struct potentialArg * potentialArgs)
+ 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/potential_c_ext/galpy_potentials.h``.
 
 3. Add the potential's function declarations to
 ``potential/potential_c_ext/galpy_potentials.h``
@@ -1120,11 +1122,17 @@ new potential (in the **_parse_pot** function).
 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``.
+``self.hasC_dxdv=True``. If you have implemented the density in C, set
+``self.hasC_dens=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.
 
+Velocity-dependent forces (e.g.,
+:ref:`ChandrasekharDynamicalFrictionForce <dynamfric_potential>`) should inherit from ``galpy.potential.DissipativeForce`` instead of from ``galpy.potential.Potential``. Because such forces are not conservative, you only need to implement the forces themselves, in the same way as for a regular ``Potential``. For dissipative forces, the force-evaluation functions (``Rforce``, etc.) need to take the velocity in cylindrical coordinates as a keyword argument: ``v=[vR,vT,vZ]``. Implementing dissipative forces in C is similar: you only need to implement the forces themselves and the forces should take the velocity in cylindrical coordinates as an additional input, e.g.,
+
+* double ChandrasekharDynamicalFrictionForceRforce(double R,double z, double phi,double t,struct potentialArg * potentialArgs,double vR,double vT,double vz)
+
 .. _addwrappot:
 
 Adding wrapper potentials to the galpy framework

doc/source/reference/potential.rst

@@ -342,7 +342,6 @@ As an example, we integrate the Sun's orbit for 10 Gyr in
 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
@@ -354,7 +353,6 @@ agree better
 >>> 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
 
@@ -365,8 +363,6 @@ We can also compare the rotation curves of these different models
 >>> legend()
 
 .. image:: ../images/mwpotentials-vcirc.png
-   :scale: 40 %
-
 
 
 

galpy/df/jeans.py

@@ -41,6 +41,7 @@ def sigmar(Pot,r,dens=None,beta=0.):
     Pot= flatten_pot(Pot)
     if dens is None:
         dens= lambda r: evaluateDensities(Pot,r*_INVSQRTTWO,r*_INVSQRTTWO,
+                                          phi=numpy.pi/4.,
                                           use_physical=False)
     if callable(beta):
         intFactor= lambda x: numpy.exp(2.*integrate.quad(lambda y: beta(y)/y,
@@ -50,7 +51,8 @@ def sigmar(Pot,r,dens=None,beta=0.):
     return numpy.sqrt(integrate.quad(lambda x: -intFactor(x)*dens(x)
                                      *evaluaterforces(Pot,
                                                       x*_INVSQRTTWO,
-                                                      x*_INVSQRTTWO,
+                                                      x*_INVSQRTTWO, 
+                                                      phi=numpy.pi/4.,
                                                       use_physical=False),
                                      r,numpy.inf)[0]/
                       dens(r)/intFactor(r))

galpy/orbit/Orbits.py

@@ -601,7 +601,7 @@ def from_name(cls,*args,**kwargs):
         _update_keys_named_objects()
         # Stack coordinate-transform parameters, so they can be changed...
         obs= numpy.array([kwargs.get('ro',None),
-                          kwargs.get('vro',None),
+                          kwargs.get('vo',None),
                           kwargs.get('zo',None),
                           kwargs.get('solarmotion',None)],
                          dtype='object')

galpy/orbit/integrateFullOrbit.py

@@ -172,8 +172,7 @@ def _parse_pot(pot,potforactions=False,potfortorus=False):
                 npot+= 1
                 pot_type.append(26)
                 stype= Sigma.get('type','exp')
-                if stype == 'exp' \
-                        or (stype == 'exp' and 'Rhole' in Sigma):
+                if stype == 'exp' and not 'Rhole' in Sigma:
                     pot_args.extend([3,0,
                                      4.*numpy.pi*Sigma.get('amp',1.)*p._amp,
                                      Sigma.get('h',1./3.)])
@@ -259,6 +258,24 @@ def _parse_pot(pot,potforactions=False,potfortorus=False):
             pot_args.extend(p._orb.z(p._orb.t,use_physical=False))
             pot_args.extend([p._amp])
             pot_args.extend([p._orb.t[0],p._orb.t[-1]]) #t_0, t_f
+        elif isinstance(p,potential.ChandrasekharDynamicalFrictionForce):
+            pot_type.append(-7)
+            wrap_npot, wrap_pot_type, wrap_pot_args= \
+                _parse_pot(p._dens_pot,
+                           potforactions=potforactions,potfortorus=potfortorus)
+            pot_args.append(wrap_npot)
+            pot_type.extend(wrap_pot_type)
+            pot_args.extend(wrap_pot_args)
+            pot_args.extend([len(p._sigmar_rs_4interp)])
+            pot_args.extend(p._sigmar_rs_4interp)
+            pot_args.extend(p._sigmars_4interp)
+            pot_args.extend([p._amp])
+            pot_args.extend([-1.,0.,0.,0.,0.,0.,0.,0.]) # for caching
+            pot_args.extend([p._ms,p._rhm,p._gamma**2.,
+                             -1 if not p._lnLambda else p._lnLambda,
+                             p._minr**2.])
+            pot_args.extend([p._sigmar_rs_4interp[0],
+                             p._sigmar_rs_4interp[-1]]) #r_0, r_f
     pot_type= numpy.array(pot_type,dtype=numpy.int32,order='C')
     pot_args= numpy.array(pot_args,dtype=numpy.float64,order='C')
     return (npot,pot_type,pot_args)