From 2d72c75fafb00a449945fa3b0df91d7ad5105e6d Mon Sep 17 00:00:00 2001
From: David Kaplan <kaplan@uwm.edu>
Date: Tue, 9 Jul 2024 09:46:21 -0500
Subject: [PATCH 01/11] type hints for derived quantities

---
 CHANGELOG-unreleased.md        |   1 +
 src/pint/derived_quantities.py | 129 +++++++++++++++++++++------------
 2 files changed, 82 insertions(+), 48 deletions(-)

diff --git a/CHANGELOG-unreleased.md b/CHANGELOG-unreleased.md
index d8ff0a4b7..77def071d 100644
--- a/CHANGELOG-unreleased.md
+++ b/CHANGELOG-unreleased.md
@@ -39,6 +39,7 @@ the released changes.
 - `pint.models.chromatic_model.ChromaticCM` for a Taylor series representation of the variable-index chromatic delay.
 - Whitened residuals (`white-res`) as a plotting axis in `pintk`
 - `TOAs.get_Tspan()` method
+- Type hints in `pint.derived_quantities`
 ### Fixed
 - `pint.utils.split_swx()` to use updated `SolarWindDispersionX()` parameter naming convention 
 - Fix #1759 by changing order of comparison
diff --git a/src/pint/derived_quantities.py b/src/pint/derived_quantities.py
index 104b73d8d..a69331f13 100644
--- a/src/pint/derived_quantities.py
+++ b/src/pint/derived_quantities.py
@@ -4,6 +4,7 @@
 import astropy.constants as const
 import astropy.units as u
 import numpy as np
+from typing import Optional, List, Tuple, Union
 
 import pint
 
@@ -33,7 +34,9 @@
 @u.quantity_input(
     p=[u.Hz, u.s], pd=[u.Hz / u.s, u.s / u.s], pdd=[u.Hz / u.s**2, u.s / u.s**2]
 )
-def p_to_f(p, pd, pdd=None):
+def p_to_f(
+    p: u.Quantity, pd: u.Quantity, pdd: Optional[u.Quantity] = None
+) -> Tuple[u.Quantity]:
     """Converts P, Pdot to F, Fdot (or vice versa)
 
     Convert period, period derivative and period second
@@ -71,7 +74,11 @@ def p_to_f(p, pd, pdd=None):
         if pdd == 0.0
         else 2.0 * pd * pd / (p**3.0) - pdd / (p * p)
     )
-    return [f, fd, fdd]
+    return (f, fd, fdd)
+
+
+# alias for the above
+f_to_p = p_to_f
 
 
 @u.quantity_input(
@@ -80,7 +87,12 @@ def p_to_f(p, pd, pdd=None):
     pdorfd=[u.Hz / u.s, u.s / u.s],
     pdorfderr=[u.Hz / u.s, u.s / u.s],
 )
-def pferrs(porf, porferr, pdorfd=None, pdorfderr=None):
+def pferrs(
+    porf: u.Quantity,
+    porferr: u.Quantity,
+    pdorfd: Optional[u.Quantity] = None,
+    pdorfderr: Optional[u.Quantity] = None,
+) -> Tuple[u.Quantity]:
     """Convert P, Pdot to F, Fdot with uncertainties (or vice versa).
 
     Calculate the period or frequency errors and
@@ -120,15 +132,16 @@ def pferrs(porf, porferr, pdorfd=None, pdorfderr=None):
         return [1.0 / porf, porferr / porf**2.0]
     forperr = porferr / porf**2.0
     fdorpderr = np.sqrt(
-        (4.0 * pdorfd**2.0 * porferr**2.0) / porf**6.0
-        + pdorfderr**2.0 / porf**4.0
+        (4.0 * pdorfd**2.0 * porferr**2.0) / porf**6.0 + pdorfderr**2.0 / porf**4.0
     )
     [forp, fdorpd] = p_to_f(porf, pdorfd)
-    return [forp, forperr, fdorpd, fdorpderr]
+    return (forp, forperr, fdorpd, fdorpderr)
 
 
-@u.quantity_input(fo=u.Hz)
-def pulsar_age(f: u.Hz, fdot: u.Hz / u.s, n=3, fo=1e99 * u.Hz):
+@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
+) -> u.Quantity:
     """Compute pulsar characteristic age
 
     Return the age of a pulsar given the spin frequency
@@ -170,8 +183,10 @@ def pulsar_age(f: u.Hz, fdot: u.Hz / u.s, n=3, fo=1e99 * u.Hz):
     return (-f / ((n - 1.0) * fdot) * (1.0 - (f / fo) ** (n - 1.0))).to(u.yr)
 
 
-@u.quantity_input(I=u.g * u.cm**2)
-def pulsar_edot(f: u.Hz, fdot: u.Hz / u.s, I=1.0e45 * u.g * u.cm**2):
+@u.quantity_input(f=u.Hz, fdot=u.Hz / u.s, I=u.g * u.cm**2)
+def pulsar_edot(
+    f: u.Quantity, fdot: u.Quantity, I: u.Quantity = 1.0e45 * u.g * u.cm**2
+) -> u.Quantity:
     """Compute pulsar spindown energy loss rate
 
     Return the pulsar `Edot` (:math:`\dot E`, in erg/s) given the spin frequency `f` and
@@ -206,8 +221,8 @@ def pulsar_edot(f: u.Hz, fdot: u.Hz / u.s, I=1.0e45 * u.g * u.cm**2):
     return (-4.0 * np.pi**2 * I * f * fdot).to(u.erg / u.s)
 
 
-@u.quantity_input
-def pulsar_B(f: u.Hz, fdot: u.Hz / u.s):
+@u.quantity_input(f=u.Hz, fdot=u.Hz / u.s)
+def pulsar_B(f: u.Quantity, fdot: u.Quantity) -> u.Quantity:
     """Compute pulsar surface magnetic field
 
     Return the estimated pulsar surface magnetic field strength
@@ -241,8 +256,8 @@ def pulsar_B(f: u.Hz, fdot: u.Hz / u.s):
     return 3.2e19 * u.G * np.sqrt(-fdot.to_value(u.Hz / u.s) / f.to_value(u.Hz) ** 3.0)
 
 
-@u.quantity_input
-def pulsar_B_lightcyl(f: u.Hz, fdot: u.Hz / u.s):
+@u.quantity_input(f=u.Hz, fdot=u.Hz / u.s)
+def pulsar_B_lightcyl(f: u.Quantity, fdot: u.Quantity) -> u.Quantity:
     """Compute pulsar magnetic field at the light cylinder
 
     Return the estimated pulsar magnetic field strength at the
@@ -283,8 +298,8 @@ def pulsar_B_lightcyl(f: u.Hz, fdot: u.Hz / u.s):
     )
 
 
-@u.quantity_input
-def mass_funct(pb: u.d, x: u.cm):
+@u.quantity_input(pb=u.d, x=u.cm)
+def mass_funct(pb: u.Quantity, x: u.Quantity) -> u.Quantity:
     """Compute binary mass function from period and semi-major axis
 
     Can handle scalar or array inputs.
@@ -324,8 +339,8 @@ def mass_funct(pb: u.d, x: u.cm):
     return fm.to(u.solMass)
 
 
-@u.quantity_input
-def mass_funct2(mp: u.Msun, mc: u.Msun, i: u.deg):
+@u.quantity_input(mp=u.Msun, mc=u.Msun, i=u.deg)
+def mass_funct2(mp: u.Quantity, mc: u.Quantity, i: u.Quantity) -> u.Quantity:
     """Compute binary mass function from masses and inclination
 
     Can handle scalar or array inputs.
