diff --git a/.github/workflows/ci_test.yml b/.github/workflows/ci_test.yml
index c3a8af4c5..79d3ca0a2 100644
--- a/.github/workflows/ci_test.yml
+++ b/.github/workflows/ci_test.yml
@@ -24,16 +24,16 @@ jobs:
       matrix:
         include:
           - os: ubuntu-latest
-            python: '3.12'
+            python: '3.13'
             tox_env: 'ephemeris_connection'
           - os: ubuntu-latest
-            python: '3.12'
+            python: '3.13'
             tox_env: 'black'
           - os: ubuntu-latest
-            python: '3.12'
-            tox_env: 'py312-test-cov'
+            python: '3.13'
+            tox_env: 'py313-test-cov'
           - os: ubuntu-latest
-            python: '3.12'
+            python: '3.13'
             tox_env: 'notebooks'
 #          - os: ubuntu-latest
 #            python: '3.8'
@@ -42,13 +42,13 @@ jobs:
 #            python: '3.10'
 #            tox_env: 'py310-test-alldeps-cov'
           - os: macos-12
-            python: '3.12'
-            tox_env: 'py312-test'
+            python: '3.13'
+            tox_env: 'py313-test'
 #          - os: windows-latest
 #            python: '3.8'
 #            tox_env: 'py38-test'
           - os: ubuntu-latest
-            python: '3.8'
+            python: '3.9'
             tox_env: 'oldestdeps'
 #          - os: ubuntu-latest
 #            python: '3.8'
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 6c2c4b8cc..82fcce346 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -6,6 +6,6 @@ repos:
     -   id: check-merge-conflict
     -   id: check-symlinks
 -   repo: https://github.com/psf/black
-    rev: 24.2.0
+    rev: 24.10.0
     hooks:
     -   id: black
diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md
index 2899b99f8..e5ca814a7 100644
--- a/CHANGELOG-unreleased.md
+++ b/CHANGELOG-unreleased.md
@@ -9,6 +9,16 @@ the released changes.
 
 ## Unreleased
 ### Changed
+- Command line scripts now automatically do `allow_tcb` and `allow_T2` while reading par files.
+- Updated the `plot_chains` function in `event_optimize` so that the subplots are a fixed size to prevent the subplots from being condensed in the case of many fit parameters.
 ### Added
+- Added an option `nbin` to `photonphase` to decide how many phase bins to use for the phaseogram
+- Added an option `linearize_model` to speed up the photon phases calculation within `event_optimize` through the designmatrix.
+- Added AIC and BIC calculation to be written in the post fit parfile from `event_optimize`
+- When TCB->TDB conversion info is missing, will print parameter name
+- `add_param` returns the name of the parameter (useful for numbered parameters)
 ### Fixed
