diff --git a/hamilflow/__init__.py b/hamilflow/__init__.py
index e69de29..1ec6289 100644
--- a/hamilflow/__init__.py
+++ b/hamilflow/__init__.py
@@ -0,0 +1 @@
+"""Hamilflow package."""
diff --git a/hamilflow/maths/__init__.py b/hamilflow/maths/__init__.py
index e69de29..9fe0c08 100644
--- a/hamilflow/maths/__init__.py
+++ b/hamilflow/maths/__init__.py
@@ -0,0 +1,5 @@
+"""Package for mathematical functions."""
+
+from .trigonometrics import acos_with_shift
+
+__all__ = ["acos_with_shift"]
diff --git a/hamilflow/maths/trigonometrics.py b/hamilflow/maths/trigonometrics.py
index 7e7a328..cddef30 100644
--- a/hamilflow/maths/trigonometrics.py
+++ b/hamilflow/maths/trigonometrics.py
@@ -1,3 +1,5 @@
+"""Trigonometric functions."""
+
 from typing import TYPE_CHECKING
 
 import numpy as np
@@ -12,6 +14,7 @@ def acos_with_shift(
     x: "Collection[float] | npt.ArrayLike",
     shift: "Collection[float] | npt.ArrayLike | None" = None,
 ) -> "npt.ArrayLike":
+    """Arccos with shift."""
     x = np.array(x, copy=False)
     value = np.arccos(x)
     shift = np.array(shift, copy=False)
diff --git a/hamilflow/models/__init__.py b/hamilflow/models/__init__.py
index e69de29..b6d452b 100644
--- a/hamilflow/models/__init__.py
+++ b/hamilflow/models/__init__.py
@@ -0,0 +1 @@
+"""Package for the models."""
diff --git a/hamilflow/models/brownian_motion.py b/hamilflow/models/brownian_motion.py
index e1bff11..2ef8309 100644
--- a/hamilflow/models/brownian_motion.py
+++ b/hamilflow/models/brownian_motion.py
@@ -1,3 +1,5 @@
+"""Main module for Brownian motion."""
+
 from functools import cached_property
 from typing import Mapping, Sequence
 
@@ -11,7 +13,7 @@
 
 
 class BrownianMotionSystem(BaseModel):
-    r"""Definition of the Brownian Motion system
+    r"""Definition of the Brownian Motion system.
 
     For consistency, we always use
     $\mathbf x$ for displacement, and
@@ -24,8 +26,8 @@ class BrownianMotionSystem(BaseModel):
     \end{align}
     $$
 
-    References:
-
+    References
+    ----------
     1. Brownian motion and random walks. [cited 13 Mar 2024].
         Available: https://web.mit.edu/8.334/www/grades/projects/projects17/OscarMickelin/brownian.html
     2. Contributors to Wikimedia projects. Brownian motion.
@@ -35,6 +37,7 @@ class BrownianMotionSystem(BaseModel):
     :cvar sigma: base standard deviation
         to be used to compute the variance
     :cvar delta_t: time granunality of the motion