@@ -369,8 +384,10 @@ def mass_funct2(mp: u.Msun, mc: u.Msun, i: u.deg):
     return (mc * np.sin(i)) ** 3.0 / (mc + mp) ** 2.0
 
 
-@u.quantity_input
-def pulsar_mass(pb: u.d, x: u.cm, mc: u.Msun, i: u.deg):
+@u.quantity_input(pb=u.d, x=u.cm, mc=u.Msun, i=u.deg)
+def pulsar_mass(
+    pb: u.Quantity, x: u.Quantity, mc: u.Quantity, i: u.Quantity
+) -> u.Quantity:
     """Compute pulsar mass from orbital parameters
 
     Return the pulsar mass (in solar mass units) for a binary.
@@ -436,8 +453,13 @@ def pulsar_mass(pb: u.d, x: u.cm, mc: u.Msun, i: u.deg):
     return ((-cb + np.sqrt(4 * massfunct * mc**3 * sini**3)) / (2 * ca)).to(u.Msun)
 
 
-@u.quantity_input(inc=u.deg, mpsr=u.solMass)
-def companion_mass(pb: u.d, x: u.cm, i=60.0 * u.deg, mp=1.4 * u.solMass):
+@u.quantity_input(pb=u.d, x=u.cm, i=u.deg, mp=u.solMass)
+def companion_mass(
+    pb: u.Quantity,
+    x: u.Quantity,
+    i: u.Quantity = 60.0 * u.deg,
+    mp: u.Quantity = 1.4 * u.solMass,
+) -> u.Quantity:
     """Commpute the companion mass from the orbital parameters
 
     Compute companion mass for a binary system from orbital mechanics,
@@ -515,17 +537,12 @@ def companion_mass(pb: u.d, x: u.cm, i=60.0 * u.deg, mp=1.4 * u.solMass):
     # delta1 is always <0
     # delta1 = 2 * b ** 3 - 9 * a * b * c + 27 * a ** 2 * d
     delta1 = (
-        -2 * massfunct**3
-        - 18 * a * mp * massfunct**2
-        - 27 * a**2 * massfunct * mp**2
+        -2 * massfunct**3 - 18 * a * mp * massfunct**2 - 27 * a**2 * massfunct * mp**2
     )
     # Q**2 is always > 0, so this is never a problem
     # in terms of complex numbers
     # Q = np.sqrt(delta1**2 - 4*delta0**3)
-    Q = np.sqrt(
-        108 * a**3 * mp**3 * massfunct**3
-        + 729 * a**4 * mp**4 * massfunct**2
-    )
+    Q = np.sqrt(108 * a**3 * mp**3 * massfunct**3 + 729 * a**4 * mp**4 * massfunct**2)
     # this could be + or - Q
     # pick the - branch since delta1 is <0 so that delta1 - Q is never near 0
     Ccubed = 0.5 * (delta1 + Q)
@@ -540,8 +557,10 @@ def companion_mass(pb: u.d, x: u.cm, i=60.0 * u.deg, mp=1.4 * u.solMass):
     return x1.to(u.Msun)
 
 
-@u.quantity_input
-def pbdot(mp: u.Msun, mc: u.Msun, pb: u.d, e: u.dimensionless_unscaled):
+@u.quantity_input(mp=u.Msun, mc=u.Msun, pb=u.d, e=u.dimensionless_unscaled)
+def pbdot(
+    mp: u.Quantity, mc: u.Quantity, pb: u.Quantity, e: Union[float, u.Quantity]
+) -> u.Quantity:
     """Post-Keplerian orbital decay pbdot, assuming general relativity.
 
     pbdot (:math:`\dot P_B`) is the change in the binary orbital period
@@ -603,8 +622,13 @@ def pbdot(mp: u.Msun, mc: u.Msun, pb: u.d, e: u.dimensionless_unscaled):
     return value.to(u.s / u.s)
 
 
-@u.quantity_input
-def gamma(mp: u.Msun, mc: u.Msun, pb: u.d, e: u.dimensionless_unscaled):
+@u.quantity_input(mp=u.Msun, mc=u.Msun, pb=u.d, e=u.dimensionless_unscaled)
+def gamma(
+    mp: u.Quantity,
+    mc: u.Quantity,
+    pb: u.Quantity,
+    e: Union[float, u.Quantity],
+) -> u.Quantity:
     """Post-Keplerian time dilation and gravitational redshift gamma, assuming general relativity.
 
     gamma (:math:`\gamma`) is the amplitude of the modification in arrival times caused by the varying
@@ -659,8 +683,13 @@ def gamma(mp: u.Msun, mc: u.Msun, pb: u.d, e: u.dimensionless_unscaled):
     return value.to(u.s)
 
 
-@u.quantity_input
-def omdot(mp: u.Msun, mc: u.Msun, pb: u.d, e: u.dimensionless_unscaled):
+@u.quantity_input(mp=u.Msun, mc=u.Msun, pb=u.d, e=u.dimensionless_unscaled)
+def omdot(
+    mp: u.Quantity,
+    mc: u.Quantity,
+    pb: u.Quantity,
+    e: Union[float, u.Quantity],
+) -> u.Quantity:
     """Post-Keplerian longitude of periastron precession rate omdot, assuming general relativity.
 
     omdot (:math:`\dot \omega`) is the relativistic advance of periastron.
@@ -714,8 +743,8 @@ def omdot(mp: u.Msun, mc: u.Msun, pb: u.d, e: u.dimensionless_unscaled):
     return value.to(u.deg / u.yr, equivalencies=u.dimensionless_angles())
 
 
-@u.quantity_input
-def sini(mp: u.Msun, mc: u.Msun, pb: u.d, x: u.cm):
+@u.quantity_input(mp=u.Msun, mc=u.Msun, pb=u.d, x=u.cm)
+def sini(mp: u.Quantity, mc: u.Quantity, pb: u.Quantity, x: u.Quantity) -> u.Quantity:
     """Post-Keplerian sine of inclination, assuming general relativity.
 
     Can handle scalar or array inputs.
@@ -768,8 +797,8 @@ def sini(mp: u.Msun, mc: u.Msun, pb: u.d, x: u.cm):
     ).decompose()
 
 
-@u.quantity_input
-def dr(mp: u.Msun, mc: u.Msun, pb: u.d):
+@u.quantity_input(mp=u.Msun, mc=u.Msun, pb=u.d)
+def dr(mp: u.Quantity, mc: u.Quantity, pb: u.Quantity) -> u.Quantity:
     """Post-Keplerian Roemer delay term
 
     dr (:math:`\delta_r`) is part of the relativistic deformation of the orbit
@@ -818,8 +847,8 @@ def dr(mp: u.Msun, mc: u.Msun, pb: u.d):
     ).decompose()
 
 
-@u.quantity_input
-def dth(mp: u.Msun, mc: u.Msun, pb: u.d):
+@u.quantity_input(mp=u.Msun, mc=u.Msun, pb=u.d)
+def dth(mp: u.Quantity, mc: u.Quantity, pb: u.Quantity) -> u.Quantity:
     """Post-Keplerian Roemer delay term
 
     dth (:math:`\delta_{\\theta}`) is part of the relativistic deformation of the orbit
@@ -868,8 +897,10 @@ def dth(mp: u.Msun, mc: u.Msun, pb: u.d):
     ).decompose()
 
 
