From eafbd898daccdf9076136ab3bf615f2ac57ad600 Mon Sep 17 00:00:00 2001
From: Jo Bovy <bovy@astro.utoronto.ca>
Date: Thu, 21 Dec 2023 14:05:06 -0500
Subject: [PATCH] Vendor scipy.integrate.romberg due to impending deprecation
 within scipy

---
 galpy/df/diskdf.py     |   4 +-
 galpy/util/quadpack.py | 119 +++++++++++++++++++++++++++++++++++++++++
 2 files changed, 121 insertions(+), 2 deletions(-)

diff --git a/galpy/df/diskdf.py b/galpy/df/diskdf.py
index ca981f601..ae2ecbb16 100644
--- a/galpy/df/diskdf.py
+++ b/galpy/df/diskdf.py
@@ -32,7 +32,7 @@
 from ..actionAngle import actionAngleAdiabatic
 from ..orbit import Orbit
 from ..potential import PowerSphericalPotential
-from ..util import conversion, save_pickles
+from ..util import conversion, quadpack, save_pickles
 from ..util.ars import ars
 from ..util.conversion import (
     _APY_LOADED,
@@ -2702,7 +2702,7 @@ def _oned_intFunc(x, twodfunc, gfun, hfun, tol, args):
     """Internal function for bovy_dblquad"""
     thisargs = copy.deepcopy(args)
     thisargs.insert(0, x)
-    return integrate.romberg(twodfunc, gfun(x), hfun(x), args=thisargs, tol=tol)
+    return quadpack.romberg(twodfunc, gfun(x), hfun(x), args=thisargs, tol=tol)
 
 
 def bovy_dblquad(func, a, b, gfun, hfun, args=(), tol=1.48e-08):
diff --git a/galpy/util/quadpack.py b/galpy/util/quadpack.py
index 5864cdf78..8ab17c96a 100644
--- a/galpy/util/quadpack.py
+++ b/galpy/util/quadpack.py
@@ -138,3 +138,122 @@ def quadrature(
             AccuracyWarning,
         )
     return val, err
+
+
+# scipy Romberg quadrature that was deprecated in 1.12.0 and associated
+# functions; same BSD-3 license as galpy
+def _difftrap(function, interval, numtraps):
+    """
+    Perform part of the trapezoidal rule to integrate a function.
+    Assume that we had called difftrap with all lower powers-of-2
+    starting with 1. Calling difftrap only returns the summation
+    of the new ordinates. It does _not_ multiply by the width
+    of the trapezoids. This must be performed by the caller.
+        'function' is the function to evaluate (must accept vector arguments).
+        'interval' is a sequence with lower and upper limits
+                   of integration.
+        'numtraps' is the number of trapezoids to use (must be a
+                   power-of-2).
+    """
+    if numtraps == 1:
+        return 0.5 * (function(interval[0]) + function(interval[1]))
+    else:
+        numtosum = numtraps / 2
+        h = float(interval[1] - interval[0]) / numtosum
+        lox = interval[0] + 0.5 * h
+        points = lox + h * numpy.arange(numtosum)
+        s = numpy.sum(function(points), axis=0)
+        return s
+
+
+def _romberg_diff(b, c, k):
+    """
+    Compute the differences for the Romberg quadrature corrections.
+    See Forman Acton's "Real Computing Made Real," p 143.
+    """
+    tmp = 4.0**k
+    return (tmp * c - b) / (tmp - 1.0)
+
+
+def romberg(
+    function,
+    a,
+    b,
+    args=(),
+    tol=1.48e-8,
+    rtol=1.48e-8,
+    divmax=10,
+    vec_func=False,
+):
+    """
+    Romberg integration of a callable function or method.
+
+    Returns the integral of `function` (a function of one variable)
+    over the interval (`a`, `b`).
+
+    If `show` is 1, the triangular array of the intermediate results
+    will be printed. If `vec_func` is True (default is False), then
+    `function` is assumed to support vector arguments.
+
+    Parameters
+    ----------
+    function : callable
+        Function to be integrated.
+    a : float
+        Lower limit of integration.
+    b : float
+        Upper limit of integration.
+
+    Returns
+    -------
+    results : float
+        Result of the integration.
+
+    Other Parameters
+    ----------------
+    args : tuple, optional
+        Extra arguments to pass to function. Each element of `args` will
+        be passed as a single argument to `func`. Default is to pass no
+        extra arguments.
+    tol, rtol : float, optional
+        The desired absolute and relative tolerances. Defaults are 1.48e-8.
+    divmax : int, optional
+        Maximum order of extrapolation. Default is 10.
+    vec_func : bool, optional
+        Whether `func` handles arrays as arguments (i.e., whether it is a
+        "vector" function). Default is False.
+
+    References
+    ----------
+    .. [1] 'Romberg's method' https://en.wikipedia.org/wiki/Romberg%27s_method
+
+    """
+    if numpy.isinf(a) or numpy.isinf(b):  # pragma: no cover
+        raise ValueError("Romberg integration only available " "for finite limits.")
+    vfunc = vectorize1(function, args, vec_func=vec_func)
+    n = 1
+    interval = [a, b]
+    intrange = b - a
+    ordsum = _difftrap(vfunc, interval, n)
+    result = intrange * ordsum
+    resmat = [[result]]
+    err = numpy.inf
+    last_row = resmat[0]
+    for i in range(1, divmax + 1):
+        n *= 2
+        ordsum += _difftrap(vfunc, interval, n)
+        row = [intrange * ordsum / n]
+        for k in range(i):
+            row.append(_romberg_diff(last_row[k], row[k], k + 1))
+        result = row[i]
+        lastresult = last_row[i - 1]
+        err = abs(result - lastresult)
+        if err < tol or err < rtol * abs(result):
+            break
+        last_row = row
+    else:  # pragma: no cover
+        warnings.warn(
+            "divmax (%d) exceeded. Latest difference = %e" % (divmax, err),
+            AccuracyWarning,
+        )
+    return result