+- Changed WAVE_OM units from 1/d to rad/d.
+- When EQUAD is created from TNEQ, has proper TCB->TDB conversion info
+- TOA selection masks will work when only TOA is the first one
 ### Removed
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3c905242f..3f24dffe3 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,12 @@ and this project, at least loosely, adheres to [Semantic Versioning](https://sem
 This file contains the released changes to the codebase. See CHANGELOG-unreleased.md for
 the unreleased changes. This file should only be changed while tagging a new version.
 
+## [1.1] 2024-11-05
+### Changed
+* Bump oldest python to 3.9
+* Change oldest dependencies: `numpy` 1.18.5 to 1.23.0; `astropy` 4.0 to 5.0.5; `scipy` 1.4.1 to 1.9.0; `matplotlib` 3.2.0 to 3.4.3
+* Update CI testing to use python 3.13
+
 ## [1.0.2] 2024-10-18
 ### Changed
 - Moved the events -> TOAs and photon weights code into the function `load_events_weights` within `event_optimize`.
diff --git a/docs/installation.rst b/docs/installation.rst
index eb780c9ea..c6e23b549 100644
--- a/docs/installation.rst
+++ b/docs/installation.rst
@@ -14,7 +14,7 @@ is more complicated (but not too much).
 Prerequisites
 -------------
 
-PINT requires Python 3.8+ [1]_
+PINT requires Python 3.9+ [1]_
 
 Your Python must have the package installation tool pip_ installed.  Also make sure your ``setuptools`` are up to date (e.g. ``pip install -U setuptools``).
 
diff --git a/pyproject.toml b/pyproject.toml
index 8b85722d9..fa5ebdbe7 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -19,22 +19,22 @@ classifiers = [
     "Operating System :: OS Independent",
     "Programming Language :: Python",
     "Programming Language :: Python :: 3",
-    "Programming Language :: Python :: 3.8",
     "Programming Language :: Python :: 3.9",
     "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
+    "Programming Language :: Python :: 3.13",
     "Topic :: Scientific/Engineering :: Astronomy",
     "Topic :: Software Development :: Libraries :: Python Modules",
 ]
-requires-python = ">=3.8"
+requires-python = ">=3.9"
 dependencies = [
-    "numpy>=1.18.5",
-    "astropy>=4.0,!=4.0.1,!=4.0.1.post1",
+    "numpy>=1.23.0",
+    "astropy>=5.0.5",
     "pyerfa",
-    "scipy>=1.4.1",
+    "scipy>=1.9.0",
     "jplephem>=2.6",
-    "matplotlib>=3.2.0",
+    "matplotlib>=3.4.3",
     "emcee>=3.0.1",
     "corner>=2.0.1",
     "uncertainties",
diff --git a/requirements.txt b/requirements.txt
index 16f76401e..293e3d07a 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,9 +1,9 @@
-numpy>=1.18.5
-astropy>=4.0,!=4.0.1,!=4.0.1.post1
+numpy>=1.23.0
+astropy>=5.0.5
 pyerfa
-scipy>=1.4.1
+scipy>=1.9.0
 jplephem>=2.6
-matplotlib>=3.2.0
+matplotlib>=3.4.3
 emcee>=3.0.1
 corner>=2.0.1
 uncertainties
diff --git a/requirements_dev.txt b/requirements_dev.txt
index a00511840..f5dfff94b 100644
--- a/requirements_dev.txt
+++ b/requirements_dev.txt
@@ -1,10 +1,9 @@
 attrs>=19.2
 pip>=9.0.1
-setuptools>=41.0
+setuptools>=61.2
 coverage>=4.3.4
 traitlets
 Sphinx==6.2.1
-astropy>=4.0,!=4.0.1,!=4.0.1.post1
 #astropy-helpers>=1.3
 sphinx-rtd-theme==1.2.2
 coveralls>=1.1
@@ -29,8 +28,7 @@ pdbpp
 tox
 pre-commit
 typed-ast>=1.5.0
-#black>=19.0a0,<20.0a0
-black~=23.0
+black~=24.0
 pygments
 ipython
 pathlib2
diff --git a/src/pint/__init__.py b/src/pint/__init__.py
index c0d511e7c..7906be01f 100644
--- a/src/pint/__init__.py
+++ b/src/pint/__init__.py
@@ -16,7 +16,6 @@
 import astropy.time as time
 import astropy.units as u
 import numpy as np
-import pkg_resources
 from astropy.units import si
 from pathlib import Path
 
@@ -106,6 +105,9 @@
     "hourangle_second": hourangle_second,
 }
 
+# define a units equivalency for gauss in cgs
+gauss_equiv = [u.Gauss, u.Hz * (u.g / u.cm) ** (1 / 2), lambda x: x, lambda x: x]
+
 import astropy.version
 
 if astropy.version.major < 4:
diff --git a/src/pint/config.py b/src/pint/config.py
index a24bec249..344bc71b1 100644
--- a/src/pint/config.py
+++ b/src/pint/config.py
@@ -1,7 +1,7 @@
 """Functions related to PINT configuration."""
 
 import os
-import pkg_resources
+import importlib.resources
 
 __all__ = ["datadir", "examplefile", "runtimefile"]
 
@@ -15,7 +15,7 @@ def datadir() -> str:
     str
         Directory of PINT data files
     """
-    return pkg_resources.resource_filename(__name__, "data/")
+    return os.path.join(importlib.resources.files("pint"), "data/")
 
 
 def examplefile(filename: str) -> str:
@@ -35,9 +35,7 @@ def examplefile(filename: str) -> str:
     This is **not** for files needed at runtime. Those are located by :func:`pint.config.runtimefile`.
     This is for files needed for the example notebooks.
     """
-    return pkg_resources.resource_filename(
-        __name__, os.path.join("data/examples/", filename)
-    )
+    return os.path.join(importlib.resources.files("pint"), f"data/examples/{filename}")
 
 
 def runtimefile(filename: str) -> str:
@@ -57,6 +55,4 @@ def runtimefile(filename: str) -> str:
     This **is**  for files needed at runtime. Files needed for the example notebooks
     are found via :func:`pint.config.examplefile`.
     """
-    return pkg_resources.resource_filename(
-        __name__, os.path.join("data/runtime/", filename)
-    )
+    return os.path.join(importlib.resources.files("pint"), f"data/runtime/{filename}")
diff --git a/src/pint/derived_quantities.py b/src/pint/derived_quantities.py
index d5852ae4b..24dcb1934 100644
--- a/src/pint/derived_quantities.py
+++ b/src/pint/derived_quantities.py
@@ -136,6 +136,14 @@ def pferrs(
     return (forp, forperr, fdorpd, fdorpderr)
 
 
+def _to_gauss(B: u.Quantity) -> u.G:
+    """Convert quantity with mass, length, and time units to Gauss.
+
+    In cgs units, magnetic field is has units (mass/length)^(1/2) / time.
+    """
+    return B.to(u.Gauss, equivalencies=[pint.gauss_equiv])
+
+
 @u.quantity_input(f=u.Hz, fdot=u.Hz / u.s, fo=u.Hz)
 def pulsar_age(
     f: u.Quantity, fdot: u.Quantity, n: int = 3, fo: u.Quantity = 1e99 * u.Hz
@@ -219,12 +227,16 @@ def pulsar_edot(
     return (-4.0 * np.pi**2 * I * f * fdot).to(u.erg / u.s)
 
 
-@u.quantity_input(f=u.Hz, fdot=u.Hz / u.s)
-def pulsar_B(f: u.Quantity, fdot: u.Quantity) -> u.G:
+@u.quantity_input(f=u.Hz, fdot=u.Hz / u.s, I=u.g * u.cm**2, R=u.km)
+def pulsar_B(
+    f: u.Quantity,
+    fdot: u.Quantity,
+    I: u.Quantity = 1.0e45 * u.g * u.cm**2,
+    R: u.Quantity = 10 * u.km,
+) -> u.G:
     r"""Compute pulsar surface magnetic field
 
-    Return the estimated pulsar surface magnetic field strength
-    given the spin frequency and frequency derivative.
+    Return the pulsar surface magnetic field strength given the spin frequency `f` and frequency derivative `fdot`.
 
     Parameters
     ----------
@@ -232,6 +244,10 @@ def pulsar_B(f: u.Quantity, fdot: u.Quantity) -> u.G:
         pulsar frequency
     fdot : astropy.units.Quantity
         frequency derivative :math:`\dot f`
+    I : astropy.units.Quantity, optional
+        pulsar moment of inertia, default of 1e45 g*cm**2
+    R : astropy.units.Quantity, optional
+        pulsar radius, default of 10 km
 
     Returns
     -------
@@ -247,15 +263,19 @@ def pulsar_B(f: u.Quantity, fdot: u.Quantity) -> u.G:
 
     Notes
     -----
-    Calculates :math:`B=3.2\times 10^{19}\,{\rm  G}\sqrt{ f \dot f^{-3}}`
+    Calculates :math:`B=\sqrt{\frac{3\,I\,c^3}{8\pi^2\,R^6}\times\frac{-\dot{f}}{f^3}}`
     """
-    # This is a hack to use the traditional formula by stripping the units.
-    # It would be nice to improve this to a  proper formula with units
-    return 3.2e19 * u.G * np.sqrt(-fdot.to_value(u.Hz / u.s) / f.to_value(u.Hz) ** 3.0)
+    factor = (3.0 * I * const.c**3) / (8.0 * np.pi**2 * R**6)
+    return _to_gauss((factor * (-fdot) / f**3) ** 0.5)
 
 
-@u.quantity_input(f=u.Hz, fdot=u.Hz / u.s)
-def pulsar_B_lightcyl(f: u.Quantity, fdot: u.Quantity) -> u.G:
+@u.quantity_input(f=u.Hz, fdot=u.Hz / u.s, I=u.g * u.cm**2, R=u.km)
+def pulsar_B_lightcyl(
+    f: u.Quantity,
+    fdot: u.Quantity,
+    I: u.Quantity = 1.0e45 * u.g * u.cm**2,
+    R=10 * u.km,
+) -> u.G:
     r"""Compute pulsar magnetic field at the light cylinder
 
     Return the estimated pulsar magnetic field strength at the
@@ -268,6 +288,10 @@ def pulsar_B_lightcyl(f: u.Quantity, fdot: u.Quantity) -> u.G:
         pulsar frequency
     fdot : astropy.units.Quantity
         frequency derivative :math:`\dot f`
+    I : astropy.units.Quantity, optional
+        pulsar moment of inertia, default of 1e45 g*cm**2
+    R : astropy.units.Quantity, optional
+        pulsar radius, default of 10 km
 
     Returns
     -------
@@ -283,17 +307,10 @@ def pulsar_B_lightcyl(f: u.Quantity, fdot: u.Quantity) -> u.G:
 
     Notes
     -----
-    Calculates :math:`B_{LC} = 2.9\times 10^8\,{\rm G} P^{-5/2} \dot P^{1/2}`
+    Calculates :math:`B_{LC} = \sqrt{\frac{-24\pi^4\,I}{c^3}\dot{f}f^3}`
     """
-    p, pd = p_to_f(f, fdot)
-    # This is a hack to use the traditional formula by stripping the units.
-    # It would be nice to improve this to a  proper formula with units
-    return (
-        2.9e8
-        * u.G
-        * p.to_value(u.s) ** (-5.0 / 2.0)
-        * np.sqrt(pd.to(u.dimensionless_unscaled).value)
-    )
+    factor = 24.0 * np.pi**4.0 * I / const.c**3.0
+    return _to_gauss((factor * (-fdot) * f**3.0) ** 0.5)
 
 
 @u.quantity_input(pb=u.d, x=u.cm)
diff --git a/src/pint/models/noise_model.py b/src/pint/models/noise_model.py
index b6ae9a1b8..5bcd9eefb 100644
--- a/src/pint/models/noise_model.py
+++ b/src/pint/models/noise_model.py
@@ -126,6 +126,7 @@ def setup(self):
                             index=tneq_par.index,
                             aliases=["T2EQUAD"],
                             description="An error term added in quadrature to the scaled (by EFAC) TOA uncertainty.",
+                            convert_tcb2tdb=False,
                         )
                     )
                 EQUAD_par = getattr(self, EQUAD_name)
@@ -167,14 +168,14 @@ def scale_toa_sigma(self, toas, warn=True):
             if equad.quantity is None:
                 continue
             mask = equad.select_toa_mask(toas)
-            if np.any(mask):
+            if len(mask) > 0:
                 sigma_scaled[mask] = np.hypot(sigma_scaled[mask], equad.quantity)
             elif warn:
                 warnings.warn(f"EQUAD {equad} has no TOAs")
         for efac_name in self.EFACs:
             efac = getattr(self, efac_name)
             mask = efac.select_toa_mask(toas)
-            if np.any(mask):
+            if len(mask) > 0:
                 sigma_scaled[mask] *= efac.quantity
             elif warn:
                 warnings.warn(f"EFAC {efac} has no TOAs")
diff --git a/src/pint/models/parameter.py b/src/pint/models/parameter.py
index ea92f4611..5e3d9e607 100644
--- a/src/pint/models/parameter.py
+++ b/src/pint/models/parameter.py
@@ -709,7 +709,7 @@ def __init__(
 
         assert (
             not convert_tcb2tdb or tcb2tdb_scale_factor is not None
-        ), "Please specify the tcb2tdb_scale_factor explicitly."
+        ), f"Please specify the tcb2tdb_scale_factor explicitly for {name}."
         self.convert_tcb2tdb = convert_tcb2tdb
         self.tcb2tdb_scale_factor = tcb2tdb_scale_factor
 
@@ -1133,7 +1133,7 @@ def __init__(
 
         assert (
             not convert_tcb2tdb or tcb2tdb_scale_factor is not None
-        ), "Please specify the tcb2tdb_scale_factor explicitly."
+        ), f"Please specify the tcb2tdb_scale_factor explicitly for {name}."
         self.convert_tcb2tdb = convert_tcb2tdb
         self.tcb2tdb_scale_factor = tcb2tdb_scale_factor
 
@@ -1334,7 +1334,7 @@ def __init__(
 
         assert (
             not convert_tcb2tdb or tcb2tdb_scale_factor is not None
-        ), "Please specify the tcb2tdb_scale_factor explicitly."
+        ), f"Please specify the tcb2tdb_scale_factor explicitly for {name}."
         self.convert_tcb2tdb = convert_tcb2tdb
         self.tcb2tdb_scale_factor = tcb2tdb_scale_factor
 
@@ -1554,7 +1554,7 @@ def __init__(
         # a function of the prefix.
         assert (
             not convert_tcb2tdb or tcb2tdb_scale_factor is not None
-        ), "Please specify the tcb2tdb_scale_factor explicitly."
+        ), f"Please specify the tcb2tdb_scale_factor explicitly for {name}."
         tcb2tdb_scale_factor_val = (
             tcb2tdb_scale_factor(self.prefix)
             if hasattr(tcb2tdb_scale_factor, "__call__")
diff --git a/src/pint/models/timing_model.py b/src/pint/models/timing_model.py
index a0b27b090..68fe886bf 100644
--- a/src/pint/models/timing_model.py
+++ b/src/pint/models/timing_model.py
@@ -3425,6 +3425,7 @@ def add_param(self, param, deriv_func=None, setup=False):
         if deriv_func is not None:
             self.register_deriv_funcs(deriv_func, param.name)
         param._parent = self
+        return param.name
 
     def remove_param(self, param):
         """Remove a parameter from the Component.
diff --git a/src/pint/models/wave.py b/src/pint/models/wave.py
index 8a902a5ca..64e6da264 100644
--- a/src/pint/models/wave.py
+++ b/src/pint/models/wave.py
@@ -15,9 +15,8 @@ class Wave(PhaseComponent):
     sine/cosine components.
 
     For consistency with the implementation in tempo2, this signal is treated
-    as a time series, but trivially converted into phase by multiplication by
-    F0, which could makes changes to PEPOCH fragile if there is strong spin
-    frequency evolution.
+    as a time series and trivially converted into phase by multiplication by
+    F0.
 
     Parameters supported:
 
@@ -42,7 +41,7 @@ def __init__(self):
             floatParameter(
                 name="WAVE_OM",
                 description="Base frequency of wave solution",
-                units="1/d",
+                units="rad/d",
                 convert_tcb2tdb=False,
             )
         )