-@u.quantity_input
-def omdot_to_mtot(omdot: u.deg / u.yr, pb: u.d, e: u.dimensionless_unscaled):
+@u.quantity_input(omdot=u.deg / u.yr, pb=u.d, e=u.dimensionless_unscaled)
+def omdot_to_mtot(
+    omdot: u.Quantity, pb: u.Quantity, e: Union[float, u.Quantity]
+) -> u.Quantity:
     """Determine total mass from Post-Keplerian longitude of periastron precession rate omdot,
     assuming general relativity.
 
@@ -931,7 +962,9 @@ def omdot_to_mtot(omdot: u.deg / u.yr, pb: u.d, e: u.dimensionless_unscaled):
 
 
 @u.quantity_input(pb=u.d, mp=u.Msun, mc=u.Msun, i=u.deg)
-def a1sini(mp, mc, pb, i=90 * u.deg):
+def a1sini(
+    mp: u.Quantity, mc: u.Quantity, pb: u.Quantity, i: u.Quantity = 90 * u.deg
+) -> u.Quantity:
     """Pulsar's semi-major axis.
 
     The full semi-major axis is given by Kepler's third law.  This is the
@@ -982,8 +1015,8 @@ def a1sini(mp, mc, pb, i=90 * u.deg):
     ).to(pint.ls)
 
 
-@u.quantity_input
-def shklovskii_factor(pmtot: u.mas / u.yr, D: u.kpc):
+@u.quantity_input(pmtot=u.mas / u.yr, D=u.kpc)
+def shklovskii_factor(pmtot: u.Quantity, D: u.Quantity) -> u.Quantity:
     """Compute magnitude of Shklovskii correction factor.
 
     Computes the Shklovskii correction factor, as defined in Eq 8.12 of Lorimer & Kramer (2005) [10]_
@@ -1020,8 +1053,8 @@ def shklovskii_factor(pmtot: u.mas / u.yr, D: u.kpc):
     return a_s
 
 
-@u.quantity_input
-def dispersion_slope(dm: pint.dmu):
+@u.quantity_input(dm=pint.dmu)
+def dispersion_slope(dm: u.Quantity) -> u.Quantity:
     """Compute the dispersion slope.
 
     This is equal to DMconst * DM.

From 6b14ddc84145ec3ab2f9364fe3d4956d3e3b889c Mon Sep 17 00:00:00 2001
From: David Kaplan <kaplan@uwm.edu>
Date: Tue, 9 Jul 2024 10:13:52 -0500
Subject: [PATCH 02/11] typo

---
 src/pint/derived_quantities.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pint/derived_quantities.py b/src/pint/derived_quantities.py
index a69331f13..d0f06aae0 100644
--- a/src/pint/derived_quantities.py
+++ b/src/pint/derived_quantities.py
@@ -140,7 +140,7 @@ def pferrs(
 
 @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
+    f: u.Quantity, fdot: u.Quantity, n: int = 3, fo: u.Quantity = 1e99 * u.Hz
 ) -> u.Quantity:
     """Compute pulsar characteristic age
 

From 5c4ce6a72442f272e37928d8d0cea92be1a4bf99 Mon Sep 17 00:00:00 2001
From: David Kaplan <kaplan@uwm.edu>
Date: Tue, 9 Jul 2024 10:22:26 -0500
Subject: [PATCH 03/11] trying to avoid black errors

---
 src/pint/derived_quantities.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pint/derived_quantities.py b/src/pint/derived_quantities.py
index d0f06aae0..41a5e53c3 100644
--- a/src/pint/derived_quantities.py
+++ b/src/pint/derived_quantities.py
@@ -631,7 +631,7 @@ def gamma(
 ) -> u.Quantity:
     """Post-Keplerian time dilation and gravitational redshift gamma, assuming general relativity.
 
-    gamma (:math:`\gamma`) is the amplitude of the modification in arrival times caused by the varying
+    gamma (:math:`\\gamma`) is the amplitude of the modification in arrival times caused by the varying
     gravitational redshift of the companion and time dilation in an elliptical orbit.  The time delay is
     :math:`\gamma \sin E`, where :math:`E` is the eccentric anomaly.
     Can handle scalar or array inputs.

From 14f3b9eb466f5d11b100c18186cf25f179bdc4f5 Mon Sep 17 00:00:00 2001
From: David Kaplan <kaplan@uwm.edu>
Date: Tue, 9 Jul 2024 10:23:32 -0500
Subject: [PATCH 04/11] trying to avoid black errors

---
 src/pint/derived_quantities.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/pint/derived_quantities.py b/src/pint/derived_quantities.py
index 41a5e53c3..8a9e0a426 100644
--- a/src/pint/derived_quantities.py
+++ b/src/pint/derived_quantities.py
@@ -147,7 +147,7 @@ def pulsar_age(
     Return the age of a pulsar given the spin frequency
     and frequency derivative.  By default, the characteristic age
     is returned (assuming a braking index `n` =3 and an initial
-    spin frequency :math:`f_0 \gg f`).  But `n` and `fo` can be set.
+    spin frequency :math:`f_0 \\gg f`).  But `n` and `fo` can be set.
 
     Parameters
     ----------
