Merge branch 'main' into numpy-docstrings

.bumpversion.cfg

@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 1.9.0.dev0
+current_version = 1.9.1.dev0
 parse = (?P<major>\d+)\.(?P<minor>\d+)\.(?P<patch>\d+)(\.(?P<release>[a-z]+)(?P<dev>\d+))?
 serialize =
 	{major}.{minor}.{patch}.{release}{dev}
@@ -14,16 +14,16 @@ values =
 [bumpversion:part:dev]
 
 [bumpversion:file:setup.py]
-search = version='{current_version}'
-replace = version='{new_version}'
+search = version="{current_version}"
+replace = version="{new_version}"
 
 [bumpversion:file:galpy/__init__.py]
 search = __version__ = "{current_version}"
 replace = __version__ = "{new_version}"
 
 [bumpversion:file:doc/source/conf.py]
-search = version = '{current_version}'
-replace = version = '{new_version}'
+search = version = "{current_version}"
+replace = version = "{new_version}"
 
 [bumpversion:file:CITATION.cff]
 search = version: {current_version}

.github/workflows/wheels.yml

@@ -88,7 +88,7 @@ jobs:
           rsync -azv gsl/lib/ /usr/local/lib/
           echo "REPAIR_LIBRARY_PATH=/usr/local/lib/" >> $GITHUB_ENV # https://github.com/pypa/cibuildwheel/issues/816#issuecomment-916197079
       - name: Build wheels
-        uses: pypa/cibuildwheel@v2.13
+        uses: pypa/cibuildwheel@v2.14
         env:
           CIBW_BUILD: ${{ matrix.python }}-${{ matrix.buildplat[1] }}_${{ matrix.buildplat[2] }}
           CIBW_BUILD_VERBOSITY: 1

.pre-commit-config.yaml

@@ -12,7 +12,7 @@ repos:
   - id: check-executables-have-shebangs
   - id: check-yaml
 - repo: https://github.com/asottile/pyupgrade
-  rev: v3.4.0
+  rev: v3.8.0
   hooks:
   - id: pyupgrade
     args: [--py38-plus]
@@ -30,10 +30,10 @@ repos:
   - id: isort
     name: isort (python)
 - repo: https://github.com/codespell-project/codespell
-  rev: "v2.2.4"
+  rev: "v2.2.5"
   hooks:
     - id: codespell