diff --git a/src/pint/pintk/pulsar.py b/src/pint/pintk/pulsar.py
index 65d8048df..ce6aa8560 100644
--- a/src/pint/pintk/pulsar.py
+++ b/src/pint/pintk/pulsar.py
@@ -151,7 +151,9 @@ def __contains__(self, key):
         return key in self.prefit_model.params
 
     def reset_model(self):
-        self.prefit_model = pint.models.get_model(self.parfile)
+        self.prefit_model = pint.models.get_model(
+            self.parfile, allow_T2=True, allow_tcb=True
+        )
         self.add_model_params()
         self.postfit_model = None
         self.postfit_resids = None
@@ -172,7 +174,9 @@ def reset_TOAs(self):
         self.update_resids()
 
     def resetAll(self):
-        self.prefit_model = pint.models.get_model(self.parfile)
+        self.prefit_model = pint.models.get_model(
+            self.parfile, allow_T2=True, allow_tcb=True
+        )
         self.postfit_model = None
         self.postfit_resids = None
         self.fitted = False
diff --git a/src/pint/plot_utils.py b/src/pint/plot_utils.py
index 45eba3977..289655f6d 100644
--- a/src/pint/plot_utils.py
+++ b/src/pint/plot_utils.py
@@ -283,7 +283,9 @@ def plot_priors(
         a, x = np.histogram(values[i], bins=bins, density=True)
         counts.append(a)
 
-    fig, axs = plt.subplots(len(keys), figsize=(8, 11), constrained_layout=True)
+    fig, axs = plt.subplots(
+        len(keys), figsize=(8, len(keys) * 1.5), constrained_layout=True
+    )
 
     for i, p in enumerate(keys):
         if i != len(keys[:-1]):
diff --git a/src/pint/scripts/compare_parfiles.py b/src/pint/scripts/compare_parfiles.py
index 0c4ca5da2..d65785041 100644
--- a/src/pint/scripts/compare_parfiles.py
+++ b/src/pint/scripts/compare_parfiles.py
@@ -83,14 +83,25 @@ def main(argv=None):
     parser.add_argument(
         "-q", "--quiet", default=0, action="count", help="Decrease output verbosity"
     )
+    parser.add_argument(
+        "--allow_tcb",
+        action="store_true",
+        help="Convert TCB par files to TDB automatically",
+    )
+    parser.add_argument(
+        "--allow_T2",
+        action="store_true",
+        help="Guess the underlying binary model when T2 is given",
+    )
 
     args = parser.parse_args(argv)
     pint.logging.setup(
         level=pint.logging.get_level(args.loglevel, args.verbosity, args.quiet)
     )
 
-    m1 = get_model(args.input1)
-    m2 = get_model(args.input2)
+    m1 = get_model(args.input1, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb)
+    m2 = get_model(args.input2, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb)
+
     print(
         m1.compare(
             m2,
diff --git a/src/pint/scripts/convert_parfile.py b/src/pint/scripts/convert_parfile.py
index 55b45415b..2975d5acc 100644
--- a/src/pint/scripts/convert_parfile.py
+++ b/src/pint/scripts/convert_parfile.py
@@ -73,6 +73,16 @@ def main(argv=None):
     parser.add_argument(
         "-q", "--quiet", default=0, action="count", help="Decrease output verbosity"
     )
+    parser.add_argument(
+        "--allow_tcb",
+        action="store_true",
+        help="Convert TCB par files to TDB automatically",
+    )
+    parser.add_argument(
+        "--allow_T2",
+        action="store_true",
+        help="Guess the underlying binary model when T2 is given",
+    )
 
     args = parser.parse_args(argv)
     pint.logging.setup(
@@ -83,7 +93,9 @@ def main(argv=None):
         return
 
     log.info(f"Reading '{args.input}'")
-    model = get_model(args.input)
+
+    model = get_model(args.input, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb)
+
     if hasattr(model, "BINARY") and args.binary is not None:
         log.info(f"Converting from {model.BINARY.value} to {args.binary}")
         if args.binary == "ELL1H":
diff --git a/src/pint/scripts/event_optimize.py b/src/pint/scripts/event_optimize.py
index e4c9dd4ff..6e26ba784 100755
--- a/src/pint/scripts/event_optimize.py
+++ b/src/pint/scripts/event_optimize.py
@@ -414,6 +414,16 @@ def __init__(
             self.model, phs, phserr
         )
         self.n_fit_params = len(self.fitvals)
+        self.M, _, _ = self.model.designmatrix(self.toas)
+        self.M = self.M.transpose() * -self.model.F0.value
+        self.phases = self.get_event_phases()
+        self.linearize_model = False
+
+    def calc_phase_matrix(self, theta):
+        d_phs = np.zeros(len(self.toas))
+        for i in range(len(theta) - 1):
+            d_phs += self.M[i + 1] * (self.fitvals[i] - theta[i])
+        return (self.phases - d_phs) % 1
 
     def get_event_phases(self):
         """
@@ -446,7 +456,10 @@ def lnposterior(self, theta):
             return -np.inf, -np.inf, -np.inf
 
         # Call PINT to compute the phases
-        phases = self.get_event_phases()
+        if self.linearize_model:
+            phases = self.calc_phase_matrix(theta)
+        else:
+            phases = self.get_event_phases()
         lnlikelihood = profile_likelihood(
             theta[-1], self.xtemp, phases, self.template, self.weights
         )
@@ -686,6 +699,23 @@ def main(argv=None):
         action="store_true",
         dest="noautocorr",
     )
+    parser.add_argument(
+        "--linearize_model",
+        help="Calculates the phase at each MCMC step using the designmatrix",
+        default=False,
+        action="store_true",
+        dest="linearize_model",
+    )
+    parser.add_argument(
+        "--allow_tcb",
+        action="store_true",
+        help="Convert TCB par files to TDB automatically",
+    )
+    parser.add_argument(
+        "--allow_T2",
+        action="store_true",
+        help="Guess the underlying binary model when T2 is given",
+    )
 
     args = parser.parse_args(argv)
     pint.logging.setup(
@@ -719,7 +749,9 @@ def main(argv=None):
     ncores = args.ncores
 
     # Read in initial model
-    modelin = pint.models.get_model(parfile)
+    modelin = pint.models.get_model(
+        parfile, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb
+    )
 
     # File name setup and clobber file check
     filepath = args.filepath or os.getcwd()
@@ -862,6 +894,9 @@ def main(argv=None):
     # This way, one walker should always be in a good position
     pos[0] = ftr.fitvals
 
+    # How phase will be calculated at each step (either with the designmatrix orthe exact phase calculation)
+    ftr.linearize_model = args.linearize_model
+
     import emcee
 
     # Setting up a backend to save the chains into an h5 file
@@ -925,7 +960,7 @@ def chains_to_dict(names, sampler):
 
     def plot_chains(chain_dict, file=False):
         npts = len(chain_dict)
-        fig, axes = plt.subplots(npts, 1, sharex=True, figsize=(8, 9))
+        fig, axes = plt.subplots(npts, 1, sharex=True, figsize=(8, npts * 1.5))
         for ii, name in enumerate(chain_dict.keys()):
             axes[ii].plot(chain_dict[name], color="k", alpha=0.3)
             axes[ii].set_ylabel(name)
@@ -950,6 +985,7 @@ def plot_chains(chain_dict, file=False):
     lnprior_samps = blobs["lnprior"]
     lnlikelihood_samps = blobs["lnlikelihood"]
     lnpost_samps = lnprior_samps + lnlikelihood_samps
+    maxpost = lnpost_samps[:][burnin:].max()
     ind = np.unravel_index(
         np.argmax(lnpost_samps[:][burnin:]), lnpost_samps[:][burnin:].shape
     )
@@ -1000,8 +1036,15 @@ def plot_chains(chain_dict, file=False):
     ]
     ftr.set_param_uncertainties(dict(zip(ftr.fitkeys[:-1], errors[:-1])))
 
+    # Calculating the AIC and BIC
+    n_params = len(ftr.model.free_params)
+    AIC = 2 * (n_params - maxpost)
+    BIC = n_params * np.log(len(ts)) - 2 * maxpost
+    ftr.model.NTOA.value = ts.ntoas
     f = open(filename + "_post.par", "w")
     f.write(ftr.model.as_parfile())
+    f.write(f"\n#The AIC is {AIC}")
+    f.write(f"\n#The BIC is {BIC}")
     f.close()
 
     # Print the best MCMC values and ranges
diff --git a/src/pint/scripts/event_optimize_MCMCFitter.py b/src/pint/scripts/event_optimize_MCMCFitter.py
index bbab8ccee..7357b1028 100755
--- a/src/pint/scripts/event_optimize_MCMCFitter.py
+++ b/src/pint/scripts/event_optimize_MCMCFitter.py
@@ -128,6 +128,17 @@ def main(argv=None):
         help="Logging level",
         dest="loglevel",
     )
+    parser.add_argument(
+        "--allow_tcb",
+        action="store_true",
+        help="Convert TCB par files to TDB automatically",
+    )
+    parser.add_argument(
+        "--allow_T2",
+        action="store_true",
+        help="Guess the underlying binary model when T2 is given",
+    )
+
     global nwalkers, nsteps, ftr
 
     args = parser.parse_args(argv)
@@ -164,7 +175,9 @@ def main(argv=None):
     wgtexp = args.wgtexp
 
     # Read in initial model
-    modelin = pint.models.get_model(parfile)
+    modelin = pint.models.get_model(
+        parfile, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb
+    )
 
     # The custom_timing version below is to manually construct the TimingModel
     # class, which allows it to be pickled. This is needed for parallelizing
diff --git a/src/pint/scripts/event_optimize_multiple.py b/src/pint/scripts/event_optimize_multiple.py
index 41316a000..0a3dc7336 100755
--- a/src/pint/scripts/event_optimize_multiple.py
+++ b/src/pint/scripts/event_optimize_multiple.py
@@ -228,6 +228,16 @@ def main(argv=None):
         help="Logging level",
         dest="loglevel",
     )
+    parser.add_argument(
+        "--allow_tcb",
+        action="store_true",
+        help="Convert TCB par files to TDB automatically",
+    )
+    parser.add_argument(
+        "--allow_T2",
+        action="store_true",
+        help="Guess the underlying binary model when T2 is given",
+    )
 
     global nwalkers, nsteps, ftr
 
@@ -261,7 +271,9 @@ def main(argv=None):
     wgtexp = args.wgtexp
 
     # Read in initial model
-    modelin = pint.models.get_model(parfile)
+    modelin = pint.models.get_model(
+        parfile, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb
+    )
 
     # Set the target coords for automatic weighting if necessary
     if "ELONG" in modelin.params:
diff --git a/src/pint/scripts/fermiphase.py b/src/pint/scripts/fermiphase.py
index b1427e61f..c5e43bb93 100755
--- a/src/pint/scripts/fermiphase.py
+++ b/src/pint/scripts/fermiphase.py
@@ -77,6 +77,16 @@ def main(argv=None):
     parser.add_argument(
         "-q", "--quiet", default=0, action="count", help="Decrease output verbosity"
     )
+    parser.add_argument(
+        "--allow_tcb",
+        action="store_true",
+        help="Convert TCB par files to TDB automatically",
+    )
+    parser.add_argument(
+        "--allow_T2",
+        action="store_true",
+        help="Guess the underlying binary model when T2 is given",
+    )
 
     args = parser.parse_args(argv)
     pint.logging.setup(
@@ -88,7 +98,10 @@ def main(argv=None):
         args.addphase = True
 
     # Read in model
-    modelin = pint.models.get_model(args.parfile)
+    modelin = pint.models.get_model(
+        args.parfile, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb
+    )
+
     if "ELONG" in modelin.params:
         tc = SkyCoord(
             modelin.ELONG.quantity,
diff --git a/src/pint/scripts/photonphase.py b/src/pint/scripts/photonphase.py
index 7d72eec18..2798e0a7f 100755
--- a/src/pint/scripts/photonphase.py
+++ b/src/pint/scripts/photonphase.py
@@ -106,12 +106,25 @@ def main(argv=None):
         help="Logging level",
         dest="loglevel",
     )
+    parser.add_argument(
+        "--nbin", help="Number of phase bins in the phaseogram", default=100, type=int
+    )
     parser.add_argument(
         "-v", "--verbosity", default=0, action="count", help="Increase output verbosity"
     )
     parser.add_argument(
         "-q", "--quiet", default=0, action="count", help="Decrease output verbosity"
     )
+    parser.add_argument(
+        "--allow_tcb",
+        action="store_true",
+        help="Convert TCB par files to TDB automatically",
+    )
+    parser.add_argument(
+        "--allow_T2",
+        action="store_true",
+        help="Guess the underlying binary model when T2 is given",
+    )
 
     args = parser.parse_args(argv)
     pint.logging.setup(
@@ -153,7 +166,10 @@ def main(argv=None):
                 "Please barycenter the event file using the official mission tools before processing with PINT"
             )
     # Read in model
-    modelin = pint.models.get_model(args.parfile)
+    modelin = pint.models.get_model(
+        args.parfile, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb
+    )
+
     use_planets = False
     if "PLANET_SHAPIRO" in modelin.params:
         if modelin.PLANET_SHAPIRO.value:
@@ -254,7 +270,7 @@ def main(argv=None):
         print("Htest : {0:.2f} ({1:.2f} sigma)".format(h, h2sig(h)))
 
     if args.plot:
-        phaseogram_binned(mjds, phases, bins=100, plotfile=args.plotfile)
+        phaseogram_binned(mjds, phases, bins=args.nbin, plotfile=args.plotfile)
 
     # Compute orbital phases for each photon TOA
     if args.addorbphase:
diff --git a/src/pint/scripts/pintbary.py b/src/pint/scripts/pintbary.py
index 4876474d8..61cdcd0a0 100755
--- a/src/pint/scripts/pintbary.py
+++ b/src/pint/scripts/pintbary.py
@@ -74,6 +74,16 @@ def main(argv=None):
     parser.add_argument(
         "-q", "--quiet", default=0, action="count", help="Decrease output verbosity"
     )
+    parser.add_argument(
+        "--allow_tcb",
+        action="store_true",
+        help="Convert TCB par files to TDB automatically",
+    )
+    parser.add_argument(
+        "--allow_T2",
+        action="store_true",
+        help="Guess the underlying binary model when T2 is given",
+    )
 
     args = parser.parse_args(argv)
     pint.logging.setup(
@@ -105,7 +115,9 @@ def main(argv=None):
     )
 
     if args.parfile is not None:
-        m = pint.models.get_model(args.parfile)
+        m = pint.models.get_model(
+            args.parfile, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb
+        )
     else:
         # Construct model by hand
         m = pint.models.StandardTimingModel
diff --git a/src/pint/scripts/pintempo.py b/src/pint/scripts/pintempo.py
index d99c772ba..ac9099008 100755
--- a/src/pint/scripts/pintempo.py
+++ b/src/pint/scripts/pintempo.py
@@ -62,6 +62,16 @@ def main(argv=None):
     parser.add_argument(
         "-q", "--quiet", default=0, action="count", help="Decrease output verbosity"
     )
+    parser.add_argument(
+        "--allow_tcb",
+        action="store_true",
+        help="Convert TCB par files to TDB automatically",
+    )
+    parser.add_argument(
+        "--allow_T2",
+        action="store_true",
+        help="Guess the underlying binary model when T2 is given",
+    )
 
     args = parser.parse_args(argv)
     pint.logging.setup(
@@ -69,7 +79,9 @@ def main(argv=None):
     )
 
     log.info("Reading model from {0}".format(args.parfile))
-    m = pint.models.get_model(args.parfile)
+    m = pint.models.get_model(
+        args.parfile, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb
+    )
 
     log.warning(m.params)
 
diff --git a/src/pint/scripts/pintpublish.py b/src/pint/scripts/pintpublish.py
index 6d5b0142f..1828a05c6 100644
--- a/src/pint/scripts/pintpublish.py
+++ b/src/pint/scripts/pintpublish.py
@@ -59,10 +59,22 @@ def main(argv=None):
         action="store_true",
         default=False,
     )
+    parser.add_argument(
+        "--allow_tcb",
+        action="store_true",
+        help="Convert TCB par files to TDB automatically",
+    )
+    parser.add_argument(
+        "--allow_T2",
+        action="store_true",
+        help="Guess the underlying binary model when T2 is given",
+    )
 
     args = parser.parse_args(argv)
 
-    model, toas = get_model_and_toas(args.parfile, args.timfile)
+    model, toas = get_model_and_toas(
+        args.parfile, args.timfile, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb
+    )
 
     output = publish(
         model,
diff --git a/src/pint/scripts/t2binary2pint.py b/src/pint/scripts/t2binary2pint.py
index 42d70f89a..bab43eef4 100644
--- a/src/pint/scripts/t2binary2pint.py
+++ b/src/pint/scripts/t2binary2pint.py
@@ -45,12 +45,17 @@ def main(argv=None):
         default=True,
         help="Whether to drop SINI if the model is DDK (True)",
     )
+    parser.add_argument(
+        "--allow_tcb",
+        action="store_true",
+        help="Convert TCB par files to TDB automatically",
+    )
 
     args = parser.parse_args(argv)
 
     mb = ModelBuilder()
 
-    model = mb(args.input_par, allow_T2=True, allow_tcb=True)
+    model = mb(args.input_par, allow_T2=True, allow_tcb=args.allow_tcb)
     model.write_parfile(args.output_par)
     print(f"Output written to {args.output_par}")
 
diff --git a/src/pint/scripts/tcb2tdb.py b/src/pint/scripts/tcb2tdb.py
index 4c427d21a..02a485d3c 100644
--- a/src/pint/scripts/tcb2tdb.py
+++ b/src/pint/scripts/tcb2tdb.py
@@ -30,11 +30,16 @@ def main(argv=None):
     )
     parser.add_argument("input_par", help="Input par file name (TCB)")
     parser.add_argument("output_par", help="Output par file name (TDB)")
+    parser.add_argument(
+        "--allow_T2",
+        action="store_true",
+        help="Guess the underlying binary model when T2 is given",
+    )
 
     args = parser.parse_args(argv)
 
     mb = ModelBuilder()
-    model = mb(args.input_par, allow_tcb=True)
+    model = mb(args.input_par, allow_tcb=True, allow_T2=args.allow_T2)
     model.write_parfile(args.output_par)
 
     log.info(f"Output written to {args.output_par}.")
diff --git a/src/pint/scripts/zima.py b/src/pint/scripts/zima.py
index 3305129f7..29b3e354c 100755
--- a/src/pint/scripts/zima.py
+++ b/src/pint/scripts/zima.py
@@ -115,6 +115,16 @@ def main(argv=None):
     parser.add_argument(
         "-q", "--quiet", default=0, action="count", help="Decrease output verbosity"
     )
+    parser.add_argument(
+        "--allow_tcb",
+        action="store_true",
+        help="Convert TCB par files to TDB automatically",
+    )
+    parser.add_argument(
+        "--allow_T2",
+        action="store_true",
+        help="Guess the underlying binary model when T2 is given",
+    )
 
     args = parser.parse_args(argv)
     pint.logging.setup(
@@ -122,7 +132,9 @@ def main(argv=None):
     )
 
     log.info("Reading model from {0}".format(args.parfile))
-    m = pint.models.get_model(args.parfile)
+    m = pint.models.get_model(
+        args.parfile, allow_T2=args.allow_T2, allow_tcb=args.allow_tcb
+    )
 
     out_format = args.format
     error = args.error * u.microsecond
diff --git a/src/pint/utils.py b/src/pint/utils.py
index da23092fc..59f0d3076 100644
--- a/src/pint/utils.py
+++ b/src/pint/utils.py
@@ -1745,11 +1745,8 @@ def _translate_wave_freqs(om: Union[float, u.Quantity], k: int) -> u.Quantity:
     astropy.units.Quantity
         WXFREQ_ quantity in units 1/d that can be used in WaveX model
     """
-    if isinstance(om, u.quantity.Quantity):
-        om.to(u.d**-1)
-    else:
-        om *= u.d**-1
-    return (om * (k + 1)) / (2.0 * np.pi)
+    om <<= u.rad / u.d
+    return (om * (k + 1)) / (2.0 * np.pi * u.rad)
 
 
 def _translate_wavex_freqs(wxfreq: Union[float, u.Quantity], k: int) -> u.Quantity:
@@ -1769,13 +1766,12 @@ def _translate_wavex_freqs(wxfreq: Union[float, u.Quantity], k: int) -> u.Quanti
     astropy.units.Quantity
         WAVEOM quantity in units 1/d that can be used in Wave model
     """
-    if isinstance(wxfreq, u.quantity.Quantity):
-        wxfreq.to(u.d**-1)
-    else:
-        wxfreq *= u.d**-1
+    wxfreq <<= u.d**-1
     if len(wxfreq) == 1:
-        return (2.0 * np.pi * wxfreq) / (k + 1.0)
-    wave_om = [((2.0 * np.pi * wxfreq[i]) / (k[i] + 1.0)) for i in range(len(wxfreq))]
+        return (2.0 * np.pi * u.rad * wxfreq) / (k + 1.0)
+    wave_om = [
+        ((2.0 * np.pi * u.rad * wxfreq[i]) / (k[i] + 1.0)) for i in range(len(wxfreq))
+    ]
     return (
         sum(wave_om) / len(wave_om)
         if np.allclose(wave_om, wave_om[0], atol=1e-3)
diff --git a/tests/datafile/vela_wave.par b/tests/datafile/vela_wave.par
index e8aba5a4c..889361b01 100644
--- a/tests/datafile/vela_wave.par
+++ b/tests/datafile/vela_wave.par
@@ -28,45 +28,18 @@ TZRMJD         55896.55312516020475
 TZRFRQ         1370.2919919999999365    
 TZRSITE        pks   
 TRES           310.453      
-EPHVER         2                           
 CLK            TT(TAI)
 MODE 1
 UNITS          TDB
 TIMEEPH        FB90
 DILATEFREQ     N
 PLANET_SHAPIRO N
-T2CMETHOD      TEMPO
 NE_SW          9.961
 CORRECT_TROPOSPHERE  N
 EPHEM          DE405
 NITS           1
 NTOA           339
 CHI2R          328088.1294 298
-JUMP -B 20CM 0 0
-JUMP -B 40CM 0 0
-JUMP -B 50CM 0 0
-JUMP -dfb3_J0437_55319_56160 1 2.2e-07 0
-JUMP -dfb3_J0437_56160_60000 1 4.5e-07 0
-JUMP -pdfb1_128_ch 1 -3.5547e-06 0
-JUMP -pdfb1_2048_ch 1 -4.68829e-05 0
-JUMP -pdfb1_512_ch 1 -1.22813e-05 0
-JUMP -pdfb1_post_2006 1 -1.3e-07 0
-JUMP -pdfb1_pre_2006 1 -1.13e-06 0
-JUMP -pdfb2_1024_MHz 1 -5.435e-06 0
-JUMP -pdfb2_256MHz_1024_ch 1 -1.1395e-05 0
-JUMP -pdfb2_256MHz_512ch 1 -4.75e-06 0
-JUMP -pdfb3_1024_256_512 1 2.45e-06 0
-JUMP -pdfb3_1024_MHz 1 1.03e-06 0
-JUMP -pdfb3_256MHz_1024ch 1 4.295e-06 0
-JUMP -pdfb3_64MHz_1024ch 1 1.494e-05 0
-JUMP -pdfb3_64MHz_512ch 1 8.9e-06 0
-JUMP -pdfb4.*_1024_[1,2]... 1 2.23e-06 0
-JUMP -pdfb4_256MHz_1024ch 1 5.05e-06 0
-JUMP -pdfb4_55319_56055_cals 1 9.27e-07 0
-JUMP -pdfb4_56110_56160_cals 1 5.41e-07 0
-JUMP -pdfb4_56160_60000_cals 1 4.25e-07 0
-JUMP -wbb256_512_128_3p_b 1 -6.2e-07 0
-JUMP -wbb_c_config 1 3.8e-07 0
 WAVEEPOCH 55305
 WAVE_OM 0.0015182579855022 0
 WAVE1 -0.21573979911255 -0.049063841960712
diff --git a/tests/test_derived_quantities.py b/tests/test_derived_quantities.py
index cf5e845b0..833e98a86 100644
--- a/tests/test_derived_quantities.py
+++ b/tests/test_derived_quantities.py
@@ -67,7 +67,7 @@ def test_Edot():
 def test_Bfield():
     # B
     assert np.isclose(
-        pulsar_B(0.033 * u.Hz, -2.0e-15 * u.Hz / u.s), 238722891596281.66 * u.G
+        pulsar_B(0.033 * u.Hz, -2.0e-15 * u.Hz / u.s), 238693670891966.22 * u.G
     )
 
 
@@ -75,7 +75,7 @@ def test_Blc():
     # B_lc
     assert np.isclose(
         pulsar_B_lightcyl(0.033 * u.Hz, -2.0e-15 * u.Hz / u.s),
-        0.07774704753236616 * u.G,
+        0.07896965114785195 * u.G,
     )
 
 
diff --git a/tests/test_model_wave.py b/tests/test_model_wave.py
index 1e57d640f..ab1707088 100644
--- a/tests/test_model_wave.py
+++ b/tests/test_model_wave.py
@@ -1,28 +1,81 @@
+from io import StringIO
 import os
 import pytest
 
 import astropy.units as u
+import numpy as np
 
-import pint.fitter
-import pint.models
+from pint.models import get_model
+from pint import toa
 import pint.residuals
-import pint.toa
 from pinttestdata import datadir
 
-# Not included in the test here, but as a sanity check I used this same
-# ephemeris to phase up Fermi data, and it looks good.
+par_nowave = """
+    PSRJ           J0835-4510
+    RAJ            08:35:20.61149
+    DECJ           -45:10:34.8751
+    F0             11.18965156782
+    PEPOCH         55305                       
+    DM             67.99
+    UNITS          TDB
+"""
 
-parfile = os.path.join(datadir, "vela_wave.par")
-timfile = os.path.join(datadir, "vela_wave.tim")
+wave_terms = """
+    WAVEEPOCH 55305
+    WAVE_OM 0.0015182579855022
+    WAVE1 -0.21573979911255 -0.049063841960712
+    WAVE2 0.62795320246729 -0.11984954655989
+    WAVE3 0.099618608456845 0.28530756908788
+    WAVE4 -0.21537340649058 0.18849486610628
+    WAVE5 0.021980474493165 -0.23566696662127
+"""
 
 
-class TestWave:
-    @classmethod
-    def setup_class(cls):
-        cls.m = pint.models.get_model(parfile)
-        cls.t = pint.toa.get_TOAs(timfile, ephem="DE405", include_bipm=False)
+def test_wave_ephem():
+    parfile = os.path.join(datadir, "vela_wave.par")
+    timfile = os.path.join(datadir, "vela_wave.tim")
+    m = get_model(parfile)
+    t = toa.get_TOAs(timfile, ephem="DE405", include_bipm=False)
+    rs = pint.residuals.Residuals(t, m).time_resids
+    assert rs.std() < 350.0 * u.us
 
-    def test_vela_rms_is_small_enough(self):
-        rs = pint.residuals.Residuals(self.t, self.m).time_resids
-        rms = rs.to(u.us).std()
-        assert rms < 350.0 * u.us
+
+def test_wave_construction():
+    m = get_model(StringIO(par_nowave + wave_terms))
+    assert np.allclose(m.WAVE_OM.quantity, 0.0015182579855022 * u.rad / u.day)
+    assert np.allclose(m.WAVE1.quantity[0], -0.21573979911255 * u.s)
+    assert np.allclose(m.WAVE2.quantity[1], -0.1198495465598 * u.s)
+
+
+def test_wave_computation():
+    m0 = get_model(StringIO(par_nowave))
+    m1 = get_model(StringIO(par_nowave + wave_terms))
+    # make some barycentric TOAs
+    tdbs = np.linspace(54500, 60000, 10)
+    ts = toa.TOAs(
+        toalist=[
+            toa.TOA(t, obs="@", freq=np.inf, error=1 * u.ms, scale="tdb") for t in tdbs
+        ]
+    )
+    ts.compute_TDBs(ephem="DE421")
+    ts.compute_posvels(ephem="DE421")
+    ph0 = m0.phase(ts)
+    ph1 = m1.phase(ts)
+    dphi = (ph1.int - ph0.int) + (ph1.frac - ph0.frac)
+    test_phase = np.zeros(len(tdbs))
+    WAVEEPOCH = 55305
+    WAVE_OM = 0.0015182579855022
+    WAVE = [
+        [-0.21573979911255, -0.049063841960712],
+        [0.62795320246729, -0.11984954655989],
+        [0.099618608456845, 0.28530756908788],
+        [-0.21537340649058, 0.18849486610628],
+        [0.021980474493165, -0.23566696662127],
+    ]
+    ph = (tdbs - WAVEEPOCH) * WAVE_OM
+    for i in range(5):
+        test_phase += WAVE[i][0] * np.sin((i + 1) * ph) + WAVE[i][1] * np.cos(
+            (i + 1) * ph
+        )
+    test_phase *= m0.F0.quantity.to(u.Hz).value
+    assert np.allclose(test_phase, dphi)
diff --git a/tox.ini b/tox.ini
index 4b7f6b153..efadd9b9c 100644
--- a/tox.ini
+++ b/tox.ini
@@ -14,7 +14,7 @@ envlist =
     codestyle
     black
     singletest
-    py{38,39,310,311,312}-test{,-alldeps,-devdeps}{,-cov}
+    py{38,39,310,311,312,313}-test{,-alldeps,-devdeps}{,-cov}
 
 skip_missing_interpreters = True
 
@@ -48,8 +48,8 @@ commands =
     cov: coverage xml -o {toxinidir}/coverage.xml
 
 depends =
-    {py38,py39,py310,py311,py312}: clean
-    report: py38,py39,py310,py312
+    {py39,py310,py311,py312,313}: clean
+    report: py39,py310,py312,py313
     docs: notebooks
 
 [testenv:singletest]
@@ -80,13 +80,13 @@ deps =
 [testenv:oldestdeps]
 description =
     Run tests on Python 3 with minimum supported versions of astropy, numpy
-basepython = python3.8
+basepython = python3.9
 deps =
-    numpy==1.18.5
+    numpy==1.23.0
     numdifftools==0.9.39
-    astropy==4.0
-    matplotlib==3.2.0
-    scipy==1.4.1
+    astropy==5.0.5
+    matplotlib==3.4.3
+    scipy==1.9.0
     pytest
     coverage
     hypothesis<=6.72.0
@@ -102,7 +102,7 @@ commands =
 
 [testenv:notebooks]
 description = update the notebooks
-basepython = python3.12
+basepython = python3.13
 deps =
     traitlets
     sphinx >= 2.2, != 5.1.0
@@ -130,7 +130,7 @@ commands = coverage erase
 [testenv:docs]
 changedir = {toxinidir}/docs
 description = invoke sphinx-build to build the HTML docs
-basepython = python3.12
+basepython = python3.13
 deps =
     traitlets
     sphinx >= 2.2, != 5.1.0