diff --git a/ares/__init__.py b/ares/__init__.py
index f52e6504b..43a52e181 100755
--- a/ares/__init__.py
+++ b/ares/__init__.py
@@ -1,13 +1,14 @@
 import os as _os
-import imp as _imp
+import importlib as _imp
+# updates to imp: https://docs.python.org/3/whatsnew/3.12.html#whatsnew312-removed-imp
 
 _HOME = _os.environ.get('HOME')
 
 # Load custom defaults    
-if _os.path.exists('{!s}/.ares/defaults.py'.format(_HOME)):
+if _os.path.exists('~/.ares/defaults.py'.format(_HOME)):
     (_f, _filename, _data) =\
-        _imp.find_module('defaults', ['{!s}/.ares/'.format(_HOME)])
-    rcParams = _imp.load_module('defaults.py', _f, _filename, _data).pf
+        _imp.util.find_spec('defaults', ['~/.ares/'.format(_HOME)])
+    rcParams = _imp.import_module('defaults.py', _f, _filename, _data).pf
 else:
     rcParams = {}
 
diff --git a/ares/analysis/MetaGalacticBackground.py b/ares/analysis/MetaGalacticBackground.py
index 0cf6244f3..0c9525ec8 100755
--- a/ares/analysis/MetaGalacticBackground.py
+++ b/ares/analysis/MetaGalacticBackground.py
@@ -14,7 +14,7 @@
 from ..util import labels
 from ..util.Pickling import read_pickle_file
 import matplotlib.pyplot as pl
-from scipy.integrate import trapz
+from scipy.integrate import   trapezoid
 from ..util.ReadData import flatten_energies
 from ..physics.Constants import erg_per_ev, J21_num, h_P, c, E_LL, E_LyA, \
     sqdeg_per_std
diff --git a/ares/analysis/MultiPhaseMedium.py b/ares/analysis/MultiPhaseMedium.py
index ab4748393..2d3d95482 100755
--- a/ares/analysis/MultiPhaseMedium.py
+++ b/ares/analysis/MultiPhaseMedium.py
@@ -15,9 +15,8 @@
 import matplotlib.pyplot as pl
 from ..util.Stats import get_nu
 from ..util.Pickling import read_pickle_file
-from scipy.misc import derivative
 from ..physics.Constants import *
-from scipy.integrate import cumtrapz
+from scipy.integrate import cumulative_trapezoid
 from scipy.interpolate import interp1d
 from ..physics import Cosmology, Hydrogen
 from ..util.SetDefaultParameterValues import *
@@ -341,7 +340,7 @@ def tau_CMB(self, include_He=True, z_HeII_EoR=3.):
 
         integrand *= sigma_T * dldz
 
-        tau = cumtrapz(integrand, self.history_asc['z'], initial=0)
+        tau = cumulative_trapezoid(integrand, self.history_asc['z'], initial=0)
 
         tau[self.history_asc['z'] > 100] = 0.0
 
@@ -803,7 +802,7 @@ def tau_post_EoR(self, include_He=True, z_HeII_EoR=3.):
 
         integrand *= sigma_T * dldz
 
-        tau = cumtrapz(integrand, ztmp, initial=0)
+        tau = cumulative_trapezoid(integrand, ztmp, initial=0)
 
         return ztmp, tau
 
diff --git a/ares/analysis/TurningPoints.py b/ares/analysis/TurningPoints.py
index 55802d94b..d9f6b17e2 100755
--- a/ares/analysis/TurningPoints.py
+++ b/ares/analysis/TurningPoints.py
@@ -11,8 +11,8 @@
 """
 
 import numpy as np
+import numdifftools as nd
 from ..util import ParameterFile
-from scipy.misc import derivative
 from scipy.optimize import minimize
 from ..physics.Constants import nu_0_mhz
 from ..util.Math import central_difference
@@ -202,9 +202,10 @@ def is_stopping_point(self, z, dTb):
                         continue
 
                 else:
+                    raise NotImplemented('havent revisited since scipy.misc.derivative Deprecation')
                     # Compute curvature at turning point (mK**2 / MHz**2)
                     nuTP = nu_0_mhz / (1. + zTP)
-                    d2 = float(derivative(lambda zz: splev(zz, Bspl_fit1),
+                    d2 = float(nd.Derivative(lambda zz: splev(zz, Bspl_fit1),
                         x0=float(zTP), n=2, dx=1e-4, order=5) * nu_0_mhz**2 / nuTP**4)
 
                     self.turning_points[TP] = (zTP, TTP, d2)
diff --git a/ares/phenom/Tanh21cm.py b/ares/phenom/Tanh21cm.py
index d2ba25aa7..d2c192524 100644
--- a/ares/phenom/Tanh21cm.py
+++ b/ares/phenom/Tanh21cm.py
@@ -12,8 +12,8 @@
 
 import time
 import numpy as np
+import numdifftools as nd
 from ..util import ParameterFile
-from scipy.misc import derivative
 from ..physics import Hydrogen, Cosmology
 from ..physics.Constants import k_B, J21_num, nu_0_mhz
 from ..physics.RateCoefficients import RateCoefficients
@@ -52,7 +52,7 @@ def __init__(self, **kwargs):
         self.hydr = Hydrogen(cosm=self.cosm, **kwargs)
 
     def dTgas_dz(self, z):
-        return derivative(self.cosm.Tgas, x0=z)
+        return nd.Derivative(self.cosm.Tgas)(z)
 
     def electron_density(self, z):
         return np.interp(z, self.cosm.thermal_history['z'],
diff --git a/ares/physics/Cosmology.py b/ares/physics/Cosmology.py
index 94e4c67fb..f6b599125 100755
--- a/ares/physics/Cosmology.py
+++ b/ares/physics/Cosmology.py
@@ -12,7 +12,7 @@
 import os
 import numpy as np
 from ..data import ARES
-from scipy.misc import derivative
+import numdifftools as nd
 from scipy.optimize import fsolve
 from scipy.integrate import quad, ode
 from ..util.Math import interp1d
@@ -462,7 +462,7 @@ def cooling_rate(self, z, T=None):
             ##s
             #func = lambda zz: np.interp(zz, self.inits['z'], self.inits['Tk'])
 
-            dTdz = derivative(self._Tgas_CosmoRec, z, dx=1e-2)
+            dTdz = nd.Derivative(self._Tgas_CosmoRec)(z)
 
             xe = np.interp(z, self.inits['z'], self.inits['xe'])
 
@@ -483,7 +483,7 @@ def cooling_rate(self, z, T=None):
             return dTdz + xe_cool * mult
 
         else:
-            return derivative(self.Tgas, z)
+            return nd.Derivative(self.Tgas)(z)
 
     def log_cooling_rate(self, z):
         if self.pf['approx_thermal_history'] == 'exp':
diff --git a/ares/physics/DustEmission.py b/ares/physics/DustEmission.py
index fdec924b1..f4c3227fd 100644
--- a/ares/physics/DustEmission.py
+++ b/ares/physics/DustEmission.py
@@ -12,7 +12,7 @@
 """
 
 import numpy as np
