Replace scipy's romberg with quad. (#412)

* Replace scipy's romberg with quad.

Romberg will be removed from scipy in v1.15.0

* Fix keyword argument name.

* Don't need error.

* Update doctests.

* Updat change log.

CHANGES.rst

@@ -12,6 +12,15 @@ sbpy.names
   rates in the LSST survey era [#406]
 
 
+Other Changes and Additions
+---------------------------
+
+sbpy.activity.gas
+^^^^^^^^^^^^^^^^^
+- Replaced calls to the deprecated function `scipy.integrate.romberg` with
+  `scipy.integrate.quad`.  [#412]
+
+
 0.5.0 (2024-08-28)
 ==================
 

docs/sbpy/activity/gas.rst

@@ -148,9 +148,9 @@ number of molecules in an aperture.  Parent and daughter data is provided via
   >>> Q = 1e28 / u.s        # water production rate
   >>> coma = gas.VectorialModel(Q, water, hydroxyl)
   >>> print(coma.column_density(10 * u.km))    # doctest: +FLOAT_CMP
-  2.8976722840952486e+17 1 / m2
+  2.8974972922255264e+17 1 / m2
   >>> print(coma.total_number(1000 * u.km))    # doctest: +FLOAT_CMP
-  6.995158827300034e+29
+  6.995084479934798e+29
 
 Production Rate calculations
 ----------------------------

sbpy/activity/gas/core.py

@@ -24,7 +24,7 @@
 try:
     import scipy
     from scipy import special
-    from scipy.integrate import quad, dblquad, romberg
+    from scipy.integrate import quad, dblquad
     from scipy.interpolate import CubicSpline, PPoly
 except ImportError:
     scipy = None
@@ -1635,8 +1635,8 @@ def _column_density_at_rho(self, rho: np.float64) -> np.float64:
         def column_density_integrand(z):
             return self._volume_density(np.sqrt(z**2 + rhosq))
 
-        c_dens = 2 * romberg(column_density_integrand, 0,
-                             z_max, rtol=0.0001, divmax=50)
+        c_dens = 2 * quad(column_density_integrand, 0,
+                          z_max, epsrel=0.0001)[0]
 
         # result is in 1/m^2
         return c_dens
@@ -1745,14 +1745,13 @@ def _calc_num_fragments_grid(self) -> np.float64:
         def vol_integrand(r, r_func):
             return r_func(r) * r**2
 
-        r_int = romberg(
+        r_int = quad(
             vol_integrand,
             0,
             max_r,
             args=(self.vmr.volume_density_interpolation,),
-            rtol=0.0001,
-            divmax=20,
-        )
+            epsrel=0.0001,
+        )[0]
         return 4 * np.pi * r_int
 
     def _column_density(self, rho) -> np.float64: