Repository: PlasmaPy/PlasmaPy

Issue:
Create a `RelativisticBody` class
Right now it's pretty inconvenient to have a particle and go between velocity, Lorentz factor, energy, and momentum.  It'd be helpful to have a `RelativisticParticle` class that we could use to take one such quantity and make it straightforward to calculate the rest.  I was thinking that this would be best in `plasmapy.formulary.relativistic`.

Some of the input parameters could be `particle`, `V`, `gamma` (for Lorentz factor), `kinetic_energy`, and `momentum`. Some of the attributes could be `rest_mass`, `relativistic_mass`, `total_energy`, `kinetic_energy`, and `mass_energy`.  There could also be an attribute like `v_over_c`, which is often referred to as `beta`, but that could potentially be confused with plasma beta.  

I'm wondering if we could even define the dunder methods for `__add__` and `__subtract__` so that we could do things like...
```Python
>>> cosmic_ray = RelativisticParticle('p+', energy = 100 * u.GeV)
>>> cosmic_ray += 50 * u.GeV  # would increase the energy of the particle by 50 GeV
```

Doing this would make #1538 a lot easier, and it's related to #1537 too.

Discussion:
codecov[bot]: # [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1540?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy) Report
> Merging [#1540](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1540?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy) (7ee192a) into [main](https://codecov.io/gh/PlasmaPy/PlasmaPy/commit/fb7c08ab37cdfd3b336627b89c4b3ead21f1d518?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy) (fb7c08a) will **increase** coverage by `0.00%`.
> The diff coverage is `97.82%`.

> :exclamation: Current head 7ee192a differs from pull request most recent head 3a2ea24. Consider uploading reports for the commit 3a2ea24 to get more accurate results

```diff
@@           Coverage Diff           @@
##             main    #1540   +/-   ##
=======================================
  Coverage   97.09%   97.09%           
=======================================
  Files          82       82           
  Lines        7666     7711   +45     
=======================================
+ Hits         7443     7487   +44     
- Misses        223      224    +1     
```


| [Impacted Files](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1540?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy) | Coverage Δ | |
|---|---|---|
| [plasmapy/formulary/relativity.py](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1540/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy#diff-cGxhc21hcHkvZm9ybXVsYXJ5L3JlbGF0aXZpdHkucHk=) | `98.50% <97.82%> (-1.50%)` | :arrow_down: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1540?src=pr&el=continue&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1540?src=pr&el=footer&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy). Last update [fb7c08a...3a2ea24](https://codecov.io/gh/PlasmaPy/PlasmaPy/pull/1540?src=pr&el=lastupdated&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=PlasmaPy).

namurphy: The essential functionality seems to be working now, but there are still a bunch of things that need to be done.

 - [ ] Allow `Quantity` arrays as inputs (needed for #1538)
 - [ ] Allow `ParticleList` inputs
 - [ ] Implement `__len__`
 - [ ] Implement `particle_input`
 - [ ] Implement `validate_quantities`
 - [ ] Require `v` to be nonnegative (?)
 - [ ] Implement `@check_relativistic`
 - [ ] Finish `__str__` and `__repr__`
 - [ ] Add `__gt__`, etc.
 - [ ] Add tests for setting attributes
 - [ ] Add tests for `Quantity` array and `ParticleList` inputs

StanczakDominik: What would you like `len` and `gt` to do here?

This reminded me that the fact we implemented gt on something else (particle?) iirc broke sorting of collections of that object. We'll have to check if that is actually the case, and if that would be the case here.
namurphy: > What would you like `len` and `gt` to do here?

Good question. In retrospect, it probably makes sense not to implement these. If we allow `RelativisticBody` to contain multiple bodies (i.e. accept a `ParticleList`) or alternatively multiple velocities, then we can implement `__len__`.  But, I don't think that'll be this PR.
namurphy: This was trickier to write than I thought it would be!  It's just about ready for review now. Some thoughts...

 - I ended up adding `dtype` as a parameter to `RelativisticBody` to account for really relativistic particles.  The default dtype is the highest precision option for NumPy, but unfortunately that's still not enough to avoid roundoff error for ultra-high-energy cosmic rays for `velocity` and `v_over_c`.  I put in an admonition.  I'm not sure if it would be worth using mpmath to account for those cases.
 - I bumped the minimum version of NumPy to 1.20, which is consistent with [NEP 29](https://numpy.org/neps/nep-0029-deprecation_policy.html) and will enable us to use `numpy.typing`, like with `DTypeLike` here.
 - I didn't add a method to calculate the gyroradius yet.  I'm still thinking about how to best handle it.  The gyroradius depends on the angle between the velocity and the magnetic field, so there would be an implicit assumption that $V ⟂ B$.
 - I used `V` as the parameter name instead of `velocity`, but `velocity` as the attribute name.  I'm tempted to use `velocity` as the parameter name too, though this would violate our convention of using `V`.  
namurphy: @pheuer — might you be interested in reviewing this PR?  Thanks!
namurphy: It turns out that data types for `numpy.longdouble`, etc. are platform dependent. I think that was causing one of my tests to fail on Windows for V/c for an ultra-high energy cosmic ray that's an electron. I updated the test to be a proton so the velocity is slower.
namurphy: I take pride<sup>1</sup> in my procrastination driven development!

#### Footnotes


<sup>1</sup> Also anxiety.

Patch:
diff --git a/changelog/1540.feature.rst b/changelog/1540.feature.rst
new file mode 100644
index 00000000..3511a360
--- /dev/null
+++ b/changelog/1540.feature.rst
@@ -0,0 +1,2 @@
+Added the |RelativisticBody| class to facilitate calculation of the
+relativistic properties of a body in motion.
diff --git a/changelog/1540.trivial.rst b/changelog/1540.trivial.rst
new file mode 100644
index 00000000..cd126945
--- /dev/null
+++ b/changelog/1540.trivial.rst
@@ -0,0 +1,1 @@
+Updated the minimum required version of NumPy to 1.20.0.
diff --git a/docs/common_links.rst b/docs/common_links.rst
index e447d85e..559ca323 100644
--- a/docs/common_links.rst
+++ b/docs/common_links.rst
@@ -6,6 +6,7 @@
 .. ------------------
 
 .. |ClassicalTransport| replace:: :class:`~plasmapy.formulary.braginskii.ClassicalTransport`
+.. |RelativisticBody| replace:: :class:`~plasmapy.formulary.relativity.RelativisticBody`
 
 .. ------------------
 .. plasmapy.particles
@@ -72,6 +73,8 @@
 .. |nan| replace:: `~numpy.nan`
 .. |ndarray| replace:: :class:`~numpy.ndarray`
 .. |array_like| replace:: :term:`numpy:array_like`
+.. |ArrayLike| replace:: `~numpy.typing.ArrayLike`
+.. |DTypeLike| replace:: `~numpy.typing.DTypeLike`
 
 .. --------------------
 .. Astropy replacements
@@ -95,6 +98,7 @@
 .. |contributor guide| replace:: :ref:`contributor guide`\
 .. |documentation guide| replace:: :ref:`documentation guide`\
 .. |glossary| replace:: :ref:`glossary`\
+.. |keyword-only| replace:: :term:`keyword-only`\
 .. |minpython| replace:: 3.8
 .. |plasma-calculator| replace:: :ref:`plasmapy-calculator`\
 .. |release guide| replace:: :ref:`release guide`\
diff --git a/plasmapy/formulary/relativity.py b/plasmapy/formulary/relativity.py
index cc970256..b0a488c4 100644
--- a/plasmapy/formulary/relativity.py
+++ b/plasmapy/formulary/relativity.py
@@ -1,12 +1,18 @@
 r"""Functions for calculating relativistic quantities (:math:`v \to c`)."""
-__all__ = ["Lorentz_factor", "relativistic_energy"]
+__all__ = ["Lorentz_factor", "relativistic_energy", "RelativisticBody"]
 
 import astropy.units as u
 import numpy as np
 
 from astropy.constants import c
+from numbers import Integral, Real
+from numpy.typing import DTypeLike
+from typing import Dict, Optional, Union
 
 from plasmapy import utils
+from plasmapy.particles._factory import _physical_particle_factory
+from plasmapy.particles.particle_class import CustomParticle, Particle, ParticleLike
+from plasmapy.particles.particle_collections import ParticleList
 from plasmapy.utils.decorators import validate_quantities
 
 
@@ -66,7 +72,7 @@ def Lorentz_factor(V: u.m / u.s):
     if not np.all((np.abs(V) <= c) | (np.isnan(V))):
         raise utils.RelativityError(
             "The Lorentz factor cannot be calculated for "
-            "speeds faster than the speed of light. "
+            "speeds faster than the speed of light."
         )
 
     if V.size > 1:
@@ -150,3 +156,351 @@ def relativistic_energy(m: u.kg, v: u.m / u.s) -> u.Joule:
     """
     γ = Lorentz_factor(v)
     return γ * m * c**2
+
+
+class RelativisticBody:
+    r"""
+    A physical body that is moving at a velocity relative to the speed
+    of light.
+
+    Parameters
+    ----------
+    particle : |ParticleLike|, |CustomParticle|, |ParticleList|, or |Quantity|
+        A representation of a particle from which to get the mass
+        of the relativistic body. If it is a |Quantity|, then it must
+        have units of mass and describe the body's rest mass.
+
+    V : |Quantity|, optional
+        The velocity of the relativistic body in units convertible to
+        m/s. The absolute magnitude of ``V`` cannot be greater than
+        :math:`c`\ .
+
+    momentum : |Quantity|, optional
+        The momentum of the relativistic body in units convertible to
+        kg·m/s.
+
+    total_energy : |Quantity|, optional, |keyword-only|
+       The sum of the mass energy and the kinetic energy in units
+       convertible to joules. Must be non-negative.
+
+    kinetic_energy : |Quantity|, optional, |keyword-only|
+       The kinetic energy of the relativistic body in units convertible
+       to joules. Must be non-negative.
+
+    v_over_c : real number or |Quantity|, optional, |keyword-only|
+       The ratio of the velocity to the speed of light. Must have an
+       absolute magnitude :math:`≤ 1`\ .
+
+    lorentz_factor : real number or |Quantity|, optional, |keyword-only|
+       The Lorentz factor of the relativistic body. Must be
+       :math:`≥ 1`\ .
+
+    Z : integer, optional, |keyword-only|
+        The charge number associated with ``particle``.
+
+    mass_numb : integer, optional, |keyword-only|
+        The mass number associated with ``particle``.
+
+    dtype : |DTypeLike|, optional, |keyword-only|, default: `numpy.longdouble`
+        The `numpy` data type to use to store the inputs.
+
+    Notes
+    -----
+    At most one of ``V``, ``momentum``, ``total_energy``,
+    ``kinetic_energy``, ``v_over_c``, and ``lorentz_factor`` must be
+    provided.
+
+    .. caution::
+
+       For ultra-high-energy cosmic rays (UHECRs), the velocity may be
+       within roundoff error of :math:`c` and :math:`\frac{V}{c}` may be
+       within roundoff error of 1.
+
+    Examples
+    --------
+    >>> import astropy.units as u
+    >>> relativistic_proton = RelativisticBody("p+", total_energy = 1 * u.GeV)
+    >>> relativistic_proton.particle
+    Particle("p+")
+    >>> relativistic_proton.velocity
+    <Quantity 1.03697...e+08 m / s>
+    >>> relativistic_proton.v_over_c
+    0.3458980898746...
+    >>> relativistic_proton.lorentz_factor
+    1.0657889247888...
+    >>> relativistic_proton.mass_energy.to("GeV")
+    <Quantity 0.93827... GeV>
+    >>> relativistic_proton.total_energy.to("GeV")
+    <Quantity 1. GeV>
+    >>> relativistic_proton.mass
+    <Quantity 1.67262...e-27 kg>
+
+    |RelativisticBody| also works with multiple particles and/or
+    velocities.
+
+    >>> particles = ["p+", "e-"]
+    >>> velocities = [2e5, 2e8] * u.m / u.s
+    >>> relativistic_particles = RelativisticBody(particles, velocities)
+    >>> relativistic_particles.momentum
+    <Quantity [3.345244...e-22, 2.445659...e-22] kg m / s>
+    """
+
+    @staticmethod
+    def _get_speed_like_input(
+        velocity_like_arguments: Dict[str, Union[u.Quantity, Real]]
+    ):
+
+        not_none_arguments = {
+            key: value
+            for key, value in velocity_like_arguments.items()
+            if value is not None
+        }
+
+        if len(not_none_arguments) != 1:
+            raise ValueError(
+                "RelativisticBody can accept no more than one of the following "
+                "arguments: V, v_over_c, momentum, total_energy, kinetic_energy, "
+                "and lorentz_factor."
+            )
+
+        return not_none_arguments or {"velocity": np.nan * u.m / u.s}
+
+    def _store_velocity_like_argument(
+        self, speed_like_input: Dict[str, Union[u.Quantity, Real]]
+    ):
+        """
+        Take the velocity-like argument and store it via the setter for
+        the corresponding attribute.
+        """
+        name = list(speed_like_input.keys())[0]
+        value = speed_like_input[name]
+        if self._dtype:
+            value = u.Quantity(value, dtype=self._dtype)
+        setattr(self, name, value)
+
+    @validate_quantities(
+        V={"can_be_inf": False, "none_shall_pass": True},
+        momentum={"can_be_inf": False, "none_shall_pass": True},
+        total_energy={"can_be_negative": False, "none_shall_pass": True},
+        kinetic_energy={"can_be_negative": False, "none_shall_pass": True},
+    )
+    def __init__(
+        self,
+        particle: Union[ParticleLike, u.Quantity],
+        V: u.m / u.s = None,
+        momentum: u.kg * u.m / u.s = None,
+        *,
+        total_energy: u.J = None,
+        kinetic_energy: u.J = None,
+        v_over_c: Optional[Real] = None,
+        lorentz_factor: Optional[Real] = None,
+        Z: Optional[Integral] = None,
+        mass_numb: Optional[Integral] = None,
+        dtype: Optional[DTypeLike] = np.longdouble,
+    ):
+
+        self._particle = _physical_particle_factory(particle, Z=Z, mass_numb=mass_numb)
+        self._dtype = dtype
+
+        velocity_like_inputs = {
+            "velocity": V,
+            "momentum": momentum,
+            "total_energy": total_energy,
+            "kinetic_energy": kinetic_energy,
+            "v_over_c": v_over_c,
+            "lorentz_factor": lorentz_factor,
+        }
+
+        speed_like_input = self._get_speed_like_input(velocity_like_inputs)
+        self._store_velocity_like_argument(speed_like_input)
+
+    def __repr__(self):
+        return f"RelativisticBody({self.particle}, {self.velocity})"
+
+    @property
+    def particle(self) -> Union[CustomParticle, Particle, ParticleList]:
+        """
+        Representation of the particle(s).
+
+        Returns
+        -------
+        |Particle|, |CustomParticle|, or |ParticleList|
+        """
+        return self._particle
+
+    @property
+    @validate_quantities
+    def mass(self) -> u.kg:
+        r"""
+        The rest mass of the body, :math:`m_0`\ .
+
+        Returns
+        -------
+        ~astropy.units.Quantity
+        """
+        return u.Quantity(self.particle.mass, dtype=self._dtype)
+
+    @property
+    @validate_quantities
+    def mass_energy(self) -> u.J:
+        r"""
+        The rest mass energy of the body, :math:`m_0 c^2`\ .
+
+        Returns
+        -------
+        ~astropy.units.Quantity
+        """
+        return self.mass * c**2
+
+    @property
+    @validate_quantities
+    def total_energy(self) -> u.J:
+        r"""
+        The sum of the rest mass energy and the kinetic energy of the
+        body,
+
+        .. math::
+
+           E_\mathrm{tot} ≡ γ m_0 c^2.
+
+        Returns
+        -------
+        ~astropy.units.Quantity
+        """
+        return np.sqrt(self.momentum**2 * c**2 + self.mass_energy**2)
+
+    @property
+    @validate_quantities
+    def kinetic_energy(self) -> u.J:
+        """
+        The kinetic energy of the body,
+
+        .. math::
+
+           E_K ≡ m_0 c^2 (γ-1).
+
+        Returns
+        -------
+        ~astropy.units.Quantity
+        """
+        return self.total_energy - self.mass_energy
+
+    @property
+    @validate_quantities
+    def v_over_c(self) -> Real:
+        r"""
+        The velocity of the body divided by the velocity of light:
+        :math:`\frac{V}{c}`\ .
+
+        Returns
+        -------
+        float
+        """
+        return (self.velocity / c).to(u.dimensionless_unscaled).value
+
+    @property
+    @validate_quantities
+    def velocity(self) -> u.m / u.s:
+        r"""
+        The velocity of the body, :math:`V`\ .
+
+        Returns
+        -------
+        ~astropy.units.Quantity
+        """
+        velocity = self.momentum / np.sqrt(self.mass**2 + self.momentum**2 / c**2)
+        return velocity.to(u.m / u.s)
+
+    @property
+    @validate_quantities
+    def lorentz_factor(self) -> Real:
+        r"""
+        The Lorentz factor of the body,
+
+        .. math::
+
+           γ ≡ \frac{1}{\sqrt{1 - \frac{V^2}{c^2}}}.
+
+        Returns
+        -------
+        float
+        """
+        return Lorentz_factor(self.velocity)
+
+    @property
+    @validate_quantities
+    def momentum(self) -> u.kg * u.m / u.s:
+        r"""
+        The magnitude of the momentum of the body,
+
+        .. math::
+
+           p ≡ γ m_0 V.
+
+        Returns
+        -------
+        ~astropy.units.Quantity
+        """
+        return getattr(self, "_momentum")
+
+    @kinetic_energy.setter
+    @validate_quantities(E_K={"can_be_negative": False})
+    def kinetic_energy(self, E_K: u.J):
+        self.total_energy = E_K + self.mass_energy
+
+    @total_energy.setter
+    @validate_quantities(E_tot={"can_be_negative": False})
+    def total_energy(self, E_tot: u.J):
+        self._momentum = np.sqrt(E_tot**2 - self.mass_energy**2) / c
+
+    @v_over_c.setter
+    def v_over_c(self, v_over_c_: Real):
+        self.velocity = v_over_c_ * c
+
+    @velocity.setter
+    @validate_quantities
+    def velocity(self, V: u.m / u.s):
+        self._momentum = (Lorentz_factor(V) * self.mass * V).to(u.kg * u.m / u.s)
+
+    @lorentz_factor.setter
+    def lorentz_factor(self, γ: Union[Real, u.Quantity]):
+        if not isinstance(γ, (Real, u.Quantity)):
+            raise TypeError("Invalid type for Lorentz factor")
+
+        if isinstance(γ, u.Quantity):
+            try:
+                γ = γ.to(u.dimensionless_unscaled).value
+            except u.UnitConversionError as exc:
+                raise u.UnitConversionError(
+                    "The Lorentz factor must be dimensionless."
+                ) from exc
+
+        if γ < 1:
+            raise ValueError("The Lorentz factor must be ≥ 1")
+
+        self.velocity = c * np.sqrt(1 - γ**-2)
+
+    @momentum.setter
+    @validate_quantities
+    def momentum(self, p: u.kg * u.m / u.s):
+        self._momentum = p.to(u.kg * u.m / u.s)
+
+    def __eq__(self, other) -> bool:
+
+        _attributes_to_compare = (
+            "particle",
+            "kinetic_energy",
+            "mass_energy",
+            "total_energy",
+            "v_over_c",
+            "momentum",
+        )
+
+        for attr in _attributes_to_compare:
+            if not hasattr(other, attr):
+                return False
+            self_value = getattr(self, attr)
+            other_value = getattr(other, attr)
+            if self_value != other_value:
+                return False
+        return True
diff --git a/plasmapy/particles/_factory.py b/plasmapy/particles/_factory.py
index 0667cf35..89eafe3d 100644
--- a/plasmapy/particles/_factory.py
+++ b/plasmapy/particles/_factory.py
@@ -24,6 +24,12 @@ def _physical_particle_factory(
     This function will select the appropriate type among |Particle|,
     |CustomParticle|, and |ParticleList|.
 
+    .. caution::
+
+       If |Quantity| instances are provided to this function as
+       positional arguments, then they must presently be in the order
+       expected by |CustomParticle|.
+
     Parameters
     ----------
     *args
@@ -49,6 +55,12 @@ def _physical_particle_factory(
     ~plasmapy.particles.particle_class.CustomParticle
     ~plasmapy.particles.particle_class.ParticleList
 
+    Notes
+    -=---
+    If ``Z`` or ``mass_numb`` is provided as keyword arguments but are
+    equal to `None`, then they will not be provided to any of the calls.
+    This is to allow |CustomParticle| instances to be created.
+
     Examples
     --------
     >>> from plasmapy.particles._factory import _physical_particle_factory
@@ -60,6 +72,15 @@ def _physical_particle_factory(
     >>> _physical_particle_factory(["p+", "e-"])
     ParticleList(['p+', 'e-'])
     """
+
+    # We need to remove Z and mass_numb from kwargs when they are `None`
+    # because they are not allowed as arguments to `CustomParticle`, and
+    # are not needed in kwargs if they are their default values. Note
+    # that this affects `not kwargs` below.
+    for parameter in ["Z", "mass_numb"]:
+        if parameter in kwargs and kwargs[parameter] is None:
+            kwargs.pop(parameter)
+
     if (
         len(args) == 1
         and not kwargs


Test patch:
diff --git a/plasmapy/formulary/tests/test_relativity.py b/plasmapy/formulary/tests/test_relativity.py
index 1fd745d1..5e8b1ef3 100644
--- a/plasmapy/formulary/tests/test_relativity.py
+++ b/plasmapy/formulary/tests/test_relativity.py
@@ -5,8 +5,14 @@
 
 from astropy import units as u
 from astropy.constants import c
+from astropy.tests.helper import assert_quantity_allclose
 
-from plasmapy.formulary.relativity import Lorentz_factor, relativistic_energy
+from plasmapy.formulary.relativity import (
+    Lorentz_factor,
+    relativistic_energy,
+    RelativisticBody,
+)
+from plasmapy.particles import CustomParticle, electron, proton
 from plasmapy.utils.exceptions import RelativityError
 
 
@@ -92,3 +98,213 @@ def test_relativistic_energy():
 
     with pytest.raises(ValueError):
         relativistic_energy(-m, v)
+
+
+proton_at_half_c_inputs = [
+    ("v_over_c", 0.5),
+    ("velocity", 0.5 * c),
+    ("lorentz_factor", 1.1547005383792517),
+    ("total_energy", 1.7358354725115025e-10 * u.J),
+    ("kinetic_energy", 2.3255785652637692e-11 * u.J),
+    ("momentum", 2.8950619440057805e-19 * u.kg * u.m / u.s),
+]
+
+
+@pytest.mark.parametrize("attribute, expected", proton_at_half_c_inputs)
+@pytest.mark.parametrize("parameter, argument", proton_at_half_c_inputs)
+def test_relativistic_body(parameter, argument, attribute, expected):
+    """
+    Test that when we create `RelativisticBody` instances for each of
+    the different velocity-like arguments, that each of the resulting
+    `RelativisticBody` attributes end up providing the correct value.
+    """
+
+    if parameter == "velocity":
+        parameter = "V"
+
+    kwargs = {"particle": proton, parameter: argument}
+
+    relativistic_body = RelativisticBody(**kwargs)
+    actual = getattr(relativistic_body, attribute)
+
+    assert_quantity_allclose(actual, expected, rtol=1e-9)
+
+
+@pytest.mark.parametrize("attr_to_set, set_value", proton_at_half_c_inputs)
+@pytest.mark.parametrize("attr_to_test, expected", proton_at_half_c_inputs)
+def test_relativistic_body_setters(attr_to_set, set_value, attr_to_test, expected):
+    """Test setting RelativisticBody attributes."""
+    relativistic_body = RelativisticBody(proton, v_over_c=0.1)
+    setattr(relativistic_body, attr_to_set, set_value)
+
+    actual = getattr(relativistic_body, attr_to_test)
+
+    assert u.isclose(expected, actual, rtol=1e-8), (
+        f"When setting {attr_to_set} to {set_value!r} in a "
+        f"RelativisticBody instance, the value of {attr_to_test} was "
+        f"expected to be {expected!r}, but was instead {actual!r}."
+    )
+
+
+@pytest.mark.parametrize("particle", [electron, proton])
+def test_relativistic_body_mass_energy(particle):
+    """Test `RelativisticBody.mass_energy`."""
+    relativistic_body = RelativisticBody(particle, v_over_c=0)
+    expected = particle.mass * c**2
+
+    actual = relativistic_body.mass_energy
+
+    assert u.isclose(actual, expected, rtol=1e-9)
+
+
+@pytest.mark.parametrize(
+    "kwargs, exception",
+    [
+        ({"V": 299792459 * (u.m / u.s)}, RelativityError),
+        ({"v_over_c": 1.00001}, RelativityError),
+        ({"total_energy": -1 * u.J}, ValueError),
+        ({"kinetic_energy": -1 * u.J}, ValueError),
+        ({"lorentz_factor": 0.99}, ValueError),
+        ({"lorentz_factor": 0}, ValueError),
+        ({"lorentz_factor": -1}, ValueError),
+        ({"lorentz_factor": -1}, ValueError),
+        ({"lorentz_factor": -1, "v_over_c": 0.5}, ValueError),
+        ({"total_energy": 1 * u.J, "momentum": 1 * u.kg * u.m / u.s}, ValueError),
+        ({}, ValueError),
+        ({"lorentz_factor": "wrong type"}, TypeError),
+        ({"lorentz_factor": 3 * u.m / u.s}, u.UnitConversionError),
+    ],
+)
+def test_relativistic_body_init_exceptions(kwargs, exception):
+    """
+    Test that `RelativisticBody` raises the appropriate exceptions
+    during instantiation.
+    """
+    with pytest.raises(exception):
+        RelativisticBody(proton, **kwargs)
+
+
+def test_relativistic_body_equality():
+    """Test that a `RelativisticBody` instance equals itself."""
+    relativistic_body = RelativisticBody(particle=proton, v_over_c=0.34)
+    assert relativistic_body == relativistic_body
+
+
+@pytest.mark.parametrize(
+    "this, that",
+    [
+        (RelativisticBody("p+", v_over_c=0.2), RelativisticBody("p+", v_over_c=0.5)),
+        (RelativisticBody("p+", v_over_c=0.2), RelativisticBody("e-", v_over_c=0.2)),
+        (RelativisticBody("p+", V=c / 2), "different type"),
+    ],
+)
+def test_relativistic_body_inequalities(this, that):
+    """Test the inequality properties of `RelativisticBody`."""
+    assert this != that
+
+
+def test_relativistic_body_inequality_with_different_velocities():
+    """
+    Test that `RelativisticBody` instances are not equal when the
+    velocity provided to them is different.
+    """
+    slower_body = RelativisticBody(particle=proton, v_over_c=0.23)
+    faster_body = RelativisticBody(particle=proton, v_over_c=0.24)
+    assert slower_body != faster_body
+
+
+def test_relativistic_body_inequality_with_different_particles():
+    """
+    Test that `RelativisticBody` instances are not equal when the
+    particles are different.
+    """
+    relativistic_proton = RelativisticBody(particle=proton, v_over_c=0.23)
+    relativistic_electron = RelativisticBody(particle=electron, v_over_c=0.23)
+    assert relativistic_proton != relativistic_electron
+
+
+@pytest.mark.parametrize(
+    "attr",
+    [
+        "lorentz_factor",
+        "v_over_c",
+        "velocity",
+        "momentum",
+        "total_energy",
+        "kinetic_energy",
+    ],
+)
+def test_relativistic_body_defined_using_mass(attr):
+    """Test that a RelativisticBody can be provided a mass as the particle."""
+    V = c / 2
+
+    relativistic_proton = RelativisticBody(proton, V=V)
+    expected = getattr(relativistic_proton, attr)
+
+    relativistic_custom_particle = RelativisticBody(proton.mass, V=V)
+    actual = getattr(relativistic_custom_particle, attr)
+
+    assert u.isclose(actual, expected, rtol=1e-9)
+
+
+@pytest.mark.xfail(reason="RelativisticBody does not yet accept nan velocities")
+def test_relativistic_body_nan_velocity():
+    """
+    Test that RelativisticBody can be created with no velocity defined,
+    and then have the velocity be nan.
+    """
+    relativistic_body = RelativisticBody("p+")
+    assert np.isnan(relativistic_body.velocity)
+
+
+def test_relativistic_body_for_custom_particle():
+    """Test that `RelativisticBody` can be created using a `CustomParticle`."""
+    mass = 1e-27 * u.kg
+    custom_particle = CustomParticle(mass=mass)
+    velocity = 0 * u.m / u.s
+    relativistic_custom_particle = RelativisticBody(
+        particle=custom_particle, V=velocity
+    )
+    assert u.isclose(relativistic_custom_particle.lorentz_factor, 1, rtol=1e-9)
+
+
+def test_relativistic_body_with_particle_list():
+    """
+    Test that `RelativisticBody` can be instantiated with multiple
+    particles.
+    """
+    particles = ["p+", "e-"]
+    relativistic_particles = RelativisticBody(particle=particles, V=0 * u.m / u.s)
+    np.testing.assert_allclose(relativistic_particles.lorentz_factor, 1)
+
+
+def test_relativistic_body_with_multiple_velocities():
+    """
+    Test that `RelativisticBody` can be instantiated with an array of
+    velocities.
+    """
+    velocities = np.array([0, 0.5]) * c
+    relativistic_particles = RelativisticBody("p+", V=velocities)
+    np.testing.assert_allclose(
+        relativistic_particles.lorentz_factor, [1, 1.1547005383792517]
+    )
+
+
+def test_relativistic_body_with_multiple_particles_and_velocities():
+    """
+    Test that `RelativisticBody` can be instantiated with multiple
+    particles and multiple velocities.
+    """
+    velocities = np.array([0.5, 0.7]) * c
+    particles = ["p+", "e-"]
+    relativistic_particles = RelativisticBody(particle=particles, V=velocities)
+    np.testing.assert_allclose(relativistic_particles.velocity, velocities)
+
+
+@pytest.mark.parametrize("function", [repr, str])
+def test_relativistic_body_into_string(function):
+    """Test that `repr` and `str` work on RelativisticBody."""
+    relativistic_body = RelativisticBody("p+", V=5.0 * u.m / u.s)
+    expected = "RelativisticBody(p+, 5.0 m / s)"
+    actual = function(relativistic_body)
+    assert actual == expected
diff --git a/plasmapy/particles/tests/test_factory.py b/plasmapy/particles/tests/test_factory.py
index 34bd0ea9..e3ac0fc6 100644
--- a/plasmapy/particles/tests/test_factory.py
+++ b/plasmapy/particles/tests/test_factory.py
@@ -9,7 +9,7 @@
 
 mass = 1e-26 * u.kg
 charge = 1e-29 * u.C
-custom_particle = CustomParticle(mass, charge)
+custom_particle = CustomParticle(mass=mass, charge=charge)
 
 test_cases = [
     ([[]], {}, ParticleList()),
@@ -18,7 +18,7 @@
     (["H"], {"Z": 1, "mass_numb": 2}, deuteron),
     (["muon"], {}, Particle("muon")),
     pytest.param([charge, mass], {}, custom_particle, marks=[pytest.mark.xfail]),
-    ([mass, charge], {}, custom_particle),
+    ([mass, charge], {"Z": None, "mass_numb": None}, custom_particle),
     ([], {"symbol": "ξ"}, CustomParticle(symbol="ξ")),
     ([[proton, electron]], {}, ParticleList([proton, electron])),
     ([], {"mass": mass}, CustomParticle(mass=mass)),
@@ -28,6 +28,8 @@
     (["H"], {"Z": 1, "mass_numb": 1}, Particle("H", Z=1, mass_numb=1)),
     ([custom_particle], {}, custom_particle),
     ([ParticleList(["p+", "e-"])], {}, ParticleList(["p+", "e-"])),
+    ([mass, charge], {}, custom_particle),
+    ([mass], {}, CustomParticle(mass=mass)),
 ]
 
 
@@ -49,3 +51,14 @@ def test_physical_particle_factory(args, kwargs, expected):
 def test_particle_factory_exceptions(args, kwargs, expected):
     with pytest.raises(expected):
         _physical_particle_factory(*args, **kwargs)
+
+
+def test_particle_factory_custom_particle_with_none_kwargs():
+    """
+    Test that when `_physical_particle_factory` is provided with a
+    `CustomParticle` along with ``Z=None`` and ``mass_numb=None`` as
+    keyword arguments, then it will return the `CustomParticle`.
+    """
+    expected = CustomParticle(mass=1.27 * u.kg, charge=1 * u.C)
+    actual = _physical_particle_factory(expected, Z=None, mass_numb=None)
+    assert expected is actual