+
     """
 
     sigma: float = Field(ge=0.0)
@@ -43,14 +46,12 @@ class BrownianMotionSystem(BaseModel):
     @computed_field  # type: ignore[misc]
     @cached_property
     def gaussian_scale(self) -> float:
-        """The scale (standard deviation) of the Gaussian term
-        in Brownian motion
-        """
+        """The scale (standard deviation) of the Gaussian term in Brownian motion."""
         return self.sigma**2 * self.delta_t
 
 
 class BrownianMotionIC(BaseModel):
-    """The initial condition for a Brownian motion
+    """The initial condition for a Brownian motion.
 
     :cvar x0: initial displacement of the particle,
         the diminsion of this initial condition determines
@@ -61,7 +62,7 @@ class BrownianMotionIC(BaseModel):
 
     @field_validator("x0")
     @classmethod
-    def check_x0_types(cls, v: float | Sequence[float]) -> float | Sequence[float]:
+    def _check_x0_types(cls, v: float | Sequence[float]) -> float | Sequence[float]:
         if not isinstance(v, (float, int, Sequence)):
             # TODO I do not think this raise can be reached
             raise ValueError(f"Value of x0 should be int/float/list of int/float: {v=}")
@@ -70,8 +71,8 @@ def check_x0_types(cls, v: float | Sequence[float]) -> float | Sequence[float]:
 
 
 class BrownianMotion:
-    r"""Brownian motion describes motion of small particles
-    with stochastic forces applied to them.
+    r"""Brownian motion describes motion of small particles with stochastic forces applied to them.
+
     The math of Brownian motion can be modeled
     with Wiener process.
 
@@ -86,8 +87,8 @@ class BrownianMotion:
     \end{align}
     $$
 
-    References:
-
+    References
+    ----------
     1. Brownian motion and random walks. [cited 13 Mar 2024].
         Available: https://web.mit.edu/8.334/www/grades/projects/projects17/OscarMickelin/brownian.html
     2. Contributors to Wikimedia projects. Brownian motion.
@@ -142,6 +143,7 @@ class BrownianMotion:
 
     :param system: the Brownian motion system definition
     :param initial_condition: the initial condition for the simulation
+
     """
 
     def __init__(
@@ -157,7 +159,7 @@ def __init__(
 
     @property
     def dim(self) -> int:
-        """Dimension of the Brownian motion"""
+        """Dimension of the Brownian motion."""
         return np.array(self.initial_condition.x0, copy=False).size
 
     @property
@@ -165,7 +167,7 @@ def _axis_names(self) -> list[str]:
         return [f"x_{i}" for i in range(self.dim)]
 
     def _trajectory(self, n_new_steps: int, seed: int) -> "npt.NDArray[np.float64]":
-        """The trajectory of the particle.
+        """Give the trajectory of the particle.
 
         We first compute the delta displacement in each step.
         With the displacement at each step, we perform a cumsum
@@ -189,7 +191,7 @@ def _trajectory(self, n_new_steps: int, seed: int) -> "npt.NDArray[np.float64]":
         return trajectory
 
     def generate_from(self, n_steps: int, seed: int = 42) -> pd.DataFrame:
-        """generate data from a set of interpretable params for this model
+        """Generate data from a set of interpretable params for this model.
 
         :param n_steps: total number of steps to be simulated, including the inital step.
         :param seed: random generator seed for the stochastic process.
@@ -200,7 +202,7 @@ def generate_from(self, n_steps: int, seed: int = 42) -> pd.DataFrame:
         return self(t=time_steps, seed=seed)
 
     def __call__(self, t: TypeTime, seed: int = 42) -> pd.DataFrame:
-        """Simulate the coordinates of the particle
+        """Simulate the coordinates of the particle.
 
         :param t: the time sequence to be used to generate data, 1-D array like
         :param seed: random generator seed for the stochastic process.
diff --git a/hamilflow/models/free_particle.py b/hamilflow/models/free_particle.py
index c36bfb5..57c25ab 100644
--- a/hamilflow/models/free_particle.py
+++ b/hamilflow/models/free_particle.py
@@ -1,3 +1,5 @@
+"""Main module for a free particle."""
+
 from functools import cached_property
 from typing import Mapping, Sequence, cast
 
@@ -13,7 +15,7 @@
 
 
 class FreeParticleIC(BaseModel):
-    """The initial condition for a free particle
+    """The initial condition for a free particle.
 
     :cvar x0: the initial displacement
     :cvar v0: the initial velocity
@@ -23,7 +25,7 @@ class FreeParticleIC(BaseModel):
     v0: float | Sequence[float] = Field()
 
     @model_validator(mode="after")
-    def check_dimensions_match(self) -> Self:
+    def _check_dimensions_match(self) -> Self:
         if (x0_seq := isinstance(self.x0, Sequence)) != isinstance(self.v0, Sequence):
             raise TypeError("x0 and v0 need both to be scalars or Sequences")
         elif x0_seq and len(cast(Sequence, self.x0)) != len(cast(Sequence, self.v0)):
@@ -46,7 +48,7 @@ def __init__(
 
     @cached_property
     def definition(self) -> dict[str, dict[str, float | list[float]]]:
-        """model params and initial conditions defined as a dictionary."""
+        """Model params and initial conditions defined as a dictionary."""
         return dict(initial_condition=self.initial_condition.model_dump())
 
     def _x(self, t: "Sequence[float] |  npt.ArrayLike") -> "npt.NDArray[np.float64]":
@@ -60,7 +62,6 @@ def __call__(self, t: "Sequence[float] |  npt.ArrayLike") -> pd.DataFrame:
 
         :param t: time(s).
         """
-
         data = self._x(t)
         columns = [f"x{i+1}" for i in range(data.shape[1])]
 
diff --git a/hamilflow/models/harmonic_oscillator.py b/hamilflow/models/harmonic_oscillator.py
index 2199e82..b3454b5 100644
--- a/hamilflow/models/harmonic_oscillator.py
+++ b/hamilflow/models/harmonic_oscillator.py
@@ -1,3 +1,5 @@
+"""Main module for undamped and damped hamornic oscillators."""
+
 from abc import ABC, abstractmethod
 from functools import cached_property
 from typing import Literal, Mapping, Sequence
@@ -9,7 +11,7 @@
 
 
 class HarmonicOscillatorSystem(BaseModel):
-    """The params for the harmonic oscillator
+    """The params for the harmonic oscillator.
 
     :cvar omega: angular frequency of the harmonic oscillator
     :cvar zeta: damping ratio
@@ -21,13 +23,13 @@ class HarmonicOscillatorSystem(BaseModel):
     @computed_field  # type: ignore[misc]
     @cached_property
     def period(self) -> float:
-        """period of the oscillator"""
+        """Period of the oscillator."""
         return 2 * np.pi / self.omega
 
     @computed_field  # type: ignore[misc]
     @cached_property
     def frequency(self) -> float:
-        """frequency of the oscillator"""
+        """Frequency of the oscillator."""
         return 1 / self.period
 
     @computed_field  # type: ignore[misc]
@@ -35,7 +37,7 @@ def frequency(self) -> float:
     def type(
         self,
     ) -> Literal["simple", "under_damped", "critical_damped", "over_damped"]:
-        """which type of harmonic oscillators"""
+        """Which type of harmonic oscillators."""
         if self.zeta == 0:
             return "simple"
         elif self.zeta < 1:
@@ -47,7 +49,7 @@ def type(
 
     @field_validator("zeta")
     @classmethod
-    def check_zeta_non_negative(cls, v: float) -> float:
+    def _check_zeta_non_negative(cls, v: float) -> float:
         if v < 0:
             raise ValueError(f"Value of zeta should be positive: {v=}")
 
@@ -55,7 +57,7 @@ def check_zeta_non_negative(cls, v: float) -> float:
 
 
 class HarmonicOscillatorIC(BaseModel):
-    """The initial condition for a harmonic oscillator
+    """The initial condition for a harmonic oscillator.
 
     :cvar x0: the initial displacement
     :cvar v0: the initial velocity
@@ -68,8 +70,7 @@ class HarmonicOscillatorIC(BaseModel):
 
 
 class HarmonicOscillatorBase(ABC):
-    r"""Base class to generate time series data
-    for a [harmonic oscillator](https://en.wikipedia.org/wiki/Harmonic_oscillator).
+    r"""Base class to generate time series data for a [harmonic oscillator](https://en.wikipedia.org/wiki/Harmonic_oscillator).
 
     :param system: all the params that defines the harmonic oscillator.
     :param initial_condition: the initial condition of the harmonic oscillator.
@@ -86,7 +87,7 @@ def __init__(
 
     @cached_property
     def definition(self) -> dict[str, dict[str, float]]:
-        """model params and initial conditions defined as a dictionary."""
+        """Model params and initial conditions defined as a dictionary."""
         return {
             "system": self.system.model_dump(),
             "initial_condition": self.initial_condition.model_dump(),
@@ -114,9 +115,7 @@ def __call__(self, n_periods: int, n_samples_per_period: int) -> pd.DataFrame:
 
 
 class SimpleHarmonicOscillator(HarmonicOscillatorBase):
-    r"""Generate time series data for a
-    [simple harmonic oscillator](https://en.wikipedia.org/wiki/Harmonic_oscillator).
-
+    r"""Generate time series data for a [simple harmonic oscillator](https://en.wikipedia.org/wiki/Harmonic_oscillator).
 
     In a one dimensional world, a mass $m$, driven by a force $F=-kx$, is described as
 
@@ -163,7 +162,7 @@ def __init__(
             )
 
     def _x(self, t: "Sequence[float] | npt.ArrayLike") -> np.ndarray:
-        r"""Solution to simple harmonic oscillators:
+        r"""Solution to simple harmonic oscillators.
 
         $$
         x(t) = x_0 \cos(\omega t + \phi).
@@ -236,7 +235,7 @@ def __init__(
             )
 
     def _x_under_damped(self, t: "Sequence[float] | npt.ArrayLike") -> npt.ArrayLike:
-        r"""Solution to under damped harmonic oscillators:
+        r"""Solution to under damped harmonic oscillators.
 
         $$
         x(t) = \left( x_0 \cos(\Omega t) + \frac{\zeta \omega x_0 + v_0}{\Omega} \sin(\Omega t) \right)
@@ -262,7 +261,7 @@ def _x_under_damped(self, t: "Sequence[float] | npt.ArrayLike") -> npt.ArrayLike
         ) * np.exp(-self.system.zeta * self.system.omega * t)
 
     def _x_critical_damped(self, t: "Sequence[float] | npt.ArrayLike") -> npt.ArrayLike:
-        r"""Solution to critical damped harmonic oscillators:
+        r"""Solution to critical damped harmonic oscillators.
 
         $$
         x(t) = \left( x_0 \cos(\Omega t) + \frac{\zeta \omega x_0 + v_0}{\Omega} \sin(\Omega t) \right)
@@ -281,7 +280,7 @@ def _x_critical_damped(self, t: "Sequence[float] | npt.ArrayLike") -> npt.ArrayL
         )
 
     def _x_over_damped(self, t: "Sequence[float] | npt.ArrayLike") -> npt.ArrayLike:
-        r"""Solution to over harmonic oscillators:
+        r"""Solution to over harmonic oscillators.
 
         $$
         x(t) = \left( x_0 \cosh(\Gamma t) + \frac{\zeta \omega x_0 + v_0}{\Gamma} \sinh(\Gamma t) \right)
@@ -325,7 +324,7 @@ def _x(self, t: "Sequence[float] | npt.ArrayLike") -> npt.ArrayLike:
 
 
 class ComplexSimpleHarmonicOscillatorIC(BaseModel):
-    """The initial condition for a complex harmonic oscillator
+    """The initial condition for a complex harmonic oscillator.
 
     :cvar x0: the initial displacements
     :cvar phi: initial phases
@@ -360,15 +359,14 @@ def __init__(
     def definition(
         self,
     ) -> dict[str, dict[str, float | tuple[float, float]]]:
-        """model params and initial conditions defined as a dictionary."""
-
+        """Model params and initial conditions defined as a dictionary."""
         return dict(
             system=self.system.model_dump(),
             initial_condition=self.initial_condition.model_dump(),
         )
 
     def _z(self, t: "Sequence[float] | npt.ArrayLike") -> npt.ArrayLike:
-        r"""Solution to complex simple harmonic oscillators:
+        r"""Solution to complex simple harmonic oscillators.
 
         $$
         x(t) = x_+ \exp(-\mathbb{i} (\omega t + \phi_+)) + x_- \exp(+\mathbb{i} (\omega t + \phi_-)).
diff --git a/hamilflow/models/harmonic_oscillator_chain.py b/hamilflow/models/harmonic_oscillator_chain.py
index e0bbc19..d6236a4 100644
--- a/hamilflow/models/harmonic_oscillator_chain.py
+++ b/hamilflow/models/harmonic_oscillator_chain.py
@@ -1,3 +1,5 @@
+"""Main module for a harmonic oscillator chain."""
+
 from functools import cached_property
 from typing import Mapping, Sequence, cast
 
@@ -11,8 +13,7 @@
 
 
 class HarmonicOscillatorsChain:
-    r"""Generate time series data for a coupled harmonic oscillator chain
-    with periodic boundary condition.
+    r"""Generate time series data for a coupled harmonic oscillator chain with periodic boundary condition.
 
     A one-dimensional circle of $N$ interacting harmonic oscillators can be described by the Lagrangian action
     $$S_L[x_i] = \int_{t_0}^{t_1}\mathbb{d} t \left\{ \sum_{i=0}^{N-1} \frac{1}{2}m \dot x_i^2 - \frac{1}{2}m\omega^2\left(x_i - x_{i+1}\right)^2 \right\}\,,$$
@@ -81,7 +82,7 @@ def definition(
         | dict[str, dict[str, float | list[float]]]
         | list[dict[str, dict[str, float | tuple[float, float]]]],
     ]:
-        """model params and initial conditions defined as a dictionary."""
+        """Model params and initial conditions defined as a dictionary."""
         return dict(
             omega=self.omega,
             n_dof=self.n_dof,
diff --git a/hamilflow/models/kepler_problem/__init__.py b/hamilflow/models/kepler_problem/__init__.py
index c669344..5896a48 100644
--- a/hamilflow/models/kepler_problem/__init__.py
+++ b/hamilflow/models/kepler_problem/__init__.py
@@ -1,3 +1,5 @@
+"""Package for the Kepler problem."""
+
 from .model import Kepler2D, Kepler2DIoM, Kepler2DSystem
 
 __all__ = ["Kepler2D", "Kepler2DIoM", "Kepler2DSystem"]
diff --git a/hamilflow/models/kepler_problem/dynamics.py b/hamilflow/models/kepler_problem/dynamics.py
index 165bc84..c978e45 100644
--- a/hamilflow/models/kepler_problem/dynamics.py
+++ b/hamilflow/models/kepler_problem/dynamics.py
@@ -1,3 +1,5 @@
+"""Exact solution of Kepler dynamics."""
+
 from typing import TYPE_CHECKING
 
 import numpy as np
@@ -37,6 +39,12 @@ def _tau_of_e_minus_u_elliptic(
 
 
 def tau_of_u_elliptic(ecc: float, u: "npt.ArrayLike") -> "npt.NDArray[np.float64]":
+    """Calculate the scaled time tau from u in the elliptic case.
+
+    :param ecc: eccentricity, 0 < ecc < 1 (unchecked)
+    :param u: convenient radial inverse
+    :return: scaled time tau
+    """
     return _approximate_at_termina(
         ecc,
         u,
@@ -47,6 +55,12 @@ def tau_of_u_elliptic(ecc: float, u: "npt.ArrayLike") -> "npt.NDArray[np.float64
 
 
 def tau_of_u_parabolic(ecc: float, u: "npt.ArrayLike") -> "npt.NDArray[np.float64]":
+    """Calculate the scaled time tau from u in the parabolic case.
+
+    :param ecc: eccentricity, ecc == 1 (unchecked, unused)
+    :param u: convenient radial inverse
+    :return: scaled time tau
+    """
     u = np.array(u, copy=False)
     return np.sqrt(1 - u) * (2 + u) / 3 / (1 + u) ** 1.5
 
@@ -81,6 +95,12 @@ def _tau_of_e_minus_u_hyperbolic(
 
 
 def tau_of_u_hyperbolic(ecc: float, u: "npt.ArrayLike") -> "npt.NDArray[np.float64]":
+    """Calculate the scaled time tau from u in the hyperbolic case.
+
+    :param ecc: eccentricity, ecc > 1 (unchecked)
+    :param u: convenient radial inverse
+    :return: scaled time tau
+    """
     return _approximate_at_termina(
         ecc,
         u,
@@ -91,11 +111,23 @@ def tau_of_u_hyperbolic(ecc: float, u: "npt.ArrayLike") -> "npt.NDArray[np.float
 
 
 def tau_of_u_prime(ecc: float, u: "npt.ArrayLike") -> "npt.NDArray[np.float64]":
+    """Calculate the first derivative of scaled time tau with respect to u.
+
+    :param ecc: eccentricity, ecc >= 0 (unchecked)
+    :param u: convenient radial inverse
+    :return: the first derivative scaled time tau with respect to u
+    """
     u = np.array(u, copy=False)
     return -1 / (1 + u) ** 2 / np.sqrt(ecc**2 - u**2)
 
 
 def tau_of_u_prime2(ecc: float, u: "npt.ArrayLike") -> "npt.NDArray[np.float64]":
+    """Calculate the second derivative of scaled time tau with respect to u.
+
+    :param ecc: eccentricity, ecc >= 0 (unchecked)
+    :param u: convenient radial inverse
+    :return: the second derivative scaled time tau with respect to u
+    """
     u = np.array(u, copy=False)
     u2 = u**2
     return (2 * ecc**2 - u - 3 * u2) / (1 + u) ** 3 / (ecc**2 - u2) ** 1.5
@@ -105,6 +137,12 @@ def esolve_u_from_tau_parabolic(
     ecc: float,
     tau: "npt.ArrayLike",
 ) -> "npt.NDArray[np.float64]":
+    """Calculate the convenient radial inverse u from tau in the parabolic case, using the exact solution.
+
+    :param ecc: eccentricity, ecc = 0 (unchecked, unused)
+    :param tau: scaled time
+    :return: convenient radial inverse u
+    """
     tau = np.array(tau, copy=False)
     tau_3 = 3 * tau
     term = 1 + tau_3**2  # 1 + 9 * tau**2
diff --git a/hamilflow/models/kepler_problem/model.py b/hamilflow/models/kepler_problem/model.py
index ae1aadf..af02417 100644
--- a/hamilflow/models/kepler_problem/model.py
+++ b/hamilflow/models/kepler_problem/model.py
@@ -1,3 +1,5 @@
+"""Main module for Kepler problem."""
+
 import math
 from functools import cached_property, partial
 from typing import TYPE_CHECKING, Any
@@ -18,7 +20,7 @@
 
 
 class Kepler2DSystem(BaseModel):
-    r"""Definition of the Kepler problem
+    r"""Definition of the Kepler problem.
 
     Potential:
 
@@ -44,7 +46,7 @@ class Kepler2DSystem(BaseModel):
 
 
 class Kepler2DIoM(BaseModel):
-    """The integrals of motion for a Kepler problem
+    """The integrals of motion for a Kepler problem.
 
     :cvar ene: the energy
     :cvar angular_mom: the angular momentum
@@ -60,17 +62,17 @@ class Kepler2DIoM(BaseModel):
     # TODO process angular momentum = 0
     @field_validator("angular_mom")
     @classmethod
-    def angular_mom_non_zero(cls, v: Any) -> float:
+    def _angular_mom_non_zero(cls, v: Any) -> float:
         if v == 0:
             raise NotImplementedError("Only non-zero angular momenta are supported")
         return v
 
 
 class Kepler2D:
-    r"""Kepler problem in two dimensional space.
+    """Kepler problem in two dimensional space.
 
-    :param system: the Kepler problem system definition
-    :param initial_condition: the initial condition for the simulation
+    :param system: the Kepler problem system specification
+    :param integrals_of_motion: the integrals of motion for the system.
     """
 
     def __init__(
@@ -104,6 +106,14 @@ def from_geometry(
         system: "Mapping[str, float]",
         geometries: "Mapping[str, bool | float]",
     ) -> "Self":
+        """Alternative initialiser from system and geometry specifications.
+
+        :param system: the Kepler problem system specification
+        :param geometries: geometric specifications
+            `positive_angular_mom`: whether the angular momentum is positive
+            `ecc`: eccentricity of the conic section
+            `parameter`: parameter of the conic section
+        """
         mass, alpha = system["mass"], system["alpha"]
         positive_angular_mom = bool(geometries["positive_angular_mom"])
         ecc, parameter = map(lambda k: float(geometries[k]), ["ecc", "parameter"])
@@ -122,35 +132,48 @@ def minimal_ene(
         angular_mom: float,
         system: "Mapping[str, float]",
     ) -> float:
+        """Minimal possible energy from the system specification and an angular momentum.
+
+        :param angular_mom: angular momentum
+        :param system: system specification
+        :return: minimal possible energy
+        """
         mass, alpha = system["mass"], system["alpha"]
         return -mass * alpha**2 / (2 * angular_mom**2)
 
     @property
     def mass(self) -> float:
+        """Mass parameter from the system specification."""
         return self.system.mass
 
     @property
     def alpha(self) -> float:
+        """Alpha parameter from the system specification."""
         return self.system.alpha
 
     @property
     def ene(self) -> float:
+        """Energy of the Kepler problem."""
         return self.integrals_of_motion.ene
 
     @property
     def angular_mom(self) -> float:
+        """Angular momentum of the Kepler problem."""
         return self.integrals_of_motion.angular_mom
 
     @property
     def t0(self) -> float:
+        """t0 of the Kepler problem."""
         return self.integrals_of_motion.t0
 
     @property
     def phi0(self) -> float:
+        """phi0 of the Kepler problem."""
         return self.integrals_of_motion.phi0
 
     @cached_property
     def period(self) -> float:
+        """Perior of the Kepler problem."""
         if self.ene >= 0:
             msg = f"Only systems with energy < 0 have a period, got {self.ene}"
             raise TypeError(msg)
@@ -159,16 +182,19 @@ def period(self) -> float:
     # FIXME is it called parameter in English?
     @cached_property
     def parameter(self) -> float:
+        """Conic section parameter of the Kepler problem."""
         return self.angular_mom**2 / self.mass / self.alpha
 
     @cached_property
     def ecc(self) -> float:
+        """Sonic section eccentricity of the Kepler problem."""
         return math.sqrt(
             1 + 2 * self.ene * self.angular_mom**2 / self.mass / self.alpha**2,
         )
 
     @cached_property
     def period_in_tau(self) -> float:
+        """Period in the scaled time tau."""
         if self.ecc >= 1:
             raise TypeError(
                 f"Only systems with 0 <= eccentricity < 1 have a period, got {self.ecc}",
@@ -177,12 +203,15 @@ def period_in_tau(self) -> float:
 
     @property
     def t_to_tau_factor(self) -> float:
+        """Scale factor from t to tau."""
         return abs(self.mass * self.alpha**2 / self.angular_mom**3)
 
     def tau(self, t: "Collection[float] | npt.ArrayLike") -> "npt.ArrayLike":
+        """Give the scaled time tau from t."""
         return (np.array(t, copy=False) - self.t0) * self.t_to_tau_factor
 
     def u_of_tau(self, tau: "Collection[float] | npt.ArrayLike") -> "npt.ArrayLike":
+        """Give the convenient radial inverse u from tau."""
         tau = np.array(tau, copy=False)
         if self.ecc == 0:
             return np.zeros(tau.shape)
@@ -196,6 +225,7 @@ def u_of_tau(self, tau: "Collection[float] | npt.ArrayLike") -> "npt.ArrayLike":
             return u_of_tau(self.ecc, tau)  # type: ignore [arg-type]
 
     def r_of_u(self, u: "Collection[float] | npt.ArrayLike") -> "npt.ArrayLike":
+        """Give the radial r from u."""
         return self.parameter / (np.array(u, copy=False) + 1)
 
     def phi_of_u_tau(
@@ -203,6 +233,7 @@ def phi_of_u_tau(
         u: "Collection[float] | npt.ArrayLike",
         tau: "Collection[float] | npt.ArrayLike",
     ) -> "npt.ArrayLike":
+        """Give the angular phi from u and tau."""
         u, tau = np.array(u, copy=False), np.array(tau, copy=False)
         if self.ecc == 0:
             phi = 2 * math.pi * tau / self.period_in_tau
@@ -215,6 +246,7 @@ def phi_of_u_tau(
         return phi + self.phi0
 
     def __call__(self, t: "Collection[float] | npt.ArrayLike") -> pd.DataFrame:
+        """Give a DataFrame of tau, u, r and phi from t."""
         tau = self.tau(t)
         u = self.u_of_tau(tau)
         r = self.r_of_u(u)
diff --git a/hamilflow/models/kepler_problem/numerics.py b/hamilflow/models/kepler_problem/numerics.py
index 015572f..87ffccd 100644
--- a/hamilflow/models/kepler_problem/numerics.py
+++ b/hamilflow/models/kepler_problem/numerics.py
@@ -1,3 +1,5 @@
+"""Numerics for the Kepler problem."""
+
 from typing import TYPE_CHECKING, Literal
 
 import numpy as np
@@ -34,6 +36,13 @@ def _u0_hyperbolic(
 
 
 def nsolve_u_from_tau_newton(ecc: float, tau: "npt.ArrayLike") -> OptimizeResult:
+    """Calculate the convenient radial inverse u from tau in the elliptic or parabolic case, using the Newton method.
+
+    :param ecc: eccentricity, ecc > 0, ecc != 1
+    :param tau: scaled time
+    :raises ValueError: when `ecc` is invalid
+    :return: numeric OptimizeResult from scipy
+    """
     tau = np.array(tau, copy=False)
     if 0 < ecc < 1:
         tau_of_u = tau_of_u_elliptic
@@ -57,6 +66,12 @@ def fprime2(u: float, tau: "npt.ArrayLike") -> "npt.NDArray[np.float64]":
 
 
 def nsolve_u_from_tau_bisect(ecc: float, tau: "npt.ArrayLike") -> list[OptimizeResult]:
+    """Calculate the convenient radial inverse u from tau in the elliptic or parabolic case, using the bisect method.
+
+    :param ecc: eccentricity, ecc > 0, ecc != 1
+    :param tau: scaled time
+    :return: numeric OptimizeResult from scipy
+    """
     tau_s = np.array(tau, copy=False).reshape(-1)
     if 0 < ecc < 1:
         tau_of_u = tau_of_u_elliptic
@@ -78,6 +93,15 @@ def u_of_tau(
     tau: "npt.ArrayLike",
     method: Literal["bisect", "newton"] = "bisect",
 ) -> "npt.NDArray[np.float64]":
+    """Calculate the convenient radial inverse u from tau, using numeric methods.
+
+    :param ecc: eccentricity, ecc >= 0
+    :param tau: scaled time
+    :param method: "newton" or "bisect" numeric methods, defaults to "bisect"
+    :raises ValueError: when `method` is invalid
+    :raises ValueError: when `ecc` is invalid
+    :return: convenient radial inverse u
+    """
     tau = np.array(tau, copy=False)
     if ecc == 0:
         return np.zeros(tau.shape)
diff --git a/hamilflow/models/pendulum.py b/hamilflow/models/pendulum.py
index 22a6864..60b9831 100644
--- a/hamilflow/models/pendulum.py
+++ b/hamilflow/models/pendulum.py
@@ -1,3 +1,5 @@
+"""Main module for a pendulum."""
+
 import math
 from functools import cached_property
 from typing import Any, Mapping, Sequence
@@ -77,6 +79,7 @@ def definition(self) -> dict[str, dict[str, Any]]:
 
     @property
     def omega0(self) -> float:
+        """Original angular frequency of the system."""
         return self.system.omega0
 
     @property
@@ -112,8 +115,7 @@ def _math_u(
         return self.omega0 * np.array(t, copy=False)
 
     def u(self, t: "Sequence[float] | npt.ArrayLike") -> "npt.NDArray[np.float64]":
-        r"""The convenient generalised coordinate $u$,
-        $\sin u \coloneqq \frac{\sin\frac{\theta}{2}}{\sin\frac{\theta_0}{2}}$.
+        r"""Give the convenient generalised coordinate $u$, $\sin u \coloneqq \frac{\sin\frac{\theta}{2}}{\sin\frac{\theta_0}{2}}$.
 
         :param t: time
         :return: $u(t) = \mathrm{am}\!\big(\omega_0 t + K(k^2), k^2\big)$, where
@@ -136,7 +138,7 @@ def theta(self, t: "Sequence[float] | npt.ArrayLike") -> "npt.NDArray[np.float64
         return 2 * np.arcsin(cn / dn * self._k)
 
     def generate_from(self, n_periods: int, n_samples_per_period: int) -> pd.DataFrame:
-        """generate the time sequence from more interpretable params
+        """Generate the time sequence from more interpretable params.
 
         :param n_periods: number of periods to include
         :param n_samples_per_period: number of samples in each period
@@ -148,7 +150,7 @@ def generate_from(self, n_periods: int, n_samples_per_period: int) -> pd.DataFra
         return self(time_steps)
 
     def __call__(self, t: TypeTime) -> pd.DataFrame:
-        """generate the variables of the pendulum in time together with the time steps
+        """Generate the variables of the pendulum in time together with the time steps.
 
         :param t: time steps
         :return: values of the variables
diff --git a/hamilflow/models/utils/__init__.py b/hamilflow/models/utils/__init__.py
index e69de29..10bcf58 100644
--- a/hamilflow/models/utils/__init__.py
+++ b/hamilflow/models/utils/__init__.py
@@ -0,0 +1 @@
+"""Package for utilities."""
diff --git a/hamilflow/models/utils/typing.py b/hamilflow/models/utils/typing.py
index 2c2da43..73daf87 100644
--- a/hamilflow/models/utils/typing.py
+++ b/hamilflow/models/utils/typing.py
@@ -1,3 +1,5 @@
+"""Typing."""
+
 from typing import Sequence, TypeVar
 
 from numpy.typing import ArrayLike
diff --git a/pyproject.toml b/pyproject.toml
index 3bb7f88..d9ccdfd 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -72,7 +72,16 @@ line-length = 88  # black default
 
 [tool.ruff.lint]
 select = ["ALL"]
-ignore = ["D", "FA", "PLR", "S", "PD", "SLF", "PT", "FBT", "C", "ANN", "EM", "TRY", "ARG", "FIX", "TD", "TID", "ERA", "INP", "E", "SIM", "RUF", "RET"]
+ignore = [
+    "D107", # undocumented-public-init: we document the class instead
+    "FA", "PLR", "S", "PD", "SLF", "PT", "FBT", "C", "ANN", "EM", "TRY", "ARG", "FIX", "TD", "TID", "ERA", "INP", "E", "SIM", "RUF", "RET"
+]
+
+
+[tool.ruff.lint.per-file-ignores]
+"docs/**/*.py" = [
+    "D100",  # undocumented-public-module: we keep tutorials clean and do not include doc strings
+]
 
 
 [tool.ruff.lint.isort]
diff --git a/tests/maths/__init__.py b/tests/maths/__init__.py
index e69de29..a3e5d76 100644
--- a/tests/maths/__init__.py
+++ b/tests/maths/__init__.py
@@ -0,0 +1 @@
+"""Package for tests for the math modules."""
diff --git a/tests/maths/test_trigonometrics.py b/tests/maths/test_trigonometrics.py
index eef704d..a1aa31c 100644
--- a/tests/maths/test_trigonometrics.py
+++ b/tests/maths/test_trigonometrics.py
@@ -1,3 +1,5 @@
+"""Tests for the trigonometrics module."""
+
 from typing import TYPE_CHECKING, ClassVar
 
 import numpy as np
@@ -11,13 +13,17 @@
 
 
 class TestAcosWithShift:
+    """Tests for arccos with shift."""
+
     _some_numbers: ClassVar[list[float]] = [x / 2 - 5 for x in range(20)]
 
     @pytest.fixture(params=[_some_numbers[0], _some_numbers])
     def phi(self, request: pytest.FixtureRequest) -> float | list[float]:
+        """Give scalar phi and a list of phi's."""
         return request.param
 
     def test_acos_with_shift(self, phi: "float | Collection[float]") -> None:
+        """Test arccos with shift."""
         phi = np.array(phi, copy=False)
         actual = np.array(acos_with_shift(np.cos(phi), phi / 2 / np.pi), copy=False)
         assert_array_almost_equal(phi, actual)
diff --git a/tests/models/__init__.py b/tests/models/__init__.py
index e69de29..541e1b4 100644
--- a/tests/models/__init__.py
+++ b/tests/models/__init__.py
@@ -0,0 +1 @@
+"""Package for tests for models."""
diff --git a/tests/models/kepler_problem/__init__.py b/tests/models/kepler_problem/__init__.py
index e69de29..44db054 100644
--- a/tests/models/kepler_problem/__init__.py
+++ b/tests/models/kepler_problem/__init__.py
@@ -0,0 +1 @@
+"""Package for tests for the Kepler problem module."""
diff --git a/tests/models/kepler_problem/conftest.py b/tests/models/kepler_problem/conftest.py
index 3071726..bb967a7 100644
--- a/tests/models/kepler_problem/conftest.py
+++ b/tests/models/kepler_problem/conftest.py
@@ -1,3 +1,5 @@
+"""Share fixtures across multiple files."""
+
 from typing import TYPE_CHECKING
 
 import numpy as np
@@ -23,11 +25,13 @@
     params=[0.0, 0.1, 0.3, 0.7, 0.9, 1.0, 1.1, 2.0, 11.0, 101.0],
 )
 def ecc(request: pytest.FixtureRequest) -> float:
+    """Give a few eccentricities."""
     return request.param
 
 
 @pytest.fixture(scope="module", params=[-0.9, False, 0.9])
 def u_s(request: pytest.FixtureRequest, ecc: float) -> "npt.ArrayLike":
+    """Give scalar u or a numpy array of u's."""
     # There are dividends sqrt(e**2 - u**2) and (u + 1),
     # hence u cannot be too close to +e / -e / -1
     f = 1 - _EPS_ECC
@@ -37,6 +41,12 @@ def u_s(request: pytest.FixtureRequest, ecc: float) -> "npt.ArrayLike":
 
 @pytest.fixture(scope="module")
 def tau_of_u(ecc: float) -> "Callable[[float, npt.ArrayLike], npt.ArrayLike]":
+    """Give function calculating tau from u.
+
+    :param ecc: eccentricity, ecc >= 0
+    :raises ValueError: when ecc is invalid
+    :return: function calculating tau from u
+    """
     if ecc == 0:
         pytest.skip("Circular case")
     elif 0 < ecc < 1:
diff --git a/tests/models/kepler_problem/test_dynamics.py b/tests/models/kepler_problem/test_dynamics.py
index 45d270c..a14a2f3 100644
--- a/tests/models/kepler_problem/test_dynamics.py
+++ b/tests/models/kepler_problem/test_dynamics.py
@@ -1,3 +1,5 @@
+"""Tests for the Kepler dynamics module."""
+
 from typing import TYPE_CHECKING
 
 import numpy as np
@@ -29,13 +31,18 @@
 
 @pytest.mark.usefixtures("ecc", "u_s", "tau_of_u")
 class TestTauOfU:
+    """Tests for calculating tau from u."""
+
     def test_const(
         self,
         ecc: float,
         tau_of_u: "Callable[[float, npt.ArrayLike], npt.ArrayLike]",
     ) -> None:
-        # There are dividends sqrt(e**2 - u**2) and (u + 1),
-        # hence u cannot be too close to +e / -e / -1
+        """Test that tau at u = ecc should be 0.
+
+        There are dividends sqrt(e**2 - u**2) and (u + 1),
+        hence u cannot be too close to +e / -e / -1.
+        """
         res = np.array(tau_of_u(ecc, ecc * (1 - _EPS_SQRT)), copy=False)
         desired = np.full(res.shape, _EPS_SQRT)
         assert_array_almost_equal(desired, res)
@@ -46,6 +53,11 @@ def test_tau_of_u(
         u_s: "npt.ArrayLike",
         tau_of_u: "Callable[[float, npt.ArrayLike], npt.ArrayLike]",
     ) -> None:
+        """Test that tau from u is consistent with the integral from its first derivative.
+
+        The integral constant has been tested in another test.
+        """
+
         def integrand(u: "npt.ArrayLike") -> "npt.ArrayLike":
             return tau_of_u_prime(ecc, u)
 
@@ -59,6 +71,11 @@ def exact_and_approx_tau_s(
         self,
         ecc: float,
     ) -> "tuple[Callable[[float, npt.NDArray[np.float64]], npt.NDArray[np.float64]], Callable[[float, npt.NDArray[np.float64]], npt.NDArray[np.float64]], Callable[[float, npt.NDArray[np.float64]], npt.NDArray[np.float64]]]":
+        """Give approximate and exact solutions for the elliptic and hyperbolic cases.
+
+        The exact solutions have removable singularities at the boundary of domain,
+        hence approximate solutions are needed.
+        """
         if ecc == 1 or ecc == 0:
             c = "Parabolic" if ecc else "Circular"
             pytest.skip(f"{c} case is exact")
@@ -84,6 +101,11 @@ def test_expansion(
         epsilon: float,
         exact_and_approx_tau_s: "tuple[Callable[[float, npt.NDArray[np.float64]], npt.NDArray[np.float64]], Callable[[float, npt.NDArray[np.float64]], npt.NDArray[np.float64]], Callable[[float, npt.NDArray[np.float64]], npt.NDArray[np.float64]]]",
     ) -> None:
+        """Test that the approximate solutions near the booundary of domain are close to the exact solution.
+
+        The exact solutions have removable singularities at the boundary of domain,
+        hence approximate solutions are needed.
+        """
         if ecc == 0.0 or ecc == 1.0:
             pytest.skip(f"Test applies to ecc > 0, ecc != 1, got {ecc}")
         factor = 1 - epsilon
diff --git a/tests/models/kepler_problem/test_model.py b/tests/models/kepler_problem/test_model.py
index e2707ca..7d5095a 100644
--- a/tests/models/kepler_problem/test_model.py
+++ b/tests/models/kepler_problem/test_model.py
@@ -1,3 +1,5 @@
+"""Tests for the Kepler model classes."""
+
 import math
 from typing import TYPE_CHECKING, ClassVar
 
@@ -20,21 +22,25 @@
 
 @pytest.fixture(params=[(1.0, 1.0), (1.0, 2.0), (2.0, 1.0), (2.0, 2.0)])
 def system_kwargs(request: pytest.FixtureRequest) -> dict[str, float]:
+    """Keyword arguments initialising a KeplerSystem."""
     return dict(zip(("alpha", "mass"), request.param))
 
 
 @pytest.fixture(params=[1.0, 2.0])
 def parameter(request: pytest.FixtureRequest) -> float:
+    """Parameter for a conic section in polar coordinates."""
     return request.param
 
 
 @pytest.fixture(params=[False, True])
 def positive_angular_mom(request: pytest.FixtureRequest) -> bool:
+    """Whether the angular momentum is positive."""
     return request.param
 
 
 @pytest.fixture()
 def kepler_system(system_kwargs: "Mapping[str, float]") -> Kepler2DSystem:
+    """Give a Kepler system from the keyword argument."""
     return Kepler2DSystem(**system_kwargs)
 
 
@@ -44,32 +50,43 @@ def geometries(
     ecc: float,
     parameter: float,
 ) -> dict[str, bool | float]:
+    """Give geometric characteristics for initialising a KeplerSystem."""
     return dict(positive_angular_mom=positive_angular_mom, ecc=ecc, parameter=parameter)
 
 
 class Test2DSystem:
+    """Tests for the class Kepler2DSystem."""
+
     def test_init(self, system_kwargs: "Mapping[str, float]") -> None:
+        """Test initialising a KeplerSystem."""
         Kepler2DSystem(**system_kwargs)
 
     @pytest.mark.parametrize(("alpha", "mass"), [(-1, 1), (1, -1), (0, 1), (1, 0)])
     def test_raise(self, alpha: int, mass: int) -> None:
+        """Test raises upon invalid alpha or mass."""
         with pytest.raises(pydantic.ValidationError):
             Kepler2DSystem(alpha=alpha, mass=mass)
 
 
 class Test2DIoM:
+    """Tests for the class Kepler2DIoM."""
+
     def test_raise(self) -> None:
+        """Test raises upon an invalid angular momentum."""
         match = "Only non-zero angular momenta are supported"
         with pytest.raises(NotImplementedError, match=match):
             Kepler2DIoM(ene=1, angular_mom=0)
 
 
 class TestKepler2D:
+    """Tests for the class Kepler2D."""
+
     def test_from_geometry(
         self,
         system_kwargs: "Mapping[str, float]",
         geometries: "Mapping[str, bool | float]",
     ) -> None:
+        """Test initialising from geometric specifications."""
         kep = Kepler2D.from_geometry(system_kwargs, geometries)
         assert_almost_equal(kep.ecc, geometries["ecc"])
         assert_almost_equal(kep.parameter, geometries["parameter"])
@@ -80,6 +97,7 @@ def test_minimal_ene(
         system_kwargs: "Mapping[str, float]",
         geometries: "Mapping[str, bool | float]",
     ) -> None:
+        """Test the value of minimal energy, and that decresing this energy raises an exception."""
         kep = Kepler2D.from_geometry(system_kwargs, geometries)
         assert_equal(kep.ene, Kepler2D.minimal_ene(kep.angular_mom, system_kwargs))
         with pytest.raises(ValueError):
@@ -91,6 +109,7 @@ def test_period_from_u(
         system_kwargs: "Mapping[str, float]",
         geometries: "Mapping[str, bool | float]",
     ) -> None:
+        """Test the periods of t and tau are consistent."""
         kep = Kepler2D.from_geometry(system_kwargs, geometries)
         if ecc >= 0 and ecc < 1:
             assert_almost_equal(kep.period_in_tau, kep.period * kep.t_to_tau_factor)
@@ -110,6 +129,7 @@ def test_phi_of_u_tau_naive(
         system_kwargs: "Mapping[str, float]",
         geometries: "Mapping[str, bool | float]",
     ) -> None:
+        """Test special values of phi from u and tau."""
         kep = Kepler2D.from_geometry(system_kwargs, geometries)
         ecc = kep.ecc  # numeric instability
         if ecc >= 0 and ecc < 1:
@@ -124,6 +144,7 @@ def test_phi_of_u_tau_naive(
 
     @pytest.fixture(params=[_some_numbers[0], _some_numbers])
     def t(self, request: pytest.FixtureRequest) -> float | list[float]:
+        """Give a scalar t and a list of t's."""
         return request.param
 
     def test_r_and_phi(
@@ -132,6 +153,7 @@ def test_r_and_phi(
         system_kwargs: "Mapping[str, float]",
         geometries: "Mapping[str, bool | float]",
     ) -> None:
+        """Test r and phi satisfies the well-known solution of r of phi."""
         kep = Kepler2D.from_geometry(system_kwargs, geometries)
         tau = kep.tau(t)
         u = kep.u_of_tau(tau)
diff --git a/tests/models/kepler_problem/test_numerics.py b/tests/models/kepler_problem/test_numerics.py
index db8db27..89a8738 100644
--- a/tests/models/kepler_problem/test_numerics.py
+++ b/tests/models/kepler_problem/test_numerics.py
@@ -1,3 +1,5 @@
+"""Tests for the Kepler numerics module."""
+
 from typing import TYPE_CHECKING, Literal
 
 import numpy as np
@@ -13,6 +15,8 @@
 
 
 class TestUOfTau:
+    """Tests for calculating u from tau."""
+
     @pytest.mark.parametrize(
         "method",
         [
@@ -32,6 +36,7 @@ def test_u_of_tau(
         u_s: "npt.ArrayLike",
         tau_of_u: "Callable[[float, npt.ArrayLike], npt.ArrayLike]",
     ) -> None:
+        """Test numeric evaluation of u from tau."""
         u_s = np.array(u_s, copy=False)
         tau = tau_of_u(ecc, u_s)
         actual = u_of_tau(ecc, tau, method)
diff --git a/tests/models/test_brownian_motion.py b/tests/models/test_brownian_motion.py
index 846468e..c89911a 100644
--- a/tests/models/test_brownian_motion.py
+++ b/tests/models/test_brownian_motion.py
@@ -1,3 +1,5 @@
+"""Tests for the Brownian motion main module."""
+
 import numpy as np
 import pandas as pd
 import pytest
@@ -14,6 +16,7 @@
     [(0.1, np.array(0.1)), (1, np.array(1.0)), ([1, 2], np.array([1, 2]))],
 )
 def test_brownian_motion_ic(x0, expected):
+    """Test BrownianMotionIC."""
     brownian_motion_ic = BrownianMotionIC(x0=x0)
 
     np.testing.assert_equal(brownian_motion_ic.x0, expected)
@@ -24,6 +27,7 @@ def test_brownian_motion_ic(x0, expected):
     [(1, 1, 1), (1, 2, 2), (2, 1, 4)],
 )
 def test_brownian_motion_system(sigma, delta_t, gaussian_scale):
+    """Test BrownianMotionSystem."""
     bms = BrownianMotionSystem(sigma=sigma, delta_t=delta_t)
 
     assert bms.gaussian_scale == gaussian_scale
@@ -37,6 +41,7 @@ def test_brownian_motion_system(sigma, delta_t, gaussian_scale):
     ],
 )
 def test_brownian_motion_system_failed_spec(sigma, delta_t):
+    """Test raises upon illegal parameters for initialising a BrownianMotionSystem."""
     with pytest.raises(ValueError):
         BrownianMotionSystem(sigma=sigma, delta_t=delta_t)
 
@@ -64,6 +69,7 @@ def test_brownian_motion_system_failed_spec(sigma, delta_t):
     ],
 )
 def test_brownian_motion_generate_from(sigma, x0, expected):
+    """Test BrownianMotion values from generate_from, comparing with calculation results of the author."""
     system = {
         "sigma": sigma,
         "delta_t": 1,
@@ -118,6 +124,7 @@ def test_brownian_motion_generate_from(sigma, x0, expected):
     ],
 )
 def test_brownian_motion(sigma, x0, t, expected):
+    """Test BrownianMotion values from __call__, comparing with calculation results of the author."""
     system = {
         "sigma": sigma,
         "delta_t": 1,
diff --git a/tests/models/test_free_particle.py b/tests/models/test_free_particle.py
index aa64afd..c9ef078 100644
--- a/tests/models/test_free_particle.py
+++ b/tests/models/test_free_particle.py
@@ -1,3 +1,5 @@
+"""Tests for the free particle main module."""
+
 from typing import Mapping, Sequence
 
 import pandas as pd
@@ -8,12 +10,15 @@
 
 
 class TestFreeParticleIC:
+    """Tests for the class FreeParticleIC."""
+
     @pytest.mark.parametrize(("x0", "v0"), [(1, 2), ((1,), (2,)), ((1, 2), (2, 3))])
     def test_constructor(
         self,
         x0: int | Sequence[int],
         v0: int | Sequence[int],
     ) -> None:
+        """Test initialising a FreeParticleIC."""
         assert FreeParticleIC(x0=x0, v0=v0)
 
     @pytest.mark.parametrize(
@@ -26,11 +31,14 @@ def test_raise(
         v0: Sequence[int],
         expected: type[Exception],
     ) -> None:
+        """Test raise upon inconsistent initial conditions."""
         with pytest.raises(expected):
             FreeParticleIC(x0=x0, v0=v0)
 
 
 class TestFreeParticle:
+    """Tests for the class FreeParticle."""
+
     @pytest.mark.parametrize(
         ("x0", "v0", "expected"),
         [
@@ -44,6 +52,7 @@ def test_definition(
         v0: int | Sequence[int],
         expected: Mapping[str, Mapping[str, int | Sequence[int]]],
     ) -> None:
+        """Test the definition property."""
         assert FreeParticle(initial_condition=dict(x0=x0, v0=v0)).definition == expected
 
     @pytest.mark.parametrize(
@@ -65,6 +74,7 @@ def test_call(
         t: int | Sequence[int],
         expected: pd.DataFrame,
     ) -> None:
+        """Test the __call__ interface."""
         assert_frame_equal(
             FreeParticle(initial_condition=dict(x0=x0, v0=v0))(t).astype(float),
             expected.astype(float),
diff --git a/tests/models/test_harmonic_oscillator.py b/tests/models/test_harmonic_oscillator.py
index 9ba262c..4c53762 100644
--- a/tests/models/test_harmonic_oscillator.py
+++ b/tests/models/test_harmonic_oscillator.py
@@ -1,3 +1,5 @@
+"""Tests for the harmonic oscillator main module."""
+
 from typing import Mapping, Sequence
 
 import numpy as np
@@ -14,19 +16,24 @@
 )
 
 
-@pytest.mark.parametrize("zeta", [(-0.5), (-2.0)])
-def test_harmonic_oscillator_system_damping_zeta(zeta):
-    with pytest.raises(ValueError):
-        HarmonicOscillatorSystem(omega=1, zeta=zeta)
+@pytest.fixture
+def omega() -> int:
+    """Give omega."""
+    return 1
+
 
+@pytest.mark.parametrize("zeta", [-0.5, -2.0])
+def test_system_damping_zeta(omega: int, zeta: float) -> None:
+    """Test raises from HarmonicOscillatorSystem upon illegal zeta."""
     with pytest.raises(ValueError):
-        SimpleHarmonicOscillator(system={"omega": 1, "zeta": zeta})
+        HarmonicOscillatorSystem(omega=omega, zeta=zeta)
 
 
-@pytest.mark.parametrize("zeta", [(0.5), (1.0)])
-def test_simple_harmonic_oscillator_instantiation(zeta):
+@pytest.mark.parametrize("zeta", [-2.0, -0.5, 0.5, 1.0])
+def test_oscillator_damping_zeta(omega: int, zeta: float) -> None:
+    """Test raises from SimpleHarmonicOscillator upon illegal zeta."""
     with pytest.raises(ValueError):
-        SimpleHarmonicOscillator(system={"omega": 1, "zeta": zeta})
+        SimpleHarmonicOscillator(system={"omega": omega, "zeta": zeta})
 
 
 @pytest.mark.parametrize(
@@ -50,6 +57,7 @@ def test_simple_harmonic_oscillator_instantiation(zeta):
     ],
 )
 def test_simple_harmonic_oscillator(omega, expected):
+    """Test SimpleHarmonicOscillator from periods and the number of samples, comparing with caculations of the author."""
     ho = SimpleHarmonicOscillator(system={"omega": omega})
 
     df = ho(n_periods=1, n_samples_per_period=10)
@@ -79,6 +87,7 @@ def test_simple_harmonic_oscillator(omega, expected):
     ],
 )
 def test_underdamped_harmonic_oscillator(omega, zeta, expected):
+    """Test under DampedHarmonicOscillator from periods and the number of samples, comparing with caculations of the author."""
     ho = DampedHarmonicOscillator(system={"omega": omega, "zeta": zeta})
 
     df = ho(n_periods=1, n_samples_per_period=10)
@@ -108,6 +117,7 @@ def test_underdamped_harmonic_oscillator(omega, zeta, expected):
     ],
 )
 def test_overdamped_harmonic_oscillator(omega, zeta, expected):
+    """Test over DampedHarmonicOscillator from periods and the number of samples, comparing with caculations of the author."""
     ho = DampedHarmonicOscillator(system={"omega": omega, "zeta": zeta})
 
     df = ho(n_periods=1, n_samples_per_period=10)
@@ -137,6 +147,7 @@ def test_overdamped_harmonic_oscillator(omega, zeta, expected):
     ],
 )
 def test_criticaldamped_harmonic_oscillator(omega, zeta, expected):
+    """Test critical DampedHarmonicOscillator from periods and the number of samples, comparing with caculations of the author."""
     ho = DampedHarmonicOscillator(system={"omega": omega, "zeta": zeta})
 
     df = ho(n_periods=1, n_samples_per_period=10)
@@ -145,13 +156,19 @@ def test_criticaldamped_harmonic_oscillator(omega, zeta, expected):
 
 
 class TestComplexHarmonicOscillatorIC:
+    """Tests for the class ComplexHarmonicOscillatorIC."""
+
     @pytest.mark.parametrize("kwargs", [dict(x0=(1, 2), phi=(2, 3)), dict(x0=(1, 2))])
     def test_ic(self, kwargs: Mapping[str, tuple[int, int]]) -> None:
+        """Test initialising ComplexSimpleHarmonicOscillatorIC."""
         assert ComplexSimpleHarmonicOscillatorIC(**kwargs)
 
 
 class TestComplexHarmonicOscillator:
+    """Tests for the class ComplexHarmonicOscillator."""
+
     def test_complex(self) -> None:
+        """Test initialising ComplexSimpleHarmonicOscillator."""
         assert ComplexSimpleHarmonicOscillator(
             dict(omega=3),
             dict(x0=(1, 2), phi=(2, 3)),
@@ -159,6 +176,7 @@ def test_complex(self) -> None:
 
     @pytest.mark.parametrize("zeta", [0.5, 1.0, 1.5])
     def test_raise(self, zeta: float) -> None:
+        """Test raises from ComplexSimpleHarmonicOscillator upon illegal zeta."""
         with pytest.raises(ValueError):
             ComplexSimpleHarmonicOscillator(
                 dict(omega=3, zeta=zeta),
@@ -167,6 +185,7 @@ def test_raise(self, zeta: float) -> None:
 
     @pytest.fixture(params=(1, (1,), [1, 2], np.array([2, 3, 5, 7, 11])))
     def times(self, request: pytest.FixtureRequest) -> int | Sequence[int]:
+        """Give scalar time, Sequences of and numpy array of times."""
         return request.param
 
     @pytest.mark.parametrize("omega", [3, 5])
@@ -179,6 +198,7 @@ def test_degenerate_real(
         phi: int,
         times: int | Sequence[int],
     ) -> None:
+        """Test the degenerate case where ComplexSimpleHarmonicOscillator reduces to SimpleHarmonicOscillator."""
         csho = ComplexSimpleHarmonicOscillator(
             dict(omega=omega),
             dict(x0=(x0, x0), phi=(phi, phi)),
diff --git a/tests/models/test_harmonic_oscillator_chain.py b/tests/models/test_harmonic_oscillator_chain.py
index 2a6b7f6..fe12fda 100644
--- a/tests/models/test_harmonic_oscillator_chain.py
+++ b/tests/models/test_harmonic_oscillator_chain.py
@@ -1,3 +1,5 @@
+"""Tests for the harmonic oscillator chain main module."""
+
 from itertools import chain, product
 from typing import Iterable, Mapping, Sequence
 
@@ -17,12 +19,16 @@
 
 
 class TestHarmonicOscillatorChain:
+    """Tests for the class harmonic oscillator chain."""
+
     @pytest.fixture(params=(1, 2))
     def omega(self, request: pytest.FixtureRequest) -> int:
+        """Give omega for the system."""
         return request.param
 
     @pytest.fixture(params=((0, 0), (0, 1), (1, 0), (1, 1)))
     def free_mode(self, request: pytest.FixtureRequest) -> dict[str, int]:
+        """Give initial conditions of the free mode of the system."""
         return dict(zip(("x0", "v0"), request.param))
 
     @pytest.fixture(
@@ -34,10 +40,15 @@ def wave_modes(
         self,
         request: pytest.FixtureRequest,
     ) -> list[dict[str, tuple[int, int]]]:
+        """Give initial conditions of the wave mode of the system.
+
+        There can be 0, 1, 2 or 3 independent wave modes.
+        """
         return request.param
 
     @pytest.fixture(params=(False, True))
     def odd_dof(self, request: pytest.FixtureRequest) -> bool:
+        """Give whether the system has an odd number of DoF."""
         return request.param
 
     @pytest.fixture()
@@ -46,10 +57,12 @@ def legal_wave_modes_and_odd_def(
         wave_modes: Iterable[Mapping[str, tuple[int, int]]],
         odd_dof: bool,
     ) -> tuple[Iterable[Mapping[str, tuple[int, int]]], bool]:
+        """Give legal wave-modes initial conditions, taking the odd_dof parameter into account."""
         return wave_modes if odd_dof else chain(wave_modes, [dict(amp=(1, 1))]), odd_dof
 
     @pytest.fixture(params=(0, 1, (0, 1)))
     def times(self, request: pytest.FixtureRequest) -> int | tuple[int]:
+        """Give a scalar time or a tuple of times."""
         return request.param
 
     def test_init(
@@ -61,6 +74,7 @@ def test_init(
             bool,
         ],
     ) -> None:
+        """Test initialising a HarmonicOscillatorsChain."""
         wave_modes, odd_dof = legal_wave_modes_and_odd_def
         assert HarmonicOscillatorsChain(omega, [free_mode, *wave_modes], odd_dof)
 
@@ -75,6 +89,7 @@ def hoc_and_zs(
         ],
         times: int | Sequence[int],
     ) -> tuple[HarmonicOscillatorsChain, np.ndarray, np.ndarray]:
+        """Give initialised harmonic oscillators chain and the solutions from times."""
         wave_modes, odd_dof = legal_wave_modes_and_odd_def
         hoc = HarmonicOscillatorsChain(omega, [free_mode, *wave_modes], odd_dof)
         return (hoc, *hoc._z(times))
@@ -83,6 +98,7 @@ def test_real(
         self,
         hoc_and_zs: tuple[HarmonicOscillatorsChain, np.ndarray, np.ndarray],
     ) -> None:
+        """Test that the solution are real."""
         _, original_zs, _ = hoc_and_zs
         assert np.all(original_zs.imag == 0.0)
 
@@ -90,6 +106,7 @@ def test_dof(
         self,
         hoc_and_zs: tuple[HarmonicOscillatorsChain, np.ndarray, np.ndarray],
     ) -> None:
+        """Test the number of degrees of freedom is consistent."""
         hoc, original_zs, _ = hoc_and_zs
         assert original_zs.shape[0] == hoc.n_dof
 
@@ -100,6 +117,7 @@ def test_raise(
         free_mode: Mapping[str, int],
         wave_mode: Mapping[str, int],
     ) -> None:
+        """Test raises when odd_dof is false but the standing-wave mode is not real."""
         ics = [free_mode, *([wave_mode] if wave_mode else [])]
         with pytest.raises(ValueError):
             HarmonicOscillatorsChain(omega, ics, False)
diff --git a/tests/models/test_pendulum.py b/tests/models/test_pendulum.py
index 6fb643c..fbaa999 100644
--- a/tests/models/test_pendulum.py
+++ b/tests/models/test_pendulum.py
@@ -1,3 +1,5 @@
+"""Tests for the pendulum main module."""
+
 import math
 from typing import Sequence
 
@@ -11,43 +13,56 @@
 
 @pytest.fixture(params=[0.3, 0.6, 1.5])
 def omega0(request: pytest.FixtureRequest) -> float:
+    """Give a few omega0's."""
     return request.param
 
 
 @pytest.fixture(params=[-1.5, -0.3, 0.2, 0.4, +1.4])
 def theta0(request: pytest.FixtureRequest) -> float:
+    """Give a few theta0's."""
     return request.param
 
 
 @pytest.fixture(params=(-1.2, [-5.5, 0.0, 0.5, 1.0], np.linspace(0, 1, 7)))
 def times(request: pytest.FixtureRequest) -> float:
+    """Give a scalar time, a list of and a numpy array of times."""
     return request.param
 
 
 class TestPendulumSystem:
+    """Tests for the class PendulumSystem."""
+
     @pytest.mark.parametrize("omega0", [-1.0, 0])
     def test_omega0_range(self, omega0: float) -> None:
+        """Test raises upon illegal omega0's."""
         with pytest.raises(ValueError):
             _ = PendulumSystem(omega0=omega0)
 
 
 class TestPendulumIC:
+    """Tests for the class PendulumIC."""
+
     @pytest.mark.parametrize("theta0", [-2.0, 2.0])
     def test_theta0_range(self, theta0: float) -> None:
+        """Test raises upon illegal theta0's."""
         with pytest.raises(ValueError):
             _ = PendulumIC(theta0=theta0)
 
     def test_k(self, theta0: float) -> None:
+        """Test the kinematic transformation from theta0 to k."""
         ic = PendulumIC(theta0=theta0)
         assert ic.k == np.sin(ic.theta0 / 2)
 
 
 class TestPendulum:
+    """Tests for the class Pendulum."""
+
     @pytest.mark.parametrize(
         ("omega0", "theta0", "period"),
         [(2 * math.pi, math.pi / 3, 1.0731820071493643751)],
     )
     def test_period_static(self, omega0: float, theta0: float, period: float) -> None:
+        """Test the calculated period from python is the same as the value from Mathematica."""
         p = Pendulum(omega0, theta0)
         assert pytest.approx(p.period) == period
 
@@ -57,6 +72,7 @@ def test_transf(
         theta0: float,
         times: "Sequence[float] | npt.ArrayLike",
     ) -> None:
+        """Test the dynamic transformation from u to theta."""
         p = Pendulum(omega0, theta0)
         arr_times = np.asarray(times)
 
@@ -70,6 +86,7 @@ def test_period_dynamic_theta(
         theta0: float,
         times: "Sequence[float] | npt.ArrayLike",
     ) -> None:
+        """Test theta of t has the theory-predicted period."""
         p = Pendulum(omega0, theta0)
         arr_times_1 = np.array(times, copy=False) + p.period
 
@@ -78,6 +95,7 @@ def test_period_dynamic_theta(
         assert_array_almost_equal(theta, theta_1)
 
     def test_generate_from(self, omega0: float, theta0: float) -> None:
+        """Test alternative UI with periods and sample numbers."""
         p = Pendulum(omega0, theta0)
 
         df = p.generate_from(n_periods=1, n_samples_per_period=10)