Merge branch 'main' into python-3.13

Author: jobovy

.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.19.0
+  rev: v3.19.1
   hooks:
   - id: pyupgrade
     args: [--py38-plus]
@@ -21,7 +21,7 @@ repos:
   hooks:
   - id: rst-backticks
 - repo: https://github.com/hadialqattan/pycln
-  rev: v2.4.0
+  rev: v2.5.0
   hooks:
     - id: pycln
 - repo: https://github.com/PyCQA/isort
@@ -36,6 +36,6 @@ repos:
       args: ["-L", "thisE,thise,mye,tE,te,hist,ro,sav,ccompiler,aas,floatIn,dOmin",
       "-x","doc/source/_static/try-galpy.js"]
 - repo: https://github.com/astral-sh/ruff-pre-commit
-  rev: v0.8.1
+  rev: v0.8.6
   hooks:
     - id: ruff-format

HISTORY.txt

@@ -1,6 +1,8 @@
 v1.10.2 (expected around 2025-03-01)
 ====================
 
+ - Implemented the IAS15 integrator of Rein & Spiegel 2015.
+
  - Combine the drift calculations in the Python leapfrog integrator (#365).
 
  - Move the checks for non-axisymmetric and dissipative potentials from internal to

README.md

@@ -6,7 +6,7 @@
 [galpy](http://www.galpy.org) is a Python package for galactic dynamics. It supports orbit integration in a variety of potentials, evaluating and sampling various distribution functions, and the calculation of action-angle coordinates for all static potentials. `galpy` is an [astropy](http://www.astropy.org/) [affiliated package](http://www.astropy.org/affiliated/) and provides full support for astropy’s [Quantity](http://docs.astropy.org/en/stable/api/astropy.units.Quantity.html) framework for variables with units.
 
 [![image](https://github.com/jobovy/galpy/actions/workflows/build.yml/badge.svg?branch=main)](https://github.com/jobovy/galpy/actions/workflows/build.yml) [![image](https://github.com/jobovy/galpy/actions/workflows/build_windows.yml/badge.svg?branch=main)](https://github.com/jobovy/galpy/actions/workflows/build_windows.yml) [![pre-commit.ci status](https://results.pre-commit.ci/badge/github/jobovy/galpy/main.svg)](https://results.pre-commit.ci/latest/github/jobovy/galpy/main)
-[![image](http://codecov.io/gh/jobovy/galpy/coverage.svg?branch=main)](http://app.codecov.io/gh/jobovy/galpy?branch=main) [![image](https://readthedocs.org/projects/galpy/badge/?version=latest)](http://docs.galpy.org/en/latest/)
+[![image](https://codecov.io/gh/jobovy/galpy/branch/main/graph/badge.svg?token=nIqjwFncfP)](https://codecov.io/gh/jobovy/galpy) [![image](https://readthedocs.org/projects/galpy/badge/?version=latest)](http://docs.galpy.org/en/latest/)
 [![image](http://img.shields.io/pypi/v/galpy.svg)](https://pypi.python.org/pypi/galpy/) [![image](https://img.shields.io/pypi/pyversions/galpy?logo=python&logoColor=white)](https://pypi.python.org/pypi/galpy/) [![image](https://anaconda.org/conda-forge/galpy/badges/version.svg)](https://anaconda.org/conda-forge/galpy) [![image](https://img.shields.io/github/commits-since/jobovy/galpy/latest)](https://github.com/jobovy/galpy/commits/main)
 [![image](http://img.shields.io/badge/license-New%20BSD-brightgreen.svg)](https://github.com/jobovy/galpy/blob/main/LICENSE) [![image](http://img.shields.io/badge/DOI-10.1088/0067%2D%2D0049/216/2/29-blue.svg)](http://dx.doi.org/10.1088/0067-0049/216/2/29) [![image](http://img.shields.io/badge/powered%20by-AstroPy-orange.svg?style=flat)](http://www.astropy.org/) [![image](https://img.shields.io/badge/join-slack-E01563.svg?style=flat&logo=slack&logoWidth=10)](https://join.slack.com/t/galpy/shared_invite/zt-p6upr4si-mX7u8MRdtm~3bW7o8NA_Ww)
 

doc/source/orbit.rst

@@ -1014,9 +1014,11 @@ the ``orbit.integrate`` method. Currently available integrators are
 * rk6_c
 * dopr54_c
 * dop853_c
+* ias15_c
 
-which are Runge-Kutta and Dormand-Prince methods. There are also a
-number of symplectic integrators available
+which are Runge-Kutta, Dormand-Prince methods and the IAS15
+integrator in `Rein & Spiegel (2014) <https://doi.org/10.1093/mnras/stu2164>`_.
+There are also a number of symplectic integrators available
 
 * leapfrog_c
 * symplec4_c
@@ -1045,6 +1047,9 @@ which are speedy and reliable. For example, compare
 >>> timeit(o.integrate(ts,mp,method='dop853_c'))
 # 4.65 ms ± 86.8 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
 
+The ``ias15_c`` method uses adaptive timestepping under the hood, and
+can be used in cases where very high precision is required during orbit integration.
+
 If the C extensions are unavailable for some reason, I recommend using
 the ``odeint`` pure-Python integrator, as it is the fastest. Using the
 same example as above

galpy/orbit/Orbits.py

@@ -1287,13 +1287,15 @@ def check_integrator(method, no_symplec=False):
             "rk6_c",
             "dopr54_c",
             "dop853_c",
+            "ias15_c",
         ]
         if no_symplec:
             symplec_methods = [
                 "leapfrog",
                 "leapfrog_c",
                 "symplec4_c",
                 "symplec6_c",
+                "ias15_c",  # practically speaking, ias15 has the same limitations as symplectic integrators in galpy
             ]
             [valid_methods.remove(symplec_method) for symplec_method in symplec_methods]
         if method.lower() not in valid_methods:
@@ -1385,6 +1387,7 @@ def integrate(
           -  'dopr54_c' for a 5-4 Dormand-Prince integrator in C
           -  'dop853' for a 8-5-3 Dormand-Prince integrator in Python
           -  'dop853_c' for a 8-5-3 Dormand-Prince integrator in C
+          -  'ias15_c' for an adaptive 15th order integrator using Gauß-Radau quadrature (see IAS15 paper) in C
 
         - 2018-10-13 - Written as parallel_map applied to regular Orbit integration - Mathew Bub (UofT)
         - 2018-12-26 - Written to use OpenMP C implementation - Bovy (UofT)
@@ -1743,7 +1746,7 @@ def integrate_dxdv(
 
         Notes
         -----
-        - Possible integration methods are:
+        - Possible integration methods are the non-symplectic ones in galpy:
 
           - 'odeint' for scipy's odeint
           -  'rk4_c' for a 4th-order Runge-Kutta integrator in C
@@ -1761,22 +1764,7 @@ def integrate_dxdv(
             raise AttributeError(
                 "integrate_dxdv is only implemented for 4D (planar) orbits"
             )
-        if method.lower() not in [
-            "odeint",
-            "dop853",
-            "rk4_c",
-            "rk6_c",
-            "dopr54_c",
-            "dop853_c",
-        ]:
-            if "leapfrog" in method.lower() or "symplec" in method.lower():
-                raise ValueError(
-                    f"{method:s} is not a valid `method for integrate_dxdv, because symplectic integrators cannot be used`"
-                )
-            else:
-                raise ValueError(
-                    f"{method:s} is not a valid `method for integrate_dxdv`"
-                )
+        self.check_integrator(method, no_symplec=True)
         pot = flatten_potential(pot)
         _check_potential_dim(self, pot)
         _check_consistent_units(self, pot)
@@ -5168,6 +5156,7 @@ def bruteSOS(
           -  'dopr54_c' for a 5-4 Dormand-Prince integrator in C
           -  'dop853' for a 8-5-3 Dormand-Prince integrator in Python
           -  'dop853_c' for a 8-5-3 Dormand-Prince integrator in C
+          -  'ias15_c' for an adaptive 15th order integrator using Gauß-Radau quadrature (see IAS15 paper) in C
 
         - 2023-05-31 - Written - Bovy (UofT)
 
@@ -5910,6 +5899,7 @@ def plotBruteSOS(
           -  'dopr54_c' for a 5-4 Dormand-Prince integrator in C
           -  'dop853' for a 8-5-3 Dormand-Prince integrator in Python
           -  'dop853_c' for a 8-5-3 Dormand-Prince integrator in C
+          -  'ias15_c' for an adaptive 15th order integrator using Gauß-Radau quadrature (see IAS15 paper) in C
 
         - 2023-05-31 - Written - Bovy (UofT)
 

galpy/orbit/integratePlanarOrbit.py

@@ -618,6 +618,8 @@ def _parse_integrator(int_method):
         int_method_c = 5
     elif int_method.lower() == "dop853_c":
         int_method_c = 6
+    elif int_method.lower() == "ias15_c":
+        int_method_c = 7
     else:
         int_method_c = 0
     return int_method_c

galpy/orbit/orbit_c_ext/integrateFullOrbit.c

@@ -14,6 +14,7 @@
 #include <bovy_symplecticode.h>
 #include <leung_dop853.h>
 #include <bovy_rk.h>
+#include <wez_ias15.h>
 #include <integrateFullOrbit.h>
 //Potentials
 #include <galpy_potentials.h>
@@ -717,6 +718,11 @@ EXPORT void integrateFullOrbit(int nobj,
     odeint_deriv_func= &evalRectDeriv;
     dim= 6;
     break;
+  case 7: //ias15
+    odeint_func= &wez_ias15;
+    odeint_deriv_func= &evalRectForce;
+    dim= 3;
+    break;
   }
 #pragma omp parallel for schedule(dynamic,ORBITS_CHUNKSIZE) private(ii,jj) num_threads(max_threads)
   for (ii=0; ii < nobj; ii++) {
@@ -883,6 +889,11 @@ void integrateOrbit_dxdv(double *yo,
     odeint_deriv_func= &evalRectDeriv_dxdv;
     dim= 12;
     break;
+  case 7: //ias15
+    odeint_func= &wez_ias15;
+    odeint_deriv_func= &evalRectForce;
+    dim= 6;
+    break;
   }
   odeint_func(odeint_deriv_func,dim,yo,nt,-9999.99,t,npot,potentialArgs,
 	      rtol,atol,result,err);

galpy/orbit/orbit_c_ext/integrateLinearOrbit.c

@@ -7,6 +7,7 @@
 #include <math.h>
 #include <bovy_symplecticode.h>
 #include <bovy_rk.h>
+#include <wez_ias15.h>
 #include <leung_dop853.h>
 #include <integrateFullOrbit.h>
 //Potentials
@@ -192,6 +193,11 @@ EXPORT void integrateLinearOrbit(int nobj,
     odeint_deriv_func= &evalLinearDeriv;
     dim= 2;
     break;
+  case 7: //ias15
+    odeint_func= &wez_ias15;
+    odeint_deriv_func= &evalLinearForce;
+    dim= 1;
+    break;
   }
 #pragma omp parallel for schedule(dynamic,ORBITS_CHUNKSIZE) private(ii) num_threads(max_threads)
   for (ii=0; ii < nobj; ii++) {

galpy/orbit/orbit_c_ext/integratePlanarOrbit.c

@@ -8,6 +8,7 @@
 #include <bovy_coords.h>
 #include <bovy_symplecticode.h>
 #include <bovy_rk.h>
+#include <wez_ias15.h>
 #include <leung_dop853.h>
 #include <integrateFullOrbit.h>
 //Potentials
@@ -674,6 +675,11 @@ EXPORT void integratePlanarOrbit(int nobj,
     odeint_deriv_func= &evalPlanarRectDeriv;
     dim= 4;
     break;
+  case 7: //ias15
+    odeint_func= &wez_ias15;
+    odeint_deriv_func= &evalPlanarRectForce;
+    dim= 2;
+    break;
   }
 #pragma omp parallel for schedule(dynamic,ORBITS_CHUNKSIZE) private(ii,jj) num_threads(max_threads)
   for (ii=0; ii < nobj; ii++) {

galpy/potential/SCFPotential.py

@@ -1,9 +1,17 @@
 import hashlib
 
 import numpy
+import scipy
 from numpy.polynomial.legendre import leggauss
+from packaging.version import parse as parse_version
 from scipy import integrate
-from scipy.special import gamma, gammaln, lpmn
+from scipy.special import gamma, gammaln
+
+_SCIPY_VERSION = parse_version(scipy.__version__)
+if _SCIPY_VERSION < parse_version("1.15"):  # pragma: no cover
+    from scipy.special import lpmn
+else:
+    from scipy.special import assoc_legendre_p_all
 
 from ..util import conversion, coords
 from ..util._optional_deps import _APY_LOADED
@@ -412,7 +420,17 @@ def _compute(self, funcTilde, R, z, phi):
         N, L, M = Acos.shape
         r, theta, phi = coords.cyl_to_spher(R, z, phi)
 
-        PP = lpmn(M - 1, L - 1, numpy.cos(theta))[0].T  ##Get the Legendre polynomials
+        ## Get the Legendre polynomials
+        if _SCIPY_VERSION < parse_version("1.15"):  # pragma: no cover
+            PP = lpmn(M - 1, L - 1, numpy.cos(theta))[0].T
+        else:
+            PP = numpy.swapaxes(
+                assoc_legendre_p_all(L - 1, M - 1, numpy.cos(theta), branch_cut=2)[
+                    0, :, :M
+                ],
+                0,
+                1,
+            ).T
         func_tilde = funcTilde(r, N, L)  ## Tilde of the function of interest
 
         func = numpy.zeros(
@@ -515,9 +533,15 @@ def _computeforce(self, R, z, phi=0, t=0):
             dPhi_dphi = self._cached_dPhi_dphi
 
         else:
-            PP, dPP = lpmn(
-                M - 1, L - 1, numpy.cos(theta)
-            )  ##Get the Legendre polynomials
+            ## Get the Legendre polynomials
+            if _SCIPY_VERSION < parse_version("1.15"):  # pragma: no cover
+                PP, dPP = lpmn(M - 1, L - 1, numpy.cos(theta))
+            else:
+                PP, dPP = assoc_legendre_p_all(
+                    L - 1, M - 1, numpy.cos(theta), branch_cut=2, diff_n=1
+                )
+                PP = numpy.swapaxes(PP[:, :M], 0, 1)
+                dPP = numpy.swapaxes(dPP[:, :M], 0, 1)
             PP = PP.T[None, :, :]
             dPP = dPP.T[None, :, :]
             phi_tilde = self._phiTilde(r, N, L)[:, :, numpy.newaxis]
@@ -924,14 +948,13 @@ def integrand(xi, costheta):
         r = _xiToR(xi, a)
         R = r * numpy.sqrt(1 - costheta**2.0)
         z = r * costheta
-        Legandre = lpmn(0, L - 1, costheta)[0].T[numpy.newaxis, :, 0]
+        if _SCIPY_VERSION < parse_version("1.15"):  # pragma: no cover
+            PP = lpmn(0, L - 1, costheta)[0].T[numpy.newaxis, :, 0]
+        else:
+            PP = assoc_legendre_p_all(L - 1, 0, costheta, branch_cut=2)[0].T
         dV = (1.0 + xi) ** 2.0 * numpy.power(1.0 - xi, -4.0)
         phi_nl = (
-            a**3
-            * (1.0 + xi) ** l
-            * (1.0 - xi) ** (l + 1.0)
-            * _C(xi, N, L)[:, :]
-            * Legandre
+            a**3 * (1.0 + xi) ** l * (1.0 - xi) ** (l + 1.0) * _C(xi, N, L)[:, :] * PP
         )
         param[0] = R
         param[1] = z
@@ -1089,15 +1112,22 @@ def integrand(xi, costheta, phi):
         r = _xiToR(xi, a)
         R = r * numpy.sqrt(1 - costheta**2.0)
         z = r * costheta
-        Legandre = lpmn(L - 1, L - 1, costheta)[0].T[numpy.newaxis, :, :]
+        if _SCIPY_VERSION < parse_version("1.15"):  # pragma: no cover
+            PP = lpmn(L - 1, L - 1, costheta)[0].T[numpy.newaxis, :, :]
+        else:
+            PP = numpy.swapaxes(
+                assoc_legendre_p_all(L - 1, L - 1, costheta, branch_cut=2)[0][:, :L],
+                0,
+                1,
+            ).T[numpy.newaxis, :, :]
         dV = (1.0 + xi) ** 2.0 * numpy.power(1.0 - xi, -4.0)
 
         phi_nl = (
             -(a**3)
             * (1.0 + xi) ** l
             * (1.0 - xi) ** (l + 1.0)
             * _C(xi, N, L)[:, :, numpy.newaxis]
-            * Legandre
+            * PP
         )
 
         return (

galpy/potential/potential_c_ext/PowerSphericalPotentialwCutoff.c

@@ -7,7 +7,7 @@
 //PowerSphericalPotentialwCutoff
 //3  arguments: amp, alpha, rc
 double mass(double r2,double alpha, double rc){
-  return 2. * M_PI * pow ( rc , 3. - alpha ) * ( gsl_sf_gamma ( 1.5 - 0.5 * alpha ) - gsl_sf_gamma_inc ( 1.5 - 0.5 * alpha , r2 / rc / rc ) );
+  return 2. * M_PI * pow ( rc , 3. - alpha ) * ( gsl_sf_gamma ( 1.5 - 0.5 * alpha ) * gsl_sf_gamma_inc_P ( 1.5 - 0.5 * alpha , r2 / rc / rc ) );
 }
 double PowerSphericalPotentialwCutoffEval(double R,double Z, double phi,
 					  double t,

galpy/util/quadpack.py

@@ -5,7 +5,6 @@
 
 import numpy
 from scipy.integrate import fixed_quad, quad
-from scipy.integrate._quadrature import vectorize1
 
 
 def _infunc(x, func, gfun, hfun, more_args, epsrel, epsabs):