-from scipy.integrate import simps
+from scipy.integrate import simpson
 from ares.physics.Constants import c, h, k_B, g_per_msun, cm_per_kpc, Lsun
 
 # T_dust parameters
@@ -376,8 +376,8 @@ def __T_dust(self, z, L_nu, tau_nu, R_dust, T_cmb):
             tmp_cmb = 8 * np.pi * h / c**2 * cmb_kappa_nu * (cmb_freqs[None, :, None])**3 \
                 / (np.exp(h * cmb_freqs[None,:,None] / k_B / T_cmb[:,None,:]) - 1)
 
-            tmp_power = simps(tmp_stellar, self.frequencies, axis = 1)
-            tmp_power += simps(tmp_cmb, cmb_freqs, axis = 1)
+            tmp_power = simpson(tmp_stellar, self.frequencies, axis = 1)
+            tmp_power += simpson(tmp_cmb, cmb_freqs, axis = 1)
 
             if self.pf.get('pop_dust_experimental'):
                 print("power =", tmp_power)
@@ -519,7 +519,7 @@ def Luminosity(self, z, wave = 3e5, band=None, idnum=None, window=1,
                 fmax = c / (8 * 1e-4)
                 fmin = c / (1000 * 1e-4)
                 freqs, luminosities = self.dust_sed(fmin, fmax, 1000)
-                luminosities = simps(luminosities[:,:,index], freqs, axis = 1)
+                luminosities = simpson(luminosities[:,:,index], freqs, axis = 1)
 
         # is not cached, we calculate everything for the given z and wave
         else:
@@ -557,7 +557,7 @@ def Luminosity(self, z, wave = 3e5, band=None, idnum=None, window=1,
                 kappa_nu = 0.1 * (nu / 1e12)**2
                 luminosities = 8 * np.pi * h / c**2 * nu[None,:]**3 * kappa_nu[None,:] \
                     / (np.exp(h * nu[None,:] / k_B / T_dust[:,None]) - 1) * (M_dust[:,None] * g_per_msun)
-                luminosities = simps(luminosities, nu, axis = 1)
+                luminosities = simpson(luminosities, nu, axis = 1)
 
 
         if idnum is not None:
diff --git a/ares/physics/ExcursionSet.py b/ares/physics/ExcursionSet.py
index 04d5140d9..049291f2c 100644
--- a/ares/physics/ExcursionSet.py
+++ b/ares/physics/ExcursionSet.py
@@ -15,9 +15,8 @@
 from .Cosmology import Cosmology
 from ..util.Math import central_difference
 from ..util.ParameterFile import ParameterFile
-from scipy.integrate import simps, quad
+from scipy.integrate import simpson, quad
 from scipy.interpolate import interp1d
-from scipy.misc import derivative
 
 two_pi = 2. * np.pi
 four_pi = 4. * np.pi
diff --git a/ares/physics/HaloMassFunction.py b/ares/physics/HaloMassFunction.py
index 333ab0cf6..0a33d5c5b 100755
--- a/ares/physics/HaloMassFunction.py
+++ b/ares/physics/HaloMassFunction.py
@@ -19,10 +19,9 @@
 from ..data import ARES
 from types import FunctionType
 from ..util import ParameterFile
-from scipy.misc import derivative
 from scipy.optimize import fsolve
 from ..util.Warnings import no_hmf
-from scipy.integrate import cumtrapz, simps
+from scipy.integrate import cumulative_trapezoid, simpson
 from ..util.PrintInfo import print_hmf
 from ..util.ProgressBar import ProgressBar
 from ..util.ParameterFile import ParameterFile
@@ -437,14 +436,14 @@ def _get_ngtm_mgtm_from_dndm(self):
                 mf_func = InterpolatedUnivariateSpline(np.log(m), np.log(dndlnm), k=1)
                 mf = mf_func(m_upper)
 
-                int_upper_n = simps(np.exp(mf), dx=m_upper[2] - m_upper[1], even='first')
-                int_upper_m = simps(np.exp(m_upper + mf), dx=m_upper[2] - m_upper[1], even='first')
+                int_upper_n = simpson(np.exp(mf), dx=m_upper[2] - m_upper[1], even='first')
+                int_upper_m = simpson(np.exp(m_upper + mf), dx=m_upper[2] - m_upper[1], even='first')
             else:
                 int_upper_n = 0
                 int_upper_m = 0
 
-            ngtm_ = np.concatenate((cumtrapz(dndlnm[::-1], dx=np.log(m[1]) - np.log(m[0]))[::-1], np.zeros(1)))
-            mgtm_ = np.concatenate((cumtrapz(m[::-1] * dndlnm[::-1], dx=np.log(m[1]) - np.log(m[0]))[::-1], np.zeros(1)))
+            ngtm_ = np.concatenate((cumulative_trapezoid(dndlnm[::-1], dx=np.log(m[1]) - np.log(m[0]))[::-1], np.zeros(1)))
+            mgtm_ = np.concatenate((cumulative_trapezoid(m[::-1] * dndlnm[::-1], dx=np.log(m[1]) - np.log(m[0]))[::-1], np.zeros(1)))
 
             ngtm.append(ngtm_ + int_upper_n)
             mgtm.append(mgtm_ + int_upper_m)
@@ -567,10 +566,10 @@ def _load_hmf(self):
                 mgtm_0 = np.trapz(self.tab_dndm[i] * self.tab_M**2,
                     x=np.log(self.tab_M))
                 self.tab_ngtm[i,:] = ngtm_0 \
-                    - cumtrapz(self.tab_dndm[i] * self.tab_M,
+                    - cumulative_trapezoid(self.tab_dndm[i] * self.tab_M,
                         x=np.log(self.tab_M), initial=0.0)
                 self.tab_mgtm[i,:] = mgtm_0 \
-                    - cumtrapz(self.tab_dndm[i] * self.tab_M**2,
+                    - cumulative_trapezoid(self.tab_dndm[i] * self.tab_M**2,
                         x=np.log(self.tab_M), initial=0.0)
 
             # Keep it positive please.
diff --git a/ares/physics/RateCoefficients.py b/ares/physics/RateCoefficients.py
index c55fcab5a..2238a6ad0 100755
--- a/ares/physics/RateCoefficients.py
+++ b/ares/physics/RateCoefficients.py
@@ -12,7 +12,7 @@
 """
 
 import numpy as np
-from scipy.misc import derivative
+import numdifftools as nd
 from ..util.Math import interp1d
 from ..util.Math import central_difference
 
@@ -71,7 +71,7 @@ def _dCollisionalIonizationRate(self):
         if not hasattr(self, '_dCollisionalIonizationRate_'):
             self._dCollisionalIonizationRate_ = {}
             for i, absorber in enumerate(self.grid.absorbers):
-                tmp = derivative(lambda T: self.CollisionalIonizationRate(i, T), self.Tarr)
+                tmp = nd.Derivative(lambda T: self.CollisionalIonizationRate(i, T))(self.Tarr)
                 self._dCollisionalIonizationRate_[i] = interp1d(self.Tarr, tmp, 
                     kind=self.interp_rc)
             
@@ -132,7 +132,7 @@ def _dRadiativeRecombinationRate(self):
         if not hasattr(self, '_dRadiativeRecombinationRate_'):
             self._dRadiativeRecombinationRate_ = {}
             for i, absorber in enumerate(self.grid.absorbers):
-                tmp = derivative(lambda T: self.RadiativeRecombinationRate(i, T), self.Tarr)
+                tmp = nd.Derivative(lambda T: self.RadiativeRecombinationRate(i, T))(self.Tarr)
                                 
                 self._dRadiativeRecombinationRate_[i] = interp1d(self.Tarr, tmp, 
                     kind=self.interp_rc)
@@ -161,7 +161,7 @@ def DielectricRecombinationRate(self, T):
     def _dDielectricRecombinationRate(self):
         if not hasattr(self, '_dDielectricRecombinationRate_'):
             self._dDielectricRecombinationRate_ = {}
-            tmp = derivative(lambda T: self.DielectricRecombinationRate(T), self.Tarr)
+            tmp = nd.Derivative(lambda T: self.DielectricRecombinationRate(T))(self.Tarr)
             self._dDielectricRecombinationRate_ = interp1d(self.Tarr, tmp, 
                 kind=self.interp_rc)
     
@@ -197,7 +197,7 @@ def _dCollisionalIonizationCoolingRate(self):
         if not hasattr(self, '_dCollisionalIonizationCoolingRate_'):
             self._dCollisionalIonizationCoolingRate_ = {}
             for i, absorber in enumerate(self.grid.absorbers):
-                tmp = derivative(lambda T: self.CollisionalExcitationCoolingRate(i, T), self.Tarr)
+                tmp = nd.Derivative(lambda T: self.CollisionalExcitationCoolingRate(i, T))(self.Tarr)
                 self._dCollisionalIonizationCoolingRate_[i] = interp1d(self.Tarr, tmp, 
                     kind=self.interp_rc)
     
@@ -233,7 +233,7 @@ def _dCollisionalExcitationCoolingRate(self):
         if not hasattr(self, '_dCollisionalExcitationCoolingRate_'):
             self._dCollisionalExcitationCoolingRate_ = {}
             for i, absorber in enumerate(self.grid.absorbers):
-                tmp = derivative(lambda T: self.CollisionalExcitationCoolingRate(i, T), self.Tarr)
+                tmp = nd.Derivative(lambda T: self.CollisionalExcitationCoolingRate(i, T))(self.Tarr)
                 self._dCollisionalExcitationCoolingRate_[i] = interp1d(self.Tarr, tmp, 
                     kind=self.interp_rc)
     
@@ -272,7 +272,7 @@ def _dRecombinationCoolingRate(self):
         if not hasattr(self, '_dRecombinationCoolingRate_'):
             self._dRecombinationCoolingRate_ = {}
             for i, absorber in enumerate(self.grid.absorbers):
-                tmp = derivative(lambda T: self.RecombinationCoolingRate(i, T), self.Tarr)
+                tmp = nd.Derivative(lambda T: self.RecombinationCoolingRate(i, T))(self.Tarr)
                 self._dRecombinationCoolingRate_[i] = interp1d(self.Tarr, tmp, 
                     kind=self.interp_rc)
     
@@ -300,7 +300,7 @@ def DielectricRecombinationCoolingRate(self, T):
     @property 
     def _dDielectricRecombinationCoolingRate(self):
         if not hasattr(self, '_dDielectricRecombinationCoolingRate_'):
-            tmp = derivative(lambda T: self.DielectricRecombinationCoolingRate(T), self.Tarr)
+            tmp = nd.Derivative(lambda T: self.DielectricRecombinationCoolingRate(T))(self.Tarr)
             self._dDielectricRecombinationCoolingRate_ = interp1d(self.Tarr, tmp, 
                 kind=self.interp_rc)
     
diff --git a/ares/populations/GalaxyAggregate.py b/ares/populations/GalaxyAggregate.py
index 62c3c8c66..0c59b3fc1 100755
--- a/ares/populations/GalaxyAggregate.py
+++ b/ares/populations/GalaxyAggregate.py
@@ -18,7 +18,7 @@
 from .Halo import HaloPopulation
 from collections import namedtuple
 from ..util.Math import interp1d
-from scipy.integrate import quad, simps
+from scipy.integrate import quad, simpson
 from ..util.Warnings import negative_SFRD
 from ..util.ParameterFile import get_pq_pars, pop_id_num
 from scipy.interpolate import interp1d as interp1d_scipy
diff --git a/ares/populations/GalaxyCohort.py b/ares/populations/GalaxyCohort.py
index 08df70d98..7155d6b37 100755
--- a/ares/populations/GalaxyCohort.py
+++ b/ares/populations/GalaxyCohort.py
@@ -13,15 +13,15 @@
 import re
 import time
 import numpy as np
+import numdifftools as nd
 from ..util import read_lit
 from inspect import ismethod
 from types import FunctionType
 from ..util import ProgressBar
 from ..analysis import ModelSet
-from scipy.misc import derivative
 from scipy.optimize import fsolve, minimize
 from ..analysis.BlobFactory import BlobFactory
-from scipy.integrate import quad, simps, cumtrapz, ode
+from scipy.integrate import quad, simpson, cumulative_trapezoid, ode
 from ..util.ParameterFile import par_info, get_pq_pars
 from ..physics.RateCoefficients import RateCoefficients
 from scipy.interpolate import RectBivariateSpline
@@ -493,7 +493,7 @@ def _get_luminosity_density(self, Emin=None, Emax=None):
                 * self.tab_focc[i] * yield_per_sfr(**kw) * ok[i]
 
             _tot = np.trapz(integrand, x=np.log(self.halos.tab_M))
-            _cumtot = cumtrapz(integrand, x=np.log(self.halos.tab_M), initial=0.0)
+            _cumtot = cumulative_trapezoid(integrand, x=np.log(self.halos.tab_M), initial=0.0)
 
             _tmp = _tot - \
                 np.interp(np.log(self._tab_Mmin[i]), np.log(self.halos.tab_M), _cumtot)
@@ -573,7 +573,7 @@ def _get_photon_density(self, z, Emin, Emax):
                 * self.tab_focc[i] * N_per_Msun * fesc * ok[i]
 
             tot = np.trapz(integrand, x=np.log(self.halos.tab_M))
-            cumtot = cumtrapz(integrand, x=np.log(self.halos.tab_M),
+            cumtot = cumulative_trapezoid(integrand, x=np.log(self.halos.tab_M),
                 initial=0.0)
 
             tab[i] = tot - \
@@ -618,7 +618,7 @@ def get_smd(self, z):
 
         if not hasattr(self, '_func_smd'):
             dtdz = np.array([self.cosm.dtdz(z) for z in self.halos.tab_z])
-            self._tab_smd = cumtrapz(self.tab_sfrd_total[-1::-1] * dtdz[-1::-1],
+            self._tab_smd = cumulative_trapezoid(self.tab_sfrd_total[-1::-1] * dtdz[-1::-1],
                 dx=np.abs(np.diff(self.halos.tab_z[-1::-1])), initial=0.)[-1::-1]
             self._func_smd = interp1d(self.halos.tab_z, self._tab_smd,
                 kind=self.pf['pop_interp_sfrd'])
@@ -688,10 +688,10 @@ def _tab_eta(self):
 
                     integ = self.halos.tab_dndlnm[i] * MAR
 
-                    p0 = simps(integ[j1-1:], x=np.log(self.halos.tab_M)[j1-1:])
-                    p1 = simps(integ[j1:], x=np.log(self.halos.tab_M)[j1:])
-                    p2 = simps(integ[j1+1:], x=np.log(self.halos.tab_M)[j1+1:])
-                    p3 = simps(integ[j1+2:], x=np.log(self.halos.tab_M)[j1+2:])
+                    p0 = simpson(integ[j1-1:], x=np.log(self.halos.tab_M)[j1-1:])
+                    p1 = simpson(integ[j1:], x=np.log(self.halos.tab_M)[j1:])
+                    p2 = simpson(integ[j1+1:], x=np.log(self.halos.tab_M)[j1+1:])
+                    p3 = simpson(integ[j1+2:], x=np.log(self.halos.tab_M)[j1+2:])
 
                     interp = interp1d(np.log(self.halos.tab_M)[j1-1:j1+3], [p0,p1,p2,p3],
                         kind=self.pf['pop_interp_MAR'])
@@ -991,7 +991,7 @@ def get_surface_density(self, z, mag=None, dz=1., dtheta=1., wave=1600.):
         Ngal = Ngal[-1::-1]
 
         # Cumulative surface density of galaxies *brighter than* Mobs
-        cgal = cumtrapz(Ngal, x=Mobs, initial=Ngal[0])
+        cgal = cumulative_trapezoid(Ngal, x=Mobs, initial=Ngal[0])
 
         if mag is not None:
             return np.interp(mag, Mobs, cgal)
@@ -1013,7 +1013,7 @@ def get_surface_density(self, z, mag=None, dz=1., dtheta=1., wave=1600.):
         # some corresponding magnitude
         assert Ngal[0] == 0, "Broaden binning range?"
         ntot = np.trapz(Ngal, x=x)
-        nltm = cumtrapz(Ngal, x=x, initial=Ngal[0])
+        nltm = cumulative_trapezoid(Ngal, x=x, initial=Ngal[0])
 
         return x, nltm
 
@@ -2004,7 +2004,7 @@ def tab_nh_active(self):
                         i1 -= 1
                     i2 = i1 + 1
 
-                    # Trapezoid here we come
+                    # trapz here we come
                     b = self._tab_logMmax[i] - self._tab_logMmin[i]
 
                     M1 = self._tab_logMmin[i]
@@ -2093,7 +2093,7 @@ def tab_sfrd_total(self):
             integrand = self.tab_sfr * self.halos.tab_dndlnm * self.tab_focc
 
             ##
-            # Use cumtrapz instead and interpolate onto Mmin, Mmax
+            # Use cumulative_trapezoid instead and interpolate onto Mmin, Mmax
             ##
             ct = 0
             self._tab_sfrd_total_ = np.zeros_like(self.halos.tab_z)
@@ -2122,7 +2122,7 @@ def tab_sfrd_total(self):
                 #        continue
 
                 tot = np.trapz(integrand[i], x=np.log(self.halos.tab_M))
-                cumtot = cumtrapz(integrand[i], x=np.log(self.halos.tab_M),
+                cumtot = cumulative_trapezoid(integrand[i], x=np.log(self.halos.tab_M),
                     initial=0.0)
 
                 above_Mmin = np.interp(np.log(self._tab_Mmin[i]),
@@ -2393,7 +2393,7 @@ def get_sfe_slope(self, z, Mh):
 
         logfst = lambda logM: np.log10(self.SFE(z=z, Mh=10**logM))
 
-        return derivative(logfst, np.log10(Mh), dx=0.01)[0]
+        return nd.Derivative(logfst)(np.log10(Mh))
 
     @property
     def _tab_Mz(self):
@@ -2622,7 +2622,7 @@ def _SAM_1z_jac(self, z, y): # pragma: no cover
 
         # Eq. 1: halo mass.
         _y1p = lambda _Mh: self.MGR(z, _Mh) * dtdz
-        y1p = derivative(_y1p, Mh)
+        y1p = nd.Derivative(_y1p)(Mh)
 
         # Eq. 2: gas mass
         if self.pf['pop_sfr'] is None:
diff --git a/ares/populations/GalaxyEnsemble.py b/ares/populations/GalaxyEnsemble.py
index 00771bf12..832cbe47f 100755
--- a/ares/populations/GalaxyEnsemble.py
+++ b/ares/populations/GalaxyEnsemble.py
@@ -25,7 +25,7 @@
 from scipy.optimize import curve_fit
 from .GalaxyCohort import GalaxyCohort
 from scipy.interpolate import interp1d
-from scipy.integrate import quad, cumtrapz
+from scipy.integrate import quad, cumulative_trapezoid
 from ..analysis.BlobFactory import BlobFactory
 from ..obs.Photometry import get_filters_from_waves
 from ..util.Stats import bin_e2c, bin_c2e, bin_samples, quantify_scatter
@@ -1062,7 +1062,7 @@ def _gen_prescribed_galaxy_histories(self, zstop=0):
         # OK. We've got a bunch of halo histories and we need to integrate them
         # to get things like stellar mass, metal mass, etc. This means we need
         # to make an assumption about how halos grow between our grid points.
-        # If we assume smooth histories, we should be trying to do a trapezoidal
+        # If we assume smooth histories, we should be trying to do a trapzal
         # integration in log-space. However, given that we often add noise to
         # MARs, in that case we should probably just assume flat MARs within
         # the timestep.
@@ -1094,7 +1094,7 @@ def _gen_prescribed_galaxy_histories(self, zstop=0):
         fml = (1. - fmr)
 
         # Integrate (crudely) mass accretion rates
-        #_Mint = cumtrapz(_MAR[:,:], dx=dt, axis=1)
+        #_Mint = cumulative_trapezoid(_MAR[:,:], dx=dt, axis=1)
         #_MAR_c = 0.5 * (np.roll(MAR, -1, axis=1) + MAR)
         #_Mint = np.cumsum(_MAR_c[:,1:] * dt, axis=1)
 
@@ -2641,7 +2641,7 @@ def get_bias(self, z, limit=None, wave=1600., cam=None, filters=None,
         integ_top = bh[ok==1] * _nh[ok==1]
         integ_bot = _nh[ok==1]
 
-        # Don't do trapz -- not a continuous curve like in GalaxyCohort.
+        # Don't dotrapz -- not a continuous curve like in GalaxyCohort.
         b = np.sum(integ_top * _Mh[ok==1]) / np.sum(integ_bot * _Mh[ok==1])
 
         return b
@@ -3676,7 +3676,7 @@ def get_surface_density(self, z, bins=None, dz=1., dtheta=1., wave=1600.,
             # some corresponding magnitude
             assert Ngal[i,0] == 0, "Broaden binning range?"
             #ntot = np.trapz(Ngal[i,:], x=x)
-            nltm[i,:] = cumtrapz(Ngal[i,:], x=x, initial=Ngal[i,0])
+            nltm[i,:] = cumulative_trapezoid(Ngal[i,:], x=x, initial=Ngal[i,0])
 
         # Can just return *maximum* number of galaxies detected,
         # regardless of band. Equivalent to requiring only single-band
diff --git a/ares/populations/Halo.py b/ares/populations/Halo.py
index efad9f8a6..7f3a3c280 100755
--- a/ares/populations/Halo.py
+++ b/ares/populations/Halo.py
@@ -15,7 +15,7 @@
 from inspect import ismethod
 from types import FunctionType
 from .Population import Population
-from scipy.integrate import cumtrapz
+from scipy.integrate import cumulative_trapezoid
 from ..util.PrintInfo import print_pop
 from scipy.interpolate import interp1d
 from ..physics.HaloModel import HaloModel
diff --git a/ares/solvers/UniformBackground.py b/ares/solvers/UniformBackground.py
index b8c1a83e9..d5041f045 100755
--- a/ares/solvers/UniformBackground.py
+++ b/ares/solvers/UniformBackground.py
@@ -23,7 +23,7 @@
 from ..physics import Hydrogen, Cosmology
 from ..populations.Composite import CompositePopulation
 from ..populations.GalaxyAggregate import GalaxyAggregate
-from scipy.integrate import quad, romberg, romb, trapz, simps
+from scipy.integrate import quad, trapezoid, simpson
 from ..physics.Constants import ev_per_hz, erg_per_ev, c, E_LyA, E_LL, dnu, h_p
 #from ..util.ReadData import flatten_energies, flatten_flux, split_flux, \
 #    flatten_emissivities
@@ -973,7 +973,7 @@ def TabulateEmissivity(self, z, E, pop):
         # Special case: delta function SED! Can't normalize a-priori without
         # knowing binning, so we do it here.
         if pop.src.is_delta:
-            # This is a little weird. Trapezoidal integration doesn't make
+            # This is a little weird. trapzal integration doesn't make
             # sense for a delta function, but it's what happens later, so
             # insert a factor of a half now so we recover all the flux we
             # should.
diff --git a/ares/sources/SynthesisModel.py b/ares/sources/SynthesisModel.py
index ea98abc04..d4f5dc6a0 100644
--- a/ares/sources/SynthesisModel.py
+++ b/ares/sources/SynthesisModel.py
@@ -16,7 +16,7 @@
 from ..util.Math import interp1d
 from ares.physics import Cosmology
 from scipy.optimize import minimize
-from scipy.integrate import cumtrapz
+from scipy.integrate import cumulative_trapezoid
 from ..util.ReadData import read_lit
 from ..physics import NebularEmission
 from ..util.ParameterFile import ParameterFile
@@ -626,7 +626,7 @@ def PhotonsPerBaryon(self, Emin, Emax, raw=True, return_all_t=False):
         if self.pf['source_ssp']:
             photons_per_b_t = photons_per_s_per_msun / self.cosm.b_per_msun
             if return_all_t:
-                return cumtrapz(photons_per_b_t, x=self.times*s_per_myr,
+                return cumulative_trapezoid(photons_per_b_t, x=self.times*s_per_myr,
                     initial=0.0)
             else:
                 return np.trapz(photons_per_b_t, x=self.times*s_per_myr)
diff --git a/ares/sources/SynthesisModelSBS.py b/ares/sources/SynthesisModelSBS.py
index 65d4c0ee3..ebc8141ab 100644
--- a/ares/sources/SynthesisModelSBS.py
+++ b/ares/sources/SynthesisModelSBS.py
@@ -15,7 +15,7 @@
 from .Source import Source
 import matplotlib.pyplot as pl
 from ..util.ReadData import read_lit
-from scipy.integrate import quad, cumtrapz
+from scipy.integrate import quad, cumulative_trapezoid
 from ..util.Stats import bin_c2e, bin_e2c
 from ..util.ParameterFile import ParameterFile
 from ..physics.Constants import h_p, c, erg_per_ev, ev_per_hz, \
@@ -348,7 +348,7 @@ def ngtm(self, m):
         return 1. - 10**np.interp(np.log10(m), np.log10(self.Ms), np.log10(self.tab_imf_cdf))
         
     def mgtm(self, m):
-        cdf_by_m = cumtrapz(self.tab_imf * self.Ms**2, x=np.log(self.Ms), initial=0.) \
+        cdf_by_m = cumulative_trapezoid(self.tab_imf * self.Ms**2, x=np.log(self.Ms), initial=0.) \
             / np.trapz(self.tab_imf * self.Ms**2, x=np.log(self.Ms))
         
         return 1. - np.interp(m, self.Ms, cdf_by_m)
@@ -439,7 +439,7 @@ def tab_imf_cdf(self):
                 self._tab_imf_cdf = np.concatenate((_tot0, _tot1, _tot2)) / _tot2[-1]
                 
             else:    
-                self._tab_imf_cdf = cumtrapz(self.tab_imf, x=self.Ms, initial=0.) \
+                self._tab_imf_cdf = cumulative_trapezoid(self.tab_imf, x=self.Ms, initial=0.) \
                     / np.trapz(self.tab_imf * self.Ms, x=np.log(self.Ms))
 
         return self._tab_imf_cdf
@@ -491,7 +491,7 @@ def avg_sn_delay(self):
     def tab_dtd_cdf(self):
         if not hasattr(self, '_tab_dtd_cdf'):
             ok = self.Ms >= 8.
-            top = cumtrapz(self.tab_life[ok==1] * self.tab_imf[ok==1] \
+            top = cumulative_trapezoid(self.tab_life[ok==1] * self.tab_imf[ok==1] \
                 * self.Ms[ok==1], x=np.log(self.Ms[ok==1]), initial=0.0)
                 
             bot = np.trapz(self.tab_life[ok==1] * self.tab_imf[ok==1] \
diff --git a/ares/static/ChemicalNetwork.py b/ares/static/ChemicalNetwork.py
index 81e0eb706..6eb8ab1d1 100755
--- a/ares/static/ChemicalNetwork.py
+++ b/ares/static/ChemicalNetwork.py
@@ -13,7 +13,7 @@
 
 import copy, sys
 import numpy as np
-from scipy.misc import derivative
+import numdifftools as nd
 from ..util.Warnings import solver_error
 from ..physics.RateCoefficients import RateCoefficients
 from ..physics.Constants import k_B, sigma_T, m_e, c, s_per_myr, erg_per_ev, h
@@ -639,8 +639,7 @@ def Jacobian(self, t, q, args): # pragma: no cover
 
         # Add in any parametric modifications?
         if self.exotic_heating:
-            J[-1,-1] += derivative(self.grid._exotic_func(z=z) * to_temp, z,
-                dx=0.05)
+            J[-1,-1] += nd.Derivative(lambda z: self.grid._exotic_func(z) * to_temp)(z)
 
         return J
 
diff --git a/ares/static/Fluctuations.py b/ares/static/Fluctuations.py
index 0f79517c4..b7697a4c3 100644
--- a/ares/static/Fluctuations.py
+++ b/ares/static/Fluctuations.py
@@ -18,7 +18,7 @@
 from scipy.special import erfinv
 from scipy.optimize import fsolve
 from scipy.interpolate import interp1d
-from scipy.integrate import quad, simps
+from scipy.integrate import quad, simpson
 from ..physics.Hydrogen import Hydrogen
 from ..physics.HaloModel import HaloModel
 from ..populations.Composite import CompositePopulation
@@ -459,8 +459,8 @@ def mean_halo_bias(self, z):
 
         return 1.0
 
-        #return simps(M_h * dndm_h * bias, x=np.log(M_h)) \
-        #    / simps(M_h * dndm_h, x=np.log(M_h))
+        #return simpson(M_h * dndm_h * bias, x=np.log(M_h)) \
+        #    / simpson(M_h * dndm_h, x=np.log(M_h))
 
     def tab_bubble_bias(self, zeta):
         if not hasattr(self, '_tab_bubble_bias'):
@@ -569,7 +569,7 @@ def mean_bubble_bias(self, z, ion=True):
     #
     #    dm_ddel = rho0 * V_i
     #
-    #    return simps(B[iM:] * dndm_b[iM:] * M_b[iM:], x=np.log(M_b[iM:]))
+    #    return simpson(B[iM:] * dndm_b[iM:] * M_b[iM:], x=np.log(M_b[iM:]))
 
     def delta_bubble_vol_weighted(self, z, ion=True):
         if not self.pf['ps_include_ion']:
@@ -610,7 +610,7 @@ def delta_bubble_vol_weighted(self, z, ion=True):
    #
    #    dm_ddel = rho0 * V_i
    #
-   #    return simps(B[iM:] * dndm_b[iM:] * M_b[iM:], x=np.log(M_b[iM:]))
+   #    return simpson(B[iM:] * dndm_b[iM:] * M_b[iM:], x=np.log(M_b[iM:]))
 
     def mean_halo_abundance(self, z, Mmin=False):
         M_h = self.halos.tab_M
@@ -1911,7 +1911,7 @@ def ExpectationValue2pt(self, z, R, term='ii', R_s=None, R3=None,
                     #Vii = all_V[0]
                     #_integrand1 = dndm * Vii
                     #
-                    #_exp_int1 = np.exp(-simps(_integrand1[iM:] * M_b[iM:],
+                    #_exp_int1 = np.exp(-simpson(_integrand1[iM:] * M_b[iM:],
                     #    x=np.log(M_b[iM:])))
                     #_P1_ii = (1. - _exp_int1)
 
diff --git a/ares/static/IntegralTables.py b/ares/static/IntegralTables.py
index b6305a585..0e9ae2005 100755
--- a/ares/static/IntegralTables.py
+++ b/ares/static/IntegralTables.py
@@ -17,7 +17,8 @@
 from ..physics.Constants import erg_per_ev
 from ..physics.SecondaryElectrons import *
 import os, re, scipy, itertools, math, copy
-from scipy.integrate import quad, trapz, simps
+# changed trapz to trapezoid (update)
+from scipy.integrate import quad, trapezoid, simpson
 
 try:
     from mpi4py import MPI
@@ -701,7 +702,7 @@ def PhiHat(self, N, absorber, donor=None, x=None, t=None, ind=None):
             c &= self.E <= self.src.Emax
             samples = np.array([integrand(E) for E in self.E[c]])[..., 0]
 
-            integral = simps(samples, self.E[c]) / erg_per_ev
+            integral = simpson(samples, self.E[c]) / erg_per_ev
 
         if not self.pf['photon_conserving']:
             integral *= self.E_th[absorber]
@@ -738,7 +739,7 @@ def PsiHat(self, N, absorber, donor=None, x=None, t=None):
             c &= self.E <= self.src.Emax
             samples = np.array([integrand(E) for E in self.E[c]])[..., 0]
 
-            integral = simps(samples, self.E[c])
+            integral = simpson(samples, self.E[c])
 
         if not self.pf['photon_conserving']:
             integral *= self.E_th[absorber]
@@ -776,7 +777,7 @@ def PhiWiggle(self, N, absorber, donor, x=None, t=None):
             c &= self.E <= self.src.Emax
             samples = np.array([integrand(E) for E in self.E[c]])[..., 0]
 
-            integral = simps(samples, self.E[c]) / erg_per_ev
+            integral = simpson(samples, self.E[c]) / erg_per_ev
 
         if not self.pf['photon_conserving']:
             integral *= self.E_th[absorber]
@@ -811,7 +812,7 @@ def PsiWiggle(self, N, absorber, donor, x=None, t=None):
             c &= self.E <= self.src.Emax
             samples = np.array([integrand(E) for E in self.E[c]])[..., 0]
 
-            integral = simps(samples, self.E[c])
+            integral = simpson(samples, self.E[c])
 
         if not self.pf['photon_conserving']:
             integral *= self.E_th[absorber]
diff --git a/ares/static/SpectralSynthesis.py b/ares/static/SpectralSynthesis.py
index b2fdb288e..c2c997e7b 100644
--- a/ares/static/SpectralSynthesis.py
+++ b/ares/static/SpectralSynthesis.py
@@ -1320,7 +1320,7 @@ def get_lum(self, wave=1600., sfh=None, tarr=None, zarr=None,
             ###
             ## Integrate over all times up to this tobs
             if batch_mode:
-                # Should really just np.sum here...using trapz assumes that
+                # Should really just np.sum here...usingtrapz assumes that
                 # the SFH is a smooth function and not a series of constant
                 # SFRs. Doesn't really matter in practice, though.
                 if not do_all_time:
diff --git a/ares/static/VolumeGlobal.py b/ares/static/VolumeGlobal.py
index 5b3a7ade0..814f05872 100755
--- a/ares/static/VolumeGlobal.py
+++ b/ares/static/VolumeGlobal.py
@@ -16,7 +16,8 @@
 from ..physics.Constants import *
 import types, os, re, sys
 from ..physics import SecondaryElectrons
-from scipy.integrate import dblquad, romb, simps, quad, trapz
+# changed trapz to trapezoid (update)
+from scipy.integrate import dblquad, romb, simpson, quad, trapezoid
 
 try:
     import h5py
@@ -526,7 +527,7 @@ def HeatingRate(self, z, species=0, popid=0, band=0, **kwargs):
                     heat = romb(integrand[0:imax] * self.E[0:imax],
                         dx=self.dlogE[0:imax])[0] * log10
                 else:
-                    heat = simps(integrand[0:imax] * self._E[popid][band][0:imax],
+                    heat = simpson(integrand[0:imax] * self._E[popid][band][0:imax],
                         x=self.logE[popid][band][0:imax]) * log10
 
             else:
@@ -539,7 +540,7 @@ def HeatingRate(self, z, species=0, popid=0, band=0, **kwargs):
                     heat = np.trapz(integrand[imin:] * self._E[popid][band][imin:],
                         x=self.logE[popid][band][imin:]) * log10
                 else:
-                    heat = simps(integrand[imin:] * self._E[popid][band][imin:],
+                    heat = simpson(integrand[imin:] * self._E[popid][band][imin:],
                         x=self.logE[popid][band][imin:]) * log10
 
         # Re-normalize, get rid of per steradian units
@@ -692,7 +693,7 @@ def IonizationRateIGM(self, z, species=0, popid=0, band=0, **kwargs):
                 ion = romb(integrand * self.E[popid][band],
                     dx=self.dlogE[popid][band])[0] * log10
             else:
-                ion = simps(integrand * self.E[popid][band],
+                ion = simpson(integrand * self.E[popid][band],
                     x=self.logE[popid][band]) * log10
 
         # Re-normalize
@@ -827,7 +828,7 @@ def SecondaryIonizationRateIGM(self, z, species=0, donor=0, popid=0,
                 ion = romb(integrand * self.E[popid][band],
                     dx=self.dlogE[popid][band])[0] * log10
             else:
-                ion = simps(integrand * self.E[popid][band],
+                ion = simpson(integrand * self.E[popid][band],
                     x=self.logE[popid][band]) * log10
 
         # Re-normalize
@@ -923,7 +924,7 @@ def SecondaryLymanAlphaFlux(self, z, species=0, popid=0, band=0,
                 e_ax = romb(integrand[0:imax] * self.E[0:imax],
                     dx=self.dlogE[0:imax])[0] * log10
             else:
-                e_ax = simps(integrand[0:imax] * self._E[popid][band][0:imax],
+                e_ax = simpson(integrand[0:imax] * self._E[popid][band][0:imax],
                     x=self.logE[popid][band][0:imax]) * log10
         else:
             imin = np.argmin(np.abs(self._E[popid][band] - pop.pf['pop_Emin']))
@@ -935,7 +936,7 @@ def SecondaryLymanAlphaFlux(self, z, species=0, popid=0, band=0,
                 e_ax = np.trapz(integrand[imin:] * self._E[popid][band][imin:],
                     x=self.logE[popid][band][imin:]) * log10
             else:
-                e_ax = simps(integrand[imin:] * self._E[popid][band][imin:],
+                e_ax = simpson(integrand[imin:] * self._E[popid][band][imin:],
                     x=self.logE[popid][band][imin:]) * log10
 
         # Re-normalize. This is essentially a photon emissivity modulo 4 pi ster
diff --git a/ares/util/Aesthetics.py b/ares/util/Aesthetics.py
index 228fa4afd..5e2ac966c 100755
--- a/ares/util/Aesthetics.py
+++ b/ares/util/Aesthetics.py
@@ -9,8 +9,8 @@
 Description:
 
 """
-
-import os, imp, re
+# relplaced imp with importlib using this guidance: https://docs.python.org/3/whatsnew/3.12.html#whatsnew312-removed-imp
+import os, importlib, re
 import numpy as np
 from matplotlib import cm
 from .ParameterFile import par_info
@@ -36,9 +36,9 @@
 
 # Load custom defaults
 HOME = os.environ.get('HOME')
-if os.path.exists('{!s}/.ares/labels.py'.format(HOME)):
-    f, filename, data = imp.find_module('labels', ['{!s}/.ares/'.format(HOME)])
-    custom_labels = imp.load_module('labels.py', f, filename, data).pf
+if os.path.exists('~/.ares/labels.py'.format(HOME)):
+    f, filename, data = importlib.util.find_spec('labels', ['~/.ares/'.format(HOME)])
+    custom_labels = importlib.import_module('labels.py', f, filename, data).pf
 else:
     custom_labels = {}
 
diff --git a/ares/util/ReadData.py b/ares/util/ReadData.py
index 49d1d0024..0530c4c71 100755
--- a/ares/util/ReadData.py
+++ b/ares/util/ReadData.py
@@ -9,9 +9,9 @@
 Description:
 
 """
-
+# replaced outdated imp module with importlib according to the following guidance: https://docs.python.org/3/whatsnew/3.12.html#whatsnew312-removed-imp
 import numpy as np
-import imp as _imp
+import importlib as _imp
 from ..data import ARES
 import os, re, sys, glob
 from .Pickling import read_pickle_file
@@ -71,8 +71,8 @@ def read_lit(prefix, path=None, verbose=True):
         print("WARNING: multiple copies of {!s} found.".format(prefix))
         print("       : precedence: CWD -> $HOME -> $ARES/input/litdata")
 
-    _f, _filename, _data = _imp.find_module(prefix, [loc])
-    mod = _imp.load_module(prefix, _f, _filename, _data)
+    _f, _filename, _data = _imp.util.find_spec(prefix, [loc])
+    mod = _imp.import_module(prefix, _f, _filename, _data)
 
     # Save this for sanity checks later
     mod.path = loc
diff --git a/ares/util/SetDefaultParameterValues.py b/ares/util/SetDefaultParameterValues.py
index 9b53e2ede..eee50672d 100644
--- a/ares/util/SetDefaultParameterValues.py
+++ b/ares/util/SetDefaultParameterValues.py
@@ -8,8 +8,8 @@
 Description: Defaults for all different kinds of parameters.
 
 """
-
-import os, imp
+# removed unnecessary imp call (outdated module)
+import os
 import numpy as np
 from ..data import ARES
 from ares import rcParams
@@ -1451,7 +1451,7 @@ def ControlParameters():
     # Real-time optical depth calculation once EoR begins
     "EoR_xavg": 1.0,        # ionized fraction indicating start of EoR (OFF by default)
     "EoR_dlogx": 0.001,
-    "EoR_approx_tau": False, # 0 = trapezoidal integration,
+    "EoR_approx_tau": False, # 0 = trapzal integration,
                              # 1 = mean ionized fraction, approx cross sections
                              # 2 = neutral approx, approx cross sections
 
@@ -1480,7 +1480,7 @@ def ControlParameters():
     "sed_prefix": None,
 
     "unsampled_integrator": 'quad',
-    "sampled_integrator": 'simps',
+    "sampled_integrator": 'simpson',
     "integrator_rtol": 1e-6,
     "integrator_atol": 1e-4,
     "integrator_divmax": 1e2,
diff --git a/input/hmf/test_fcoll.py b/input/hmf/test_fcoll.py
index 12e61a75f..635fe6c33 100644
--- a/input/hmf/test_fcoll.py
+++ b/input/hmf/test_fcoll.py
@@ -13,7 +13,7 @@
 import ares
 import numpy as np
 import matplotlib.pyplot as pl
-from scipy.integrate import simps
+from scipy.integrate import simpson
 
 pop = ares.populations.HaloPopulation(pop_sfr_model='fcoll', pop_Mmin=1e8,
     hmf_interp='linear')
@@ -38,7 +38,7 @@
     dndlnm = dndm * M
     
     #fcoll_mgtm2 = np.trapz(dndlnm, x=np.log(M)) / pop.halos.MF.mean_density0
-    fcoll_mgtm2 = simps(dndlnm[ok], x=np.log(M[ok])) / pop.halos.MF.mean_density0
+    fcoll_mgtm2 = simpson(dndlnm[ok], x=np.log(M[ok])) / pop.halos.MF.mean_density0
         
     print('{0!s} {1!s} {2!s}'.format(z, fcoll_mgtm1, fcoll_mgtm2))#, fcoll
     
diff --git a/input/litdata/leitherer1999.py b/input/litdata/leitherer1999.py
index 44b322a30..5629f0e89 100755
--- a/input/litdata/leitherer1999.py
+++ b/input/litdata/leitherer1999.py
@@ -12,7 +12,7 @@
 import numpy as np
 from ares.data import ARES
 from ares.physics import Cosmology
-from scipy.integrate import cumtrapz
+from scipy.integrate import cumulative_trapezoid
 from scipy.interpolate import interp1d, RectBivariateSpline
 from ares.physics.Constants import h_p, c, erg_per_ev, g_per_msun, s_per_yr, \
     s_per_myr, m_H
diff --git a/requirements.txt b/requirements.txt
index 9c436d14d..6f6c4acc4 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -11,4 +11,5 @@ docutils==0.17.1
 cached_property>=1.5.2<2.0
 camb>=1.3<2.0
 hmf>=3.1<4.0
+numdifftools<1
 -e git+https://bitbucket.org/ktausch/distpy/#egg=distpy