-      args: ["-L", "thisE,thise,mye,tE,te,hist,ro,sav,ccompiler",
+      args: ["-L", "thisE,thise,mye,tE,te,hist,ro,sav,ccompiler,aas",
       "-x","doc/source/_static/try-galpy.js"]
 - repo: https://github.com/psf/black
   rev: 23.3.0

CITATION.cff

@@ -5,7 +5,7 @@ authors:
   given-names: "Jo"
   orcid: "https://orcid.org/0000-0001-6855-442X"
 title: "galpy"
-version: 1.9.0.dev0
+version: 1.9.1.dev0
 url: "https://github.com/jobovy/galpy"
 preferred-citation:
   type: article

HISTORY.txt

@@ -1,4 +1,11 @@
-v1.9.0 (XXXX-XX-XX)
+v1.9.1 (expected around 2023-11-01)
+===================
+
+- Allow vector inputs of solar parameters to Orbit initialization: ro, zo, vo,
+  and solarmotion (#595). Useful when sampling over the uncertainty in the solar
+  parameters.
+
+v1.9.0 (2023-07-02)
 ===================
 
 - Add specialized integration method to determine surfaces of sections of orbits and

README.md

@@ -48,8 +48,8 @@ PYTHON VERSIONS AND DEPENDENCIES
 
 `galpy` supports Python 3. Specifically, galpy supports Python 3.8, 3.9, 3.10,
 and 3.11. GitHub Actions CI builds regularly check support for
-Python 3.10 (and of 3.8, and 3.9 using a more limited, core set of tests)
-on Linux and Windows (and 3.10 on Mac OS). Python 2.7 is no longer supported.
+Python 3.11 (and of 3.8, 3.9, and 3.11 using a more limited, core set of tests)
+on Linux and Windows (and 3.11 on Mac OS). Python 2.7 is no longer supported.
 
 This package requires [Numpy](https://numpy.org/),
 [Scipy](http://www.scipy.org/), and

doc/source/conf.py

@@ -116,7 +116,7 @@ def linkcode_resolve(domain, info):
 # built documents.
 #
 # The short X.Y version.
-version = "1.9.0.dev0"
+version = "1.9.1.dev0"
 # The full version, including alpha/beta/rc tags.
 release = version
 on_rtd = os.environ.get("READTHEDOCS", None) == "True"

doc/source/installation.rst

@@ -367,8 +367,8 @@ technique or its variants as usual.
 
 .. _install_pyodide:
 
-**NEW IN v1.8** Using ``galpy`` in web applications
-++++++++++++++++++++++++++++++++++++++++++++++++++++
+Using ``galpy`` in web applications
++++++++++++++++++++++++++++++++++++
 
 ``galpy`` can be compiled to `WebAssembly <https://webassembly.org/>`__ using the `emscripten <https://emscripten.org/>`__ compiler. In particular, ``galpy`` is part of the `pyodide <https://pyodide.org/en/stable/>`__ Python distribution for the browser, meaning that ``galpy`` can be used on websites without user installation and it still runs at the speed of a compiled language. This powers, for example, the :ref:`Try galpy <try_galpy>` interactive session on this documentation's home page. Thus, it is easy to, e.g., build web-based, interactive galactic-dynamics examples or tutorials without requiring users to install the scientific Python stack and ``galpy`` itself.
 

doc/source/orbit.rst

@@ -221,6 +221,8 @@ multiple objects. Specifically, the initial conditions can be:
    * None: assumed to be the Sun; if None occurs in a list it is assumed to be the Sun *and all other items in the list are assumed to be [ra,dec,...]*; cannot be combined with Quantity lists
    * lists of scalar phase-space coordinates arranged as in the first bullet above (so things like [R,vR,...] where R,vR are scalars in internal units
 
+The solar parameters ``ro``, ``zo``, ``vo``, and ``solarmotion`` can also be specified as arrays with the same shape as the resulting ``Orbit`` instance. This can be useful, for example, when sampling over the uncertainty in the Sun's position and velocity in the Milky Way.
+
 .. TIP::
    For multiple object initialization using an array or SkyCoord, arbitrary input shapes are supported.
 
@@ -253,6 +255,18 @@ and with a SkyCoord:
 As before, you can use the ``SkyCoord`` Galactocentric frame
 specification here.
 
+An example initialization with array input for ``ro`` is:
+
+>>> vxvvs= numpy.array([[1.,0.1,1.,0.1,-0.2,1.5] for ii in range(2)])
+>>> orbits= Orbit(vxvvs,ro=numpy.array([8.,8.5]))
+>>> print(orbits.R())
+# [ 8.   8.5]
+
+Array inputs can also be used for ``zo``, ``vo``, and ``solarmotion``. When using arrays,
+their shapes need to agree with the shape of the resulting ``Orbit`` instance (for
+``solarmotion``, the shape needs to be ``[3,*shape_orbit]``).
+
+
 ``Orbit`` instances containing multiple objects act like numpy arrays
 in many ways, but have some subtly different behaviors for some
 functions. For example, one can do:
@@ -295,8 +309,8 @@ an array with shape ``(6,5,ntimes)`` here.
 
 .. _orbfromname:
 
-**UPDATED IN v1.8** Initialization from an object's name
-*********************************************************
+Initialization from an object's name
+************************************
 
 A convenience method, ``Orbit.from_name``, is also available to initialize
 orbits from the name of an object. For example, for the star `Lacaille 8760 <https://en.wikipedia.org/wiki/Lacaille_8760>`__:
@@ -509,8 +523,8 @@ this can be overwritten). A simple example is
 
 .. _orbintegration-noninertial:
 
-**NEW in v1.8** Orbit integration in non-inertial frames
----------------------------------------------------------
+Orbit integration in non-inertial frames
+----------------------------------------
 
 The default assumption in ``galpy`` is that the frame that an orbit is
 integrated in is an inertial one. However, ``galpy`` also supports
@@ -1043,6 +1057,8 @@ same example as above
 >>> timeit(o.integrate(ts,mp,method='dop853'))
 # 1.61 s ± 218 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
 
+.. _orbitsos:
+
 **NEW in v1.9** Surfaces of section
 -----------------------------------
 
@@ -1410,7 +1426,7 @@ the Milky Way.
 
 .. _orbit-example-barycentric-acceleration-LMC:
 
-**NEW in v1.8** Example: Including the Milky Way center's barycentric acceleration due to the Large Magellanic Cloud in orbit integrations
+Example: Including the Milky Way center's barycentric acceleration due to the Large Magellanic Cloud in orbit integrations
 ---------------------------------------------------------------------------------------------------------------------------
 
 Observations over the last few decades have revealed that the Large Magellanic

doc/source/potential.rst

@@ -430,8 +430,8 @@ you find any problems with this.
 
 .. _scf_potential_docs:
 
-**UPDATED IN v1.8** General density/potential pairs with basis-function expansions
-------------------------------------------------------------------------------------
+General density/potential pairs with basis-function expansions
+--------------------------------------------------------------
 
 ``galpy`` allows for the potential and forces of general,
 time-independent density functions to be computed by expanding the

doc/source/streamdf.rst

@@ -509,7 +509,7 @@ for the difference in :math:`v_Y` as a function of unperturbed :math:`X`:
 
 .. _streamspray-tutorial:
 
-**NEW in v1.8** Particle-spray modeling of streams with ``streamspraydf``
+Particle-spray modeling of streams with ``streamspraydf``
 ---------------------------------------------------------
 
 ``galpy`` also contains an implementation of the particle-spray method

doc/source/whatsnew.rst

@@ -5,6 +5,62 @@ This page gives some of the key improvements in each galpy
 version. See the ``HISTORY.txt`` file in the galpy source for full
 details on what is new and different in each version.
 
+v1.9
++++++
+
+Version 1.9 contains two major new additions since version 1.8, a bunch of small
+updates and bug fixes, and a few deprecations. The major new additions are:
+
+* Specialized support for calculating and displaying surfaces of section for orbits.
+  Surfaces of section for 3D and 2D potentials can now be computed using a dedicated
+  integration technique ``Orbit.SOS`` or using a brute-force technique ``Orbit.bruteSOS``
+  (for orbits for which te dedicated technique fails, e.g., orbits in a rotating frame).
+  There is also support for directly plotting surfaces of section in ``Orbit.plotSOS``
+  and ``Orbit.plotBruteSOS``. See :ref:`orbitsos` for more info.
+
+* There is now general support for action-angle and reverse action-angle transformations
+  for 1D potentials using the ``galpy.actionAngle.actionAngleVertical`` and
+  ``galpy.actionAngle.actionAngleVerticalInverse`` classes. The inverse transformation is
+  computing using a robust variation of the torus-mapper.
+
+Important deprecations are:
+
+* ``galpy.util.bovy_coords``, ``galpy.util.bovy_plot``, and
+  ``galpy.util.bovy_conversion`` have now been removed in favor of their ``bovy_``
+  lacking versions.
+
+* ``phiforce`` and associated Potential methods and functions have been removed in favor
+  of ``phitorque`` and associated methods and functions.
+
+Other user-facing improvements and additions are:
+
+* Potential classes, methods, and functions:
+
+  * Made the potential input explicitly positional-only for all ``galpy.potential``
+    functions to avoid errors when specifying it as a keyword argument.
+
+  * Added general support for DissipativeForce instances in 2D.
+
+  * Implemented NonInertialFrameForce in 2D.
+
+  * Allow potentials' density and surface density to be plotted on physical axes.
+
+* New and improved ``Orbit`` methods:
+
+  * Added ``Orbit.animate3d`` to display a 3D animation of an integrated orbit
+    with an optional Milky-Way representation at the origin when plotting x,y,z.
+
+  * Improved the performance of Orbit.animate performance by using webgl, some UI
+    tweaks. Also fixed using Orbit.animate in jupyterlab and retrolab.
+
+* Improvements to spherical distribution functions:
+
+  * Made it possible to use an interpSphericalPotential as the potential in spherical
+    distribution functions.
+
+  * Added a method, ``dMdE``, to calculate the differential energy distribution of
+    spherical distribution functions.
+
 v1.8
 +++++
 

galpy/__init__.py

@@ -1,4 +1,4 @@
-__version__ = "1.9.0.dev0"
+__version__ = "1.9.1.dev0"
 # Check whether a new version is available
 import datetime
 import http.client

galpy/actionAngle/actionAngleAdiabatic.py

@@ -152,9 +152,9 @@ def _evaluate(self, *args, **kwargs):
                 for ii in range(len(R)):
                     targs = (R[ii], vR[ii], vT[ii], z[ii], vz[ii])
                     tjr, tlz, tjz = self(*targs, **copy.copy(kwargs))
-                    ojr[ii] = tjr
-                    ojz[ii] = tjz
-                    olz[ii] = tlz
+                    ojr[ii] = tjr[0]
+                    ojz[ii] = tjz[0]
+                    olz[ii] = tlz[0]
                 return (ojr, olz, ojz)
             else:
                 if kwargs.get("_justjr", False):
@@ -180,7 +180,11 @@ def _evaluate(self, *args, **kwargs):
                     )
                 else:
                     axiJ = self._aAS(R[0], vR[0], vT[0], 0.0, 0.0, _Jz=Jz)
-                    return (axiJ[0], axiJ[1], Jz)
+                    return (
+                        numpy.atleast_1d(axiJ[0]),
+                        numpy.atleast_1d(axiJ[1]),
+                        numpy.atleast_1d(Jz),
+                    )
 
     def _EccZmaxRperiRap(self, *args, **kwargs):
         """
@@ -254,10 +258,10 @@ def _EccZmaxRperiRap(self, *args, **kwargs):
                     tecc, tzmax, trperi, trap = self._EccZmaxRperiRap(
                         *targs, **copy.copy(kwargs)
                     )
-                    oecc[ii] = tecc
-                    ozmax[ii] = tzmax
-                    orperi[ii] = trperi
-                    orap[ii] = trap
+                    oecc[ii] = tecc[0]
+                    ozmax[ii] = tzmax[0]
+                    orperi[ii] = trperi[0]
+                    orap[ii] = trap[0]
                 return (oecc, ozmax, orperi, orap)
             else:
                 if _dim(self._pot) == 3:

galpy/actionAngle/actionAngleAdiabaticGrid.py

@@ -142,7 +142,7 @@ def __init__(
                             numpy.sqrt(2.0 * y[jj] * self._EzZmaxs[ii]),
                             _justjz=True,
                             **kwargs
-                        )[2]
+                        )[2][0]
                         if jj == nEz - 1:
                             jzEzzmax[ii] = jz[ii, jj]
         for ii in range(nR):
@@ -261,7 +261,7 @@ def __init__(
                                 0.0,
                                 _justjr=True,
                                 **kwargs
-                            )[0]
+                            )[0][0]
                         except UnboundError:  # pragma: no cover
                             raise
                         if jj == 0:

galpy/actionAngle/actionAngleStaeckel.py

@@ -200,9 +200,9 @@ def _evaluate(self, *args, **kwargs):
                     except (TypeError, IndexError):
                         tkwargs["delta"] = delta
                     tjr, tlz, tjz = self(*targs, **tkwargs)
-                    ojr[ii] = tjr
-                    ojz[ii] = tjz
-                    olz[ii] = tlz
+                    ojr[ii] = tjr[0]
+                    ojz[ii] = tjz[0]
+                    olz[ii] = tlz[0]
                 return (ojr, olz, ojz)
             else:
                 # Set up the actionAngleStaeckelSingle object
@@ -481,7 +481,7 @@ def _uminumaxvmin(self, *args, **kwargs):
         HISTORY:
            2017-12-12 - Written - Bovy (UofT)
         """
-        delta = kwargs.pop("delta", self._delta)
+        delta = numpy.atleast_1d(kwargs.pop("delta", self._delta))
         if len(args) == 5:  # R,vR.vT, z, vz
             R, vR, vT, z, vz = args
         elif len(args) == 6:  # R,vR.vT, z, vz, phi
@@ -552,19 +552,16 @@ def _uminumaxvmin(self, *args, **kwargs):
                 for ii in range(len(R)):
                     targs = (R[ii], vR[ii], vT[ii], z[ii], vz[ii])
                     tkwargs = copy.copy(kwargs)
-                    try:
-                        tkwargs["delta"] = delta[ii]
-                    except TypeError:
-                        tkwargs["delta"] = delta
+                    tkwargs["delta"] = delta[ii] if len(delta) > 1 else delta[0]
                     tumin, tumax, tvmin = self._uminumaxvmin(*targs, **tkwargs)
-                    oumin[ii] = tumin
-                    oumax[ii] = tumax
-                    ovmin[ii] = tvmin
+                    oumin[ii] = tumin[0]
+                    oumax[ii] = tumax[0]
+                    ovmin[ii] = tvmin[0]
                 return (oumin, oumax, ovmin)
             else:
                 # Set up the actionAngleStaeckelSingle object
                 aASingle = actionAngleStaeckelSingle(
-                    R[0], vR[0], vT[0], z[0], vz[0], pot=self._pot, delta=delta
+                    R[0], vR[0], vT[0], z[0], vz[0], pot=self._pot, delta=delta[0]
                 )
                 umin, umax = aASingle.calcUminUmax()
                 vmin = aASingle.calcVmin()
@@ -975,8 +972,8 @@ def calcUminUmax(self, **kwargs):
                     try:
                         umin = optimize.brentq(
                             _JRStaeckelIntegrandSquared,
-                            rstart,
-                            self._ux - eps,
+                            numpy.atleast_1d(rstart)[0],
+                            numpy.atleast_1d(self._ux)[0] - eps,
                             (
                                 E,
                                 L,
@@ -1011,8 +1008,8 @@ def calcUminUmax(self, **kwargs):
                 )
                 umax = optimize.brentq(
                     _JRStaeckelIntegrandSquared,
-                    self._ux + eps,
-                    rend,
+                    numpy.atleast_1d(self._ux)[0] + eps,
+                    numpy.atleast_1d(rend)[0],
                     (
                         E,
                         L,