@@ -633,7 +633,7 @@ def gamma(
 
     gamma (:math:`\\gamma`) is the amplitude of the modification in arrival times caused by the varying
     gravitational redshift of the companion and time dilation in an elliptical orbit.  The time delay is
-    :math:`\gamma \sin E`, where :math:`E` is the eccentric anomaly.
+    :math:`\\gamma \sin E`, where :math:`E` is the eccentric anomaly.
     Can handle scalar or array inputs.
 
     Parameters
@@ -664,7 +664,7 @@ def gamma(
     Calculates
 
     .. math::
-        \gamma = T_{\odot}^{2/3} \\left(\\frac{P_b}{2\pi}\\right)^{1/3} e \\frac{m_c(m_p+2m_c)}{(m_p+m_c)^{4/3}}
+        \\gamma = T_{\odot}^{2/3} \\left(\\frac{P_b}{2\pi}\\right)^{1/3} e \\frac{m_c(m_p+2m_c)}{(m_p+m_c)^{4/3}}
 
     with :math:`T_\odot = GM_\odot c^{-3}`.
 

From d9039a47146d531aa2387b8d1002a0b31d78bc7b Mon Sep 17 00:00:00 2001
From: David Kaplan <kaplan@uwm.edu>
Date: Tue, 9 Jul 2024 10:29:34 -0500
Subject: [PATCH 05/11] trying to avoid black errors

---
 src/pint/derived_quantities.py | 94 +++++++++++++++++-----------------
 1 file changed, 47 insertions(+), 47 deletions(-)

diff --git a/src/pint/derived_quantities.py b/src/pint/derived_quantities.py
index 8a9e0a426..5c82ec9aa 100644
--- a/src/pint/derived_quantities.py
+++ b/src/pint/derived_quantities.py
@@ -93,7 +93,7 @@ def pferrs(
     pdorfd: Optional[u.Quantity] = None,
     pdorfderr: Optional[u.Quantity] = None,
 ) -> Tuple[u.Quantity]:
-    """Convert P, Pdot to F, Fdot with uncertainties (or vice versa).
+    r"""Convert P, Pdot to F, Fdot with uncertainties (or vice versa).
 
     Calculate the period or frequency errors and
     the Pdot or fdot errors from the opposite ones.
@@ -142,12 +142,12 @@ def pferrs(
 def pulsar_age(
     f: u.Quantity, fdot: u.Quantity, n: int = 3, fo: u.Quantity = 1e99 * u.Hz
 ) -> u.Quantity:
-    """Compute pulsar characteristic age
+    r"""Compute pulsar characteristic age
 
     Return the age of a pulsar given the spin frequency
     and frequency derivative.  By default, the characteristic age
     is returned (assuming a braking index `n` =3 and an initial
-    spin frequency :math:`f_0 \\gg f`).  But `n` and `fo` can be set.
+    spin frequency :math:`f_0 \gg f`).  But `n` and `fo` can be set.
 
     Parameters
     ----------
@@ -178,7 +178,7 @@ def pulsar_age(
 
     .. math::
 
-        \\tau = \\frac{f}{(n-1)\dot f}\\left(1-\\left(\\frac{f}{f_0}\\right)^{n-1}\\right)
+        \tau = \frac{f}{(n-1)\dot f}\left(1-\left(\frac{f}{f_0}\right)^{n-1}\right)
     """
     return (-f / ((n - 1.0) * fdot) * (1.0 - (f / fo) ** (n - 1.0))).to(u.yr)
 
@@ -187,7 +187,7 @@ def pulsar_age(
 def pulsar_edot(
     f: u.Quantity, fdot: u.Quantity, I: u.Quantity = 1.0e45 * u.g * u.cm**2
 ) -> u.Quantity:
-    """Compute pulsar spindown energy loss rate
+    r"""Compute pulsar spindown energy loss rate
 
     Return the pulsar `Edot` (:math:`\dot E`, in erg/s) given the spin frequency `f` and
     frequency derivative `fdot`. The NS moment of inertia is assumed to be
@@ -223,7 +223,7 @@ def pulsar_edot(
 
 @u.quantity_input(f=u.Hz, fdot=u.Hz / u.s)
 def pulsar_B(f: u.Quantity, fdot: u.Quantity) -> u.Quantity:
-    """Compute pulsar surface magnetic field
+    r"""Compute pulsar surface magnetic field
 
     Return the estimated pulsar surface magnetic field strength
     given the spin frequency and frequency derivative.
@@ -249,7 +249,7 @@ def pulsar_B(f: u.Quantity, fdot: u.Quantity) -> u.Quantity:
 
     Notes
     -----
-    Calculates :math:`B=3.2\\times 10^{19}\\,{\\rm  G}\\sqrt{ f \dot f^{-3}}`
+    Calculates :math:`B=3.2\times 10^{19}\,{\rm  G}\sqrt{ f \dot 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
@@ -258,7 +258,7 @@ def pulsar_B(f: u.Quantity, fdot: u.Quantity) -> u.Quantity:
 
 @u.quantity_input(f=u.Hz, fdot=u.Hz / u.s)
 def pulsar_B_lightcyl(f: u.Quantity, fdot: u.Quantity) -> u.Quantity:
-    """Compute pulsar magnetic field at the light cylinder
+    r"""Compute pulsar magnetic field at the light cylinder
 
     Return the estimated pulsar magnetic field strength at the
     light cylinder given the spin frequency and
@@ -285,7 +285,7 @@ def pulsar_B_lightcyl(f: u.Quantity, fdot: u.Quantity) -> u.Quantity:
 
     Notes
     -----
-    Calculates :math:`B_{LC} = 2.9\\times 10^8\\,{\\rm G} P^{-5/2} \dot P^{1/2}`
+    Calculates :math:`B_{LC} = 2.9\times 10^8\,{\rm G} P^{-5/2} \dot P^{1/2}`
     """
     p, pd = p_to_f(f, fdot)
     # This is a hack to use the traditional formula by stripping the units.
@@ -300,7 +300,7 @@ def pulsar_B_lightcyl(f: u.Quantity, fdot: u.Quantity) -> u.Quantity:
 
 @u.quantity_input(pb=u.d, x=u.cm)
 def mass_funct(pb: u.Quantity, x: u.Quantity) -> u.Quantity:
-    """Compute binary mass function from period and semi-major axis
+    r"""Compute binary mass function from period and semi-major axis
 
     Can handle scalar or array inputs.
 
@@ -329,7 +329,7 @@ def mass_funct(pb: u.Quantity, x: u.Quantity) -> u.Quantity:
 
     .. math::
 
-        f(m_p, m_c) = \\frac{4\pi^2 x^3}{G P_b^2}
+        f(m_p, m_c) = \frac{4\pi^2 x^3}{G P_b^2}
 
     See [1]_
 
@@ -341,7 +341,7 @@ def mass_funct(pb: u.Quantity, x: u.Quantity) -> u.Quantity:
 
 @u.quantity_input(mp=u.Msun, mc=u.Msun, i=u.deg)
 def mass_funct2(mp: u.Quantity, mc: u.Quantity, i: u.Quantity) -> u.Quantity:
-    """Compute binary mass function from masses and inclination
+    r"""Compute binary mass function from masses and inclination
 
     Can handle scalar or array inputs.
 
@@ -374,7 +374,7 @@ def mass_funct2(mp: u.Quantity, mc: u.Quantity, i: u.Quantity) -> u.Quantity:
     Calculates
 
     .. math::
-        f(m_p, m_c) = \\frac{m_c^3\sin^3 i}{(m_c + m_p)^2}
+        f(m_p, m_c) = \frac{m_c^3\sin^3 i}{(m_c + m_p)^2}
 
     See [2]_
 
@@ -388,7 +388,7 @@ def mass_funct2(mp: u.Quantity, mc: u.Quantity, i: u.Quantity) -> u.Quantity:
 def pulsar_mass(
     pb: u.Quantity, x: u.Quantity, mc: u.Quantity, i: u.Quantity
 ) -> u.Quantity:
-    """Compute pulsar mass from orbital parameters
+    r"""Compute pulsar mass from orbital parameters
 
     Return the pulsar mass (in solar mass units) for a binary.
     Can handle scalar or array inputs.
@@ -460,7 +460,7 @@ def companion_mass(
     i: u.Quantity = 60.0 * u.deg,
     mp: u.Quantity = 1.4 * u.solMass,
 ) -> u.Quantity:
-    """Commpute the companion mass from the orbital parameters
+    r"""Commpute the companion mass from the orbital parameters
 
     Compute companion mass for a binary system from orbital mechanics,
     not Shapiro delay.
@@ -503,9 +503,9 @@ def companion_mass(
     :math:`a M_c^3 + b M_c^2 + c M_c + d = 0`
 
     - :math:`a = \sin^3(inc)`
-    - :math:`b = -{\\rm massfunct}`
+    - :math:`b = -{\rm massfunct}`
     - :math:`c = -2 M_p {\\rm massfunct}`
-    - :math:`d = -{\\rm massfunct} M_p^2`
+    - :math:`d = -{\rm massfunct} M_p^2`
 
     To solve it we can use a direct calculation of the cubic roots [3]_.
 
@@ -561,7 +561,7 @@ def companion_mass(
 def pbdot(
     mp: u.Quantity, mc: u.Quantity, pb: u.Quantity, e: Union[float, u.Quantity]
 ) -> u.Quantity:
-    """Post-Keplerian orbital decay pbdot, assuming general relativity.
+    r"""Post-Keplerian orbital decay pbdot, assuming general relativity.
 
     pbdot (:math:`\dot P_B`) is the change in the binary orbital period
     due to emission of gravitational waves.
@@ -595,13 +595,13 @@ def pbdot(
     Calculates
 
     .. math::
-        \dot P_b = -\\frac{192\pi}{5}T_{\odot}^{5/3} \\left(\\frac{P_b}{2\pi}\\right)^{-5/3}
-        f(e)\\frac{m_p m_c}{(m_p+m_c)^{1/3}}
+        \dot P_b = -\frac{192\pi}{5}T_{\odot}^{5/3} \left(\frac{P_b}{2\pi}\right)^{-5/3}
+        f(e)\frac{m_p m_c}{(m_p+m_c)^{1/3}}
 
     with
 
     .. math::
-        f(e)=\\frac{1+(73/24)e^2+(37/96)e^4}{(1-e^2)^{7/2}}
+        f(e)=\frac{1+(73/24)e^2+(37/96)e^4}{(1-e^2)^{7/2}}
 
     and :math:`T_\odot = GM_\odot c^{-3}`.
 
@@ -629,11 +629,11 @@ def gamma(
     pb: u.Quantity,
     e: Union[float, u.Quantity],
 ) -> u.Quantity:
-    """Post-Keplerian time dilation and gravitational redshift gamma, assuming general relativity.
+    r"""Post-Keplerian time dilation and gravitational redshift gamma, assuming general relativity.
 
-    gamma (:math:`\\gamma`) is the amplitude of the modification in arrival times caused by the varying
+    gamma (:math:`\gamma`) is the amplitude of the modification in arrival times caused by the varying
     gravitational redshift of the companion and time dilation in an elliptical orbit.  The time delay is
-    :math:`\\gamma \sin E`, where :math:`E` is the eccentric anomaly.
+    :math:`\gamma \sin E`, where :math:`E` is the eccentric anomaly.
     Can handle scalar or array inputs.
 
     Parameters
@@ -664,7 +664,7 @@ def gamma(
     Calculates
 
     .. math::
-        \\gamma = T_{\odot}^{2/3} \\left(\\frac{P_b}{2\pi}\\right)^{1/3} e \\frac{m_c(m_p+2m_c)}{(m_p+m_c)^{4/3}}
+        \gamma = T_{\odot}^{2/3} \left(\frac{P_b}{2\pi}\right)^{1/3} e \frac{m_c(m_p+2m_c)}{(m_p+m_c)^{4/3}}
 
     with :math:`T_\odot = GM_\odot c^{-3}`.
 
@@ -690,7 +690,7 @@ def omdot(
     pb: u.Quantity,
     e: Union[float, u.Quantity],
 ) -> u.Quantity:
-    """Post-Keplerian longitude of periastron precession rate omdot, assuming general relativity.
+    r"""Post-Keplerian longitude of periastron precession rate omdot, assuming general relativity.
 
     omdot (:math:`\dot \omega`) is the relativistic advance of periastron.
     Can handle scalar or array inputs.
@@ -724,8 +724,8 @@ def omdot(
 
     .. math::
 
-        \dot \omega = 3T_{\odot}^{2/3} \\left(\\frac{P_b}{2\pi}\\right)^{-5/3}
-        \\frac{1}{1-e^2}(m_p+m_c)^{2/3}
+        \dot \omega = 3T_{\odot}^{2/3} \left(\frac{P_b}{2\pi}\right)^{-5/3}
+        \frac{1}{1-e^2}(m_p+m_c)^{2/3}
 
     with :math:`T_\odot = GM_\odot c^{-3}`.
 
@@ -745,7 +745,7 @@ def omdot(
 
 @u.quantity_input(mp=u.Msun, mc=u.Msun, pb=u.d, x=u.cm)
 def sini(mp: u.Quantity, mc: u.Quantity, pb: u.Quantity, x: u.Quantity) -> u.Quantity:
-    """Post-Keplerian sine of inclination, assuming general relativity.
+    r"""Post-Keplerian sine of inclination, assuming general relativity.
 
     Can handle scalar or array inputs.
 
@@ -777,8 +777,8 @@ def sini(mp: u.Quantity, mc: u.Quantity, pb: u.Quantity, x: u.Quantity) -> u.Qua
 
     .. math::
 
-        s = T_{\odot}^{-1/3} \\left(\\frac{P_b}{2\pi}\\right)^{-2/3}
-        \\frac{(m_p+m_c)^{2/3}}{m_c}
+        s = T_{\odot}^{-1/3} \left(\frac{P_b}{2\pi}\right)^{-2/3}
+        \frac{(m_p+m_c)^{2/3}}{m_c}
 
     with :math:`T_\odot = GM_\odot c^{-3}`.
 
@@ -799,7 +799,7 @@ def sini(mp: u.Quantity, mc: u.Quantity, pb: u.Quantity, x: u.Quantity) -> u.Qua
 
 @u.quantity_input(mp=u.Msun, mc=u.Msun, pb=u.d)
 def dr(mp: u.Quantity, mc: u.Quantity, pb: u.Quantity) -> u.Quantity:
-    """Post-Keplerian Roemer delay term
+    r"""Post-Keplerian Roemer delay term
 
     dr (:math:`\delta_r`) is part of the relativistic deformation of the orbit
 
@@ -829,8 +829,8 @@ def dr(mp: u.Quantity, mc: u.Quantity, pb: u.Quantity) -> u.Quantity:
 
     .. math::
 
-        \delta_r = T_{\odot}^{2/3} \\left(\\frac{P_b}{2\pi}\\right)^{2/3}
-        \\frac{3 m_p^2+6 m_p m_c +2m_c^2}{(m_p+m_c)^{4/3}}
+        \delta_r = T_{\odot}^{2/3} \left(\frac{P_b}{2\pi}\right)^{2/3}
+        \frac{3 m_p^2+6 m_p m_c +2m_c^2}{(m_p+m_c)^{4/3}}
 
     with :math:`T_\odot = GM_\odot c^{-3}`.
 
@@ -849,9 +849,9 @@ def dr(mp: u.Quantity, mc: u.Quantity, pb: u.Quantity) -> u.Quantity:
 
 @u.quantity_input(mp=u.Msun, mc=u.Msun, pb=u.d)
 def dth(mp: u.Quantity, mc: u.Quantity, pb: u.Quantity) -> u.Quantity:
-    """Post-Keplerian Roemer delay term
+    r"""Post-Keplerian Roemer delay term
 
-    dth (:math:`\delta_{\\theta}`) is part of the relativistic deformation of the orbit
+    dth (:math:`\delta_{\theta}`) is part of the relativistic deformation of the orbit
 
     Parameters
     ----------
@@ -879,8 +879,8 @@ def dth(mp: u.Quantity, mc: u.Quantity, pb: u.Quantity) -> u.Quantity:
 
     .. math::
 
-        \delta_{\\theta} = T_{\odot}^{2/3} \\left(\\frac{P_b}{2\pi}\\right)^{2/3}
-        \\frac{3.5 m_p^2+6 m_p m_c +2m_c^2}{(m_p+m_c)^{4/3}}
+        \delta_{\theta} = T_{\odot}^{2/3} \left(\frac{P_b}{2\pi}\right)^{2/3}
+        \frac{3.5 m_p^2+6 m_p m_c +2m_c^2}{(m_p+m_c)^{4/3}}
 
     with :math:`T_\odot = GM_\odot c^{-3}`.
 
@@ -901,7 +901,7 @@ def dth(mp: u.Quantity, mc: u.Quantity, pb: u.Quantity) -> u.Quantity:
 def omdot_to_mtot(
     omdot: u.Quantity, pb: u.Quantity, e: Union[float, u.Quantity]
 ) -> u.Quantity:
-    """Determine total mass from Post-Keplerian longitude of periastron precession rate omdot,
+    r"""Determine total mass from Post-Keplerian longitude of periastron precession rate omdot,
     assuming general relativity.
 
     omdot (:math:`\dot \omega`) is the relativistic advance of periastron.  It relates to the total
@@ -935,10 +935,10 @@ def omdot_to_mtot(
 
     .. math::
 
-        \dot \omega = 3T_{\odot}^{2/3} \\left(\\frac{P_b}{2\pi}\\right)^{-5/3}
-        \\frac{1}{1-e^2}(m_p+m_c)^{2/3}
+        \dot \omega = 3T_{\odot}^{2/3} \left(\frac{P_b}{2\pi}\right)^{-5/3}
+        \frac{1}{1-e^2}(m_p+m_c)^{2/3}
 
-    to calculate :math:`m_{\\rm tot} = m_p + m_c`,
+    to calculate :math:`m_{\rm tot} = m_p + m_c`,
     with :math:`T_\odot = GM_\odot c^{-3}`.
 
     More details in :ref:`Timing Models`.  Also see [9]_.
@@ -965,7 +965,7 @@ def omdot_to_mtot(
 def a1sini(
     mp: u.Quantity, mc: u.Quantity, pb: u.Quantity, i: u.Quantity = 90 * u.deg
 ) -> u.Quantity:
-    """Pulsar's semi-major axis.
+    r"""Pulsar's semi-major axis.
 
     The full semi-major axis is given by Kepler's third law.  This is the
     projection (:math:`\sin i`) of just the pulsar's orbit (:math:`m_c/(m_p+m_c)`
@@ -1001,8 +1001,8 @@ def a1sini(
 
     .. math::
 
-        \\frac{a_p \sin i}{c} = \\frac{m_c \sin i}{(m_p+m_c)^{2/3}}
-        G^{1/3}\\left(\\frac{P_b}{2\pi}\\right)^{2/3}
+        \frac{a_p \sin i}{c} = \frac{m_c \sin i}{(m_p+m_c)^{2/3}}
+        G^{1/3}\left(\frac{P_b}{2\pi}\right)^{2/3}
 
     More details in :ref:`Timing Models`.  Also see [8]_
 
@@ -1017,7 +1017,7 @@ def a1sini(
 
 @u.quantity_input(pmtot=u.mas / u.yr, D=u.kpc)
 def shklovskii_factor(pmtot: u.Quantity, D: u.Quantity) -> u.Quantity:
-    """Compute magnitude of Shklovskii correction factor.
+    r"""Compute magnitude of Shklovskii correction factor.
 
     Computes the Shklovskii correction factor, as defined in Eq 8.12 of Lorimer & Kramer (2005) [10]_
     This is the factor by which :math:`\dot P /P` is increased due to the transverse velocity.
@@ -1026,7 +1026,7 @@ def shklovskii_factor(pmtot: u.Quantity, D: u.Quantity) -> u.Quantity:
 
     .. math::
 
-        \dot P_{\\rm intrinsic} = \dot P_{\\rm observed} - a_s P
+        \dot P_{\rm intrinsic} = \dot P_{\rm observed} - a_s P
 
     Parameters
     ----------

From 98d31fc1cb1b7c1aef35b7cd846d23a7cca51546 Mon Sep 17 00:00:00 2001
From: David Kaplan <kaplan@uwm.edu>
Date: Tue, 9 Jul 2024 10:30:53 -0500
Subject: [PATCH 06/11] trying to avoid black errors

---
 src/pint/derived_quantities.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/pint/derived_quantities.py b/src/pint/derived_quantities.py
index 5c82ec9aa..83cedddae 100644
--- a/src/pint/derived_quantities.py
+++ b/src/pint/derived_quantities.py
@@ -37,7 +37,7 @@
 def p_to_f(
     p: u.Quantity, pd: u.Quantity, pdd: Optional[u.Quantity] = None
 ) -> Tuple[u.Quantity]:
-    """Converts P, Pdot to F, Fdot (or vice versa)
+    r"""Converts P, Pdot to F, Fdot (or vice versa)
 
     Convert period, period derivative and period second
     derivative (if supplied) to the equivalent frequency counterparts.
@@ -504,7 +504,7 @@ def companion_mass(
 
     - :math:`a = \sin^3(inc)`
     - :math:`b = -{\rm massfunct}`
-    - :math:`c = -2 M_p {\\rm massfunct}`
+    - :math:`c = -2 M_p {\rm massfunct}`
     - :math:`d = -{\rm massfunct} M_p^2`
 
     To solve it we can use a direct calculation of the cubic roots [3]_.

From 52fb27c57bd382b98308c7ddf6e8e5a9d78afc8c Mon Sep 17 00:00:00 2001
From: David Kaplan <kaplan@uwm.edu>
Date: Tue, 9 Jul 2024 10:32:23 -0500
Subject: [PATCH 07/11] slight black version change?

---
 src/pint/derived_quantities.py | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/src/pint/derived_quantities.py b/src/pint/derived_quantities.py
index 83cedddae..6c267c2c7 100644
--- a/src/pint/derived_quantities.py
+++ b/src/pint/derived_quantities.py
@@ -132,7 +132,8 @@ def pferrs(
         return [1.0 / porf, porferr / porf**2.0]
     forperr = porferr / porf**2.0
     fdorpderr = np.sqrt(
-        (4.0 * pdorfd**2.0 * porferr**2.0) / porf**6.0 + pdorfderr**2.0 / porf**4.0
+        (4.0 * pdorfd**2.0 * porferr**2.0) / porf**6.0
+        + pdorfderr**2.0 / porf**4.0
     )
     [forp, fdorpd] = p_to_f(porf, pdorfd)
     return (forp, forperr, fdorpd, fdorpderr)
@@ -537,12 +538,17 @@ def companion_mass(
     # delta1 is always <0
     # delta1 = 2 * b ** 3 - 9 * a * b * c + 27 * a ** 2 * d
     delta1 = (
-        -2 * massfunct**3 - 18 * a * mp * massfunct**2 - 27 * a**2 * massfunct * mp**2
+        -2 * massfunct**3
+        - 18 * a * mp * massfunct**2
+        - 27 * a**2 * massfunct * mp**2
     )
     # Q**2 is always > 0, so this is never a problem
     # in terms of complex numbers
     # Q = np.sqrt(delta1**2 - 4*delta0**3)
-    Q = np.sqrt(108 * a**3 * mp**3 * massfunct**3 + 729 * a**4 * mp**4 * massfunct**2)
+    Q = np.sqrt(
+        108 * a**3 * mp**3 * massfunct**3
+        + 729 * a**4 * mp**4 * massfunct**2
+    )
     # this could be + or - Q
     # pick the - branch since delta1 is <0 so that delta1 - Q is never near 0
     Ccubed = 0.5 * (delta1 + Q)

From 56d1f66b1237cc8b0a62936fe558e70053efccc0 Mon Sep 17 00:00:00 2001
From: David Kaplan <kaplan@uwm.edu>
Date: Thu, 11 Jul 2024 10:03:02 -0500
Subject: [PATCH 08/11] changed return types to units rather than classes for
 compatibility

---
 src/pint/derived_quantities.py | 54 +++++++++++++++-------------------
 1 file changed, 24 insertions(+), 30 deletions(-)

diff --git a/src/pint/derived_quantities.py b/src/pint/derived_quantities.py
index 6c267c2c7..a15ac86ac 100644
--- a/src/pint/derived_quantities.py
+++ b/src/pint/derived_quantities.py
@@ -132,8 +132,7 @@ def pferrs(
         return [1.0 / porf, porferr / porf**2.0]
     forperr = porferr / porf**2.0
     fdorpderr = np.sqrt(
-        (4.0 * pdorfd**2.0 * porferr**2.0) / porf**6.0
-        + pdorfderr**2.0 / porf**4.0
+        (4.0 * pdorfd**2.0 * porferr**2.0) / porf**6.0 + pdorfderr**2.0 / porf**4.0
     )
     [forp, fdorpd] = p_to_f(porf, pdorfd)
     return (forp, forperr, fdorpd, fdorpderr)
@@ -142,7 +141,7 @@ def pferrs(
 @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
-) -> u.Quantity:
+) -> u.yr:
     r"""Compute pulsar characteristic age
 
     Return the age of a pulsar given the spin frequency
@@ -187,7 +186,7 @@ def pulsar_age(
 @u.quantity_input(f=u.Hz, fdot=u.Hz / u.s, I=u.g * u.cm**2)
 def pulsar_edot(
     f: u.Quantity, fdot: u.Quantity, I: u.Quantity = 1.0e45 * u.g * u.cm**2
-) -> u.Quantity:
+) -> u.erg / u.s:
     r"""Compute pulsar spindown energy loss rate
 
     Return the pulsar `Edot` (:math:`\dot E`, in erg/s) given the spin frequency `f` and
@@ -223,7 +222,7 @@ def pulsar_edot(
 
 
 @u.quantity_input(f=u.Hz, fdot=u.Hz / u.s)
-def pulsar_B(f: u.Quantity, fdot: u.Quantity) -> u.Quantity:
+def pulsar_B(f: u.Quantity, fdot: u.Quantity) -> u.G:
     r"""Compute pulsar surface magnetic field
 
     Return the estimated pulsar surface magnetic field strength
@@ -258,7 +257,7 @@ def pulsar_B(f: u.Quantity, fdot: u.Quantity) -> u.Quantity:
 
 
 @u.quantity_input(f=u.Hz, fdot=u.Hz / u.s)
-def pulsar_B_lightcyl(f: u.Quantity, fdot: u.Quantity) -> u.Quantity:
+def pulsar_B_lightcyl(f: u.Quantity, fdot: u.Quantity) -> u.G:
     r"""Compute pulsar magnetic field at the light cylinder
 
     Return the estimated pulsar magnetic field strength at the
@@ -300,7 +299,7 @@ def pulsar_B_lightcyl(f: u.Quantity, fdot: u.Quantity) -> u.Quantity:
 
 
 @u.quantity_input(pb=u.d, x=u.cm)
-def mass_funct(pb: u.Quantity, x: u.Quantity) -> u.Quantity:
+def mass_funct(pb: u.Quantity, x: u.Quantity) -> u.Msun:
     r"""Compute binary mass function from period and semi-major axis
 
     Can handle scalar or array inputs.
@@ -341,7 +340,7 @@ def mass_funct(pb: u.Quantity, x: u.Quantity) -> u.Quantity:
 
 
 @u.quantity_input(mp=u.Msun, mc=u.Msun, i=u.deg)
-def mass_funct2(mp: u.Quantity, mc: u.Quantity, i: u.Quantity) -> u.Quantity:
+def mass_funct2(mp: u.Quantity, mc: u.Quantity, i: u.Quantity) -> u.Msun:
     r"""Compute binary mass function from masses and inclination
 
     Can handle scalar or array inputs.
@@ -386,9 +385,7 @@ def mass_funct2(mp: u.Quantity, mc: u.Quantity, i: u.Quantity) -> u.Quantity:
 
 
 @u.quantity_input(pb=u.d, x=u.cm, mc=u.Msun, i=u.deg)
-def pulsar_mass(
-    pb: u.Quantity, x: u.Quantity, mc: u.Quantity, i: u.Quantity
-) -> u.Quantity:
+def pulsar_mass(pb: u.Quantity, x: u.Quantity, mc: u.Quantity, i: u.Quantity) -> u.Msun:
     r"""Compute pulsar mass from orbital parameters
 
     Return the pulsar mass (in solar mass units) for a binary.
@@ -460,7 +457,7 @@ def companion_mass(
     x: u.Quantity,
     i: u.Quantity = 60.0 * u.deg,
     mp: u.Quantity = 1.4 * u.solMass,
-) -> u.Quantity:
+) -> u.Msun:
     r"""Commpute the companion mass from the orbital parameters
 
     Compute companion mass for a binary system from orbital mechanics,
@@ -538,17 +535,12 @@ def companion_mass(
     # delta1 is always <0
     # delta1 = 2 * b ** 3 - 9 * a * b * c + 27 * a ** 2 * d
     delta1 = (
-        -2 * massfunct**3
-        - 18 * a * mp * massfunct**2
-        - 27 * a**2 * massfunct * mp**2
+        -2 * massfunct**3 - 18 * a * mp * massfunct**2 - 27 * a**2 * massfunct * mp**2
     )
     # Q**2 is always > 0, so this is never a problem
     # in terms of complex numbers
     # Q = np.sqrt(delta1**2 - 4*delta0**3)
-    Q = np.sqrt(
-        108 * a**3 * mp**3 * massfunct**3
-        + 729 * a**4 * mp**4 * massfunct**2
-    )
+    Q = np.sqrt(108 * a**3 * mp**3 * massfunct**3 + 729 * a**4 * mp**4 * massfunct**2)
     # this could be + or - Q
     # pick the - branch since delta1 is <0 so that delta1 - Q is never near 0
     Ccubed = 0.5 * (delta1 + Q)
@@ -566,7 +558,7 @@ def companion_mass(
 @u.quantity_input(mp=u.Msun, mc=u.Msun, pb=u.d, e=u.dimensionless_unscaled)
 def pbdot(
     mp: u.Quantity, mc: u.Quantity, pb: u.Quantity, e: Union[float, u.Quantity]
-) -> u.Quantity:
+) -> u.dimensionless_unscaled:
     r"""Post-Keplerian orbital decay pbdot, assuming general relativity.
 
     pbdot (:math:`\dot P_B`) is the change in the binary orbital period
@@ -634,7 +626,7 @@ def gamma(
     mc: u.Quantity,
     pb: u.Quantity,
     e: Union[float, u.Quantity],
-) -> u.Quantity:
+) -> u.s:
     r"""Post-Keplerian time dilation and gravitational redshift gamma, assuming general relativity.
 
     gamma (:math:`\gamma`) is the amplitude of the modification in arrival times caused by the varying
@@ -695,7 +687,7 @@ def omdot(
     mc: u.Quantity,
     pb: u.Quantity,
     e: Union[float, u.Quantity],
-) -> u.Quantity:
+) -> u.deg / u.yr:
     r"""Post-Keplerian longitude of periastron precession rate omdot, assuming general relativity.
 
     omdot (:math:`\dot \omega`) is the relativistic advance of periastron.
@@ -750,7 +742,9 @@ def omdot(
 
 
 @u.quantity_input(mp=u.Msun, mc=u.Msun, pb=u.d, x=u.cm)
-def sini(mp: u.Quantity, mc: u.Quantity, pb: u.Quantity, x: u.Quantity) -> u.Quantity:
+def sini(
+    mp: u.Quantity, mc: u.Quantity, pb: u.Quantity, x: u.Quantity
+) -> u.dimensionless_unscaled:
     r"""Post-Keplerian sine of inclination, assuming general relativity.
 
     Can handle scalar or array inputs.
@@ -804,7 +798,7 @@ def sini(mp: u.Quantity, mc: u.Quantity, pb: u.Quantity, x: u.Quantity) -> u.Qua
 
 
 @u.quantity_input(mp=u.Msun, mc=u.Msun, pb=u.d)
-def dr(mp: u.Quantity, mc: u.Quantity, pb: u.Quantity) -> u.Quantity:
+def dr(mp: u.Quantity, mc: u.Quantity, pb: u.Quantity) -> u.dimensionless_angles:
     r"""Post-Keplerian Roemer delay term
 
     dr (:math:`\delta_r`) is part of the relativistic deformation of the orbit
@@ -854,7 +848,7 @@ def dr(mp: u.Quantity, mc: u.Quantity, pb: u.Quantity) -> u.Quantity:
 
 
 @u.quantity_input(mp=u.Msun, mc=u.Msun, pb=u.d)
-def dth(mp: u.Quantity, mc: u.Quantity, pb: u.Quantity) -> u.Quantity:
+def dth(mp: u.Quantity, mc: u.Quantity, pb: u.Quantity) -> u.dimensionless_unscaled:
     r"""Post-Keplerian Roemer delay term
 
     dth (:math:`\delta_{\theta}`) is part of the relativistic deformation of the orbit
@@ -906,7 +900,7 @@ def dth(mp: u.Quantity, mc: u.Quantity, pb: u.Quantity) -> u.Quantity:
 @u.quantity_input(omdot=u.deg / u.yr, pb=u.d, e=u.dimensionless_unscaled)
 def omdot_to_mtot(
     omdot: u.Quantity, pb: u.Quantity, e: Union[float, u.Quantity]
-) -> u.Quantity:
+) -> u.Msun:
     r"""Determine total mass from Post-Keplerian longitude of periastron precession rate omdot,
     assuming general relativity.
 
@@ -967,10 +961,10 @@ def omdot_to_mtot(
     ).to(u.Msun, equivalencies=u.dimensionless_angles())
 
 
-@u.quantity_input(pb=u.d, mp=u.Msun, mc=u.Msun, i=u.deg)
+@u.quantity_input(mp=u.Msun, mc=u.Msun, pb=u.d, i=u.deg)
 def a1sini(
     mp: u.Quantity, mc: u.Quantity, pb: u.Quantity, i: u.Quantity = 90 * u.deg
-) -> u.Quantity:
+) -> pint.ls:
     r"""Pulsar's semi-major axis.
 
     The full semi-major axis is given by Kepler's third law.  This is the
@@ -1022,7 +1016,7 @@ def a1sini(
 
 
 @u.quantity_input(pmtot=u.mas / u.yr, D=u.kpc)
-def shklovskii_factor(pmtot: u.Quantity, D: u.Quantity) -> u.Quantity:
+def shklovskii_factor(pmtot: u.Quantity, D: u.Quantity) -> 1 / u.s:
     r"""Compute magnitude of Shklovskii correction factor.
 
     Computes the Shklovskii correction factor, as defined in Eq 8.12 of Lorimer & Kramer (2005) [10]_
@@ -1060,7 +1054,7 @@ def shklovskii_factor(pmtot: u.Quantity, D: u.Quantity) -> u.Quantity:
 
 
 @u.quantity_input(dm=pint.dmu)
-def dispersion_slope(dm: u.Quantity) -> u.Quantity:
+def dispersion_slope(dm: u.Quantity) -> 1 / u.s:
     """Compute the dispersion slope.
 
     This is equal to DMconst * DM.

From 1ae4557f9b6b595e4b4a2dc5e124751b5dc32a40 Mon Sep 17 00:00:00 2001
From: David Kaplan <kaplan@uwm.edu>
Date: Thu, 11 Jul 2024 10:08:19 -0500
Subject: [PATCH 09/11] still weird black version differences

---
 src/pint/derived_quantities.py | 12 +++++++++---
 1 file changed, 9 insertions(+), 3 deletions(-)

diff --git a/src/pint/derived_quantities.py b/src/pint/derived_quantities.py
index a15ac86ac..273777354 100644
--- a/src/pint/derived_quantities.py
+++ b/src/pint/derived_quantities.py
@@ -132,7 +132,8 @@ def pferrs(
         return [1.0 / porf, porferr / porf**2.0]
     forperr = porferr / porf**2.0
     fdorpderr = np.sqrt(
-        (4.0 * pdorfd**2.0 * porferr**2.0) / porf**6.0 + pdorfderr**2.0 / porf**4.0
+        (4.0 * pdorfd**2.0 * porferr**2.0) / porf**6.0
+        + pdorfderr**2.0 / porf**4.0
     )
     [forp, fdorpd] = p_to_f(porf, pdorfd)
     return (forp, forperr, fdorpd, fdorpderr)
@@ -535,12 +536,17 @@ def companion_mass(
     # delta1 is always <0
     # delta1 = 2 * b ** 3 - 9 * a * b * c + 27 * a ** 2 * d
     delta1 = (
-        -2 * massfunct**3 - 18 * a * mp * massfunct**2 - 27 * a**2 * massfunct * mp**2
+        -2 * massfunct**3
+        - 18 * a * mp * massfunct**2
+        - 27 * a**2 * massfunct * mp**2
     )
     # Q**2 is always > 0, so this is never a problem
     # in terms of complex numbers
     # Q = np.sqrt(delta1**2 - 4*delta0**3)
-    Q = np.sqrt(108 * a**3 * mp**3 * massfunct**3 + 729 * a**4 * mp**4 * massfunct**2)
+    Q = np.sqrt(
+        108 * a**3 * mp**3 * massfunct**3
+        + 729 * a**4 * mp**4 * massfunct**2
+    )
     # this could be + or - Q
     # pick the - branch since delta1 is <0 so that delta1 - Q is never near 0
     Ccubed = 0.5 * (delta1 + Q)

From f87f75f46b5d81732cbfe6e8ddd3e8e5b21687ea Mon Sep 17 00:00:00 2001
From: David Kaplan <kaplan@uwm.edu>
Date: Thu, 11 Jul 2024 11:01:27 -0500
Subject: [PATCH 10/11] removed some output types

---
 src/pint/derived_quantities.py | 8 +++-----
 1 file changed, 3 insertions(+), 5 deletions(-)

diff --git a/src/pint/derived_quantities.py b/src/pint/derived_quantities.py
index 273777354..e1e45cc55 100644
--- a/src/pint/derived_quantities.py
+++ b/src/pint/derived_quantities.py
@@ -4,7 +4,7 @@
 import astropy.constants as const
 import astropy.units as u
 import numpy as np
-from typing import Optional, List, Tuple, Union
+from typing import Optional, Union
 
 import pint
 
@@ -34,9 +34,7 @@
 @u.quantity_input(
     p=[u.Hz, u.s], pd=[u.Hz / u.s, u.s / u.s], pdd=[u.Hz / u.s**2, u.s / u.s**2]
 )
-def p_to_f(
-    p: u.Quantity, pd: u.Quantity, pdd: Optional[u.Quantity] = None
-) -> Tuple[u.Quantity]:
+def p_to_f(p: u.Quantity, pd: u.Quantity, pdd: Optional[u.Quantity] = None):
     r"""Converts P, Pdot to F, Fdot (or vice versa)
 
     Convert period, period derivative and period second
@@ -92,7 +90,7 @@ def pferrs(
     porferr: u.Quantity,
     pdorfd: Optional[u.Quantity] = None,
     pdorfderr: Optional[u.Quantity] = None,
-) -> Tuple[u.Quantity]:
+):
     r"""Convert P, Pdot to F, Fdot with uncertainties (or vice versa).
 
     Calculate the period or frequency errors and

From 80528f4a712a1221e7ed4b5359fe4cfe7c5e80f1 Mon Sep 17 00:00:00 2001
From: David Kaplan <kaplan@uwm.edu>
Date: Thu, 11 Jul 2024 11:58:11 -0500
Subject: [PATCH 11/11] fixed output units

---
 src/pint/derived_quantities.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/pint/derived_quantities.py b/src/pint/derived_quantities.py
index e1e45cc55..5d69de328 100644
--- a/src/pint/derived_quantities.py
+++ b/src/pint/derived_quantities.py
@@ -802,7 +802,7 @@ def sini(
 
 
 @u.quantity_input(mp=u.Msun, mc=u.Msun, pb=u.d)
-def dr(mp: u.Quantity, mc: u.Quantity, pb: u.Quantity) -> u.dimensionless_angles:
+def dr(mp: u.Quantity, mc: u.Quantity, pb: u.Quantity) -> u.dimensionless_unscaled:
     r"""Post-Keplerian Roemer delay term
 
     dr (:math:`\delta_r`) is part of the relativistic deformation of the orbit