Repository: jorgepiloto/lamberthub (language: python)

Issue:
Implement gooding1990 solver
💻 **Solver request: gooding1990**

This algorithm was originally devised by Gooding in 1990. Many of the work is built on top of Lancaster's one, but an accurate guess is made so fast convergence to the solution is achieved. In addition, the root finder used is Halley's method, so this solver is expected to be a high robust one. Gooding's routine was also extended to the multi-revolution case, so a full algorithm is available.

We should implement this under the name of `gooding1990`.  Because this algorithm is the first of the library, let us create a sub-package called `universal` in which all the universal-formulae solvers are included in the sense of Python modules.

📝 **References**

The algorithm appeared for the first time in 1988 Gooding's publication, see the BibTeX citation:

```
@techreport{gooding1988solution,
  title={On the solution of Lambert's orbital boundary-value problem},
  author={Gooding, RH},
  year={1988},
  institution={ROYAL AEROSPACE ESTABLISHMENT FARNBOROUGH (UNITED KINGDOM)}
}
```
However, the implementation we must follow is the one comming from the revised revision, which was published two years later in 1990:

```
@article{gooding1990,
  title={A procedure for the solution of Lambert's orbital boundary-value problem},
  author={Gooding, RH},
  journal={Celestial Mechanics and Dynamical Astronomy},
  volume={48},
  number={2},
  pages={145--165},
  year={1990},
  publisher={Springer}
}
```

Source code under Fortran77 is attached at the end of this article, so we just need to update it to modern Python 🚀 



Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index a7f6598..d9c0661 100644
--- a/README.md
+++ b/README.md
@@ -33,7 +33,7 @@ print([solver.__name__ for solver in ALL_SOLVERS])
 At this moment, the following algorithms are available:
 
 ```bash
->>> []
+>>> ['gooding1990']
 ```
 
 ## How can I use a solver?
diff --git a/docs/source/how_to_guides/how_to_use.md b/docs/source/how_to_guides/how_to_use.md
index 21ca96a..058dd19 100644
--- a/docs/source/how_to_guides/how_to_use.md
+++ b/docs/source/how_to_guides/how_to_use.md
@@ -24,17 +24,11 @@ that is: `authorYYYY`.
 To answer this question, simply run the following code snippet or refer to the
 official package API reference:
 
-```python
+```{code-cell}
 from lamberthub import ALL_SOLVERS
 print([solver.__name__ for solver in ALL_SOLVERS])
 ```
 
-and these are the available solvers for current version:
-
-```bash
->>> []
-```
-
 In addition, `lamberthub` provides other lists holding algorithms which
 present particular features such us multi-revolutions or high-robustness. These
 macros are listed down:
@@ -124,6 +118,38 @@ The orbit is seen to be prograde (inclination is less than $90^{\circ}$) and
 direct $M=0$. Remember that when $M=0$, there is only one possible solution, so
 the `low_path` flag will not play any role in this problem.
 
+To solve for the problem, we first import a solver. Let us use Gooding's one:
+
+```{code-cell}
+from lamberthub import gooding1990
+```
+
+Next, we specify the initial conditions of the problem:
+
+```{code-cell}
+# Import NumPy for declaring position vectors
+import numpy as np
+
+# Initial conditions for the problem
+mu_sun = 39.47692641  # [AU ** 3 / year ** 2]
+r1 = np.array([0.159321004, 0.579266185, 0.052359607])  # [AU]
+r2 = np.array([0.057594337, 0.605750797, 0.068345246])  # [AU]
+tof = 0.010794065  # [year]
+```
+
+Finally, we just need to solve for it. Notice that, as explained before, the
+default value for the `prograde` flag is `True`, which matches the one from
+problem's statement.
+
+```{code-cell}
+# Solving the problem
+v1, v2 = gooding1990(mu_sun, r1, r2, tof)
+
+# Let us print the results
+print(f"Initial velocity: {v1} [AU / years]\nFinal velocity: {v2} [AU / years]")
+```
+
+previous values are the same ones comming from the original example.
 
 [^1]: Directly taken from *An Introduction to the Mathematics and Methods of
   Astrodynamics, revised edition*, by R.H. Battin, problem 7-12.
diff --git a/src/lamberthub/__init__.py b/src/lamberthub/__init__.py
index cf719d7..ecb8772 100644
--- a/src/lamberthub/__init__.py
+++ b/src/lamberthub/__init__.py
@@ -1,15 +1,17 @@
 """ A collection of Lambert's problem solvers """
 
+from lamberthub.universal_solvers.gooding import gooding1990
+
 __version__ = "0.1.dev0"
 
-ALL_SOLVERS = []
+ALL_SOLVERS = [gooding1990]
 """ A list holding all lamberthub available solvers """
 
-ZERO_REV_SOLVERS = []
+ZERO_REV_SOLVERS = [gooding1990]
 """ A list holding all direct-revolution lamberthub solvers """
 
-MULTI_REV_SOLVERS = []
+MULTI_REV_SOLVERS = [gooding1990]
 """ A list holding all multi-revolution lamberthub solvers """
 
-ROBUST_SOLVERS = []
+ROBUST_SOLVERS = [gooding1990]
 """ A list holding all robust lamberthub solvers """
diff --git a/src/lamberthub/universal_solvers/__init__.py b/src/lamberthub/universal_solvers/__init__.py
new file mode 100644
index 0000000..9274111
--- /dev/null
+++ b/src/lamberthub/universal_solvers/__init__.py
@@ -0,0 +1,1 @@
+""" A sub-package holding all universal formulae based solvers """
diff --git a/src/lamberthub/universal_solvers/gooding.py b/src/lamberthub/universal_solvers/gooding.py
new file mode 100644
index 0000000..7479d88
--- /dev/null
+++ b/src/lamberthub/universal_solvers/gooding.py
@@ -0,0 +1,658 @@
+""" A module hosting all algorithms devised by Gooding  """
+
+import numpy as np
+from numpy.linalg import norm
+
+from lamberthub.utils.angles import get_orbit_normal_vector, get_transfer_angle
+from lamberthub.utils.assertions import (
+    assert_parameters_are_valid,
+    assert_transfer_angle_not_zero,
+)
+
+
+def gooding1990(
+    mu,
+    r1,
+    r2,
+    tof,
+    M=0,
+    prograde=True,
+    low_path=True,
+    maxiter=35,
+    atol=1e-5,
+    rtol=1e-7,
+    full_output=False,
+):
+    r"""
+    Lambert's problem solver using the method proposed by R. H. Gooding in 1990.
+
+    Parameters
+    ----------
+    mu: float
+        Gravitational parameter, equivalent to :math:`GM` of attractor body.
+    r1: numpy.array
+        Initial position vector.
+    r2: numpy.array
+        Final position vector.
+    M: int
+        Number of revolutions. Must be equal or greater than 0 value.
+    prograde: bool
+        If `True`, specifies prograde motion. Otherwise, retrograde motion is imposed.
+    low_path: bool
+        If two solutions are available, it selects between high or low path.
+    maxiter: int
+        Maximum number of iterations.
+    atol: float
+        Absolute tolerance.
+    rtol: float
+        Relative tolerance.
+    full_output: bool
+        If True, the number of iterations is also returned.
+
+    Returns
+    -------
+    v1: numpy.array
+        Initial velocity vector.
+    v2: numpy.array
+        Final velocity vector.
+    numiter: list
+        Number of iterations.
+
+    Notes
+    -----
+    This module holds the Lambert's problem solver devised by R. H. Gooding in
+    his technical report [1]_ originally published in 1988. However, the
+    implementation corresponds to the one proposed by the author a couple of
+    years later in his article [2]_ from 1990. Some improvements to the
+    originally algorithm were also made by Klumpp in his performance comparison
+    [3]_ between Lamberts problem solvers.  Those have been added to this code
+    to prevent failures for particular inputs.  The result is a fully working
+    algorithm for both single and multi-revolution orbits.
+
+    The code has been kept as close as possible to the original FORTRAN-77  one.
+    However, some statements (like "goto line" ones) have been deprecated as
+    they introduce "spaghetti code". Since the original implementation imposed a
+    relative tolerance together with the number of iterations, these parameters
+    have been modified so the user can freely choose their values.
+
+    References
+    ----------
+    .. [1] Gooding, R. H. (1988). On the solution of Lambert's orbital
+           boundary-value problem. ROYAL AEROSPACE ESTABLISHMENT FARNBOROUGH (UNITED
+           KINGDOM).
+
+    .. [2] Gooding, R. H. (1990). A procedure for the solution of Lambert's orbital
+           boundary-value problem. Celestial Mechanics and Dynamical Astronomy, 48(2),
+           145-165.
+
+    .. [3] Klumpp,  A. (1999). Performance  Comparison  of  Lambert  and  Kepler
+           Algorithms,  Interoffice  Memorandum, JPL.
+
+    """
+
+    # Check that input parameters are safe
+    assert_parameters_are_valid(mu, r1, r2, tof, M)
+
+    # Norm of the initial and final position vectors.
+    r1_norm, r2_norm = [norm(r) for r in [r1, r2]]
+    i_r1, i_r2 = [r / r_norm for r, r_norm in zip([r1, r2], [r1_norm, r2_norm])]
+
+    # Compute the real transfer angle according to sense of motion. Raise an
+    # exception its value is found to be the zero.
+    theta = get_transfer_angle(r1, r2, prograde)
+    assert_transfer_angle_not_zero(theta)
+
+    # Include additional revolutions if necessary.
+    dtheta = 2.0 * np.pi * M + theta
+
+    # Compute a vector normal to orbit plane, parallel to angular momentum. Need
+    # to specify orbit motion sense.
+    i_h = get_orbit_normal_vector(r1, r2, prograde)
+
+    # Compute the tangential unitary vectors at initial and final position vectors.
+    i_t1, i_t2 = [np.cross(i_h, i_r) for i_r in [i_r1, i_r2]]
+
+    # Compute amount of solutions and the velocity components.
+    vr1, vt1, vr2, vt2, numiter = vlamb(
+        mu,
+        r1_norm,
+        r2_norm,
+        dtheta,
+        tof,
+        low_path,
+        maxiter,
+        atol,
+        rtol,
+    )
+
+    # Final velocity vectors are the vector summation of radial and tangential.
+    v1 = vr1 * i_r1 + vt1 * i_t1
+    v2 = vr2 * i_r2 + vt2 * i_t2
+
+    return (v1, v2, numiter) if full_output is True else (v1, v2)
+
+
+def tlamb(m, q, qsqfm1, x, n):
+    """
+    Auxiliary routine for computing the non-dimensional time of flight as
+    function of the number of revolutions, the transfer parameter and the
+    independent variable.
+
+    Parameters
+    ----------
+    m: float
+        Number of revolutions.
+    q: float
+        The transfer angle parameter.
+    qsqfm1:
+        Equivalent to :math:`1-q^2`.
+    x: float
+        The independent variable.
+    n: float
+        Number of output parameters to be returned.
+
+    Returns
+    -------
+    t: float
+        Non-dimensional time evaluated at :math:`x`.
+    dt: float
+        First derivative of the non-dimensional time evaluated at :math:`x`.
+    d2t: float
+        Second derivative of the non-dimensional time evaluated at :math:`x`.
+    d3t: float
+        Third derivative of the non-dimensional time evaluated at :math:`x`.
+
+    """
+
+    # Define necessary parameters.
+    sw = 0.4
+    lm1 = n == -1
+    l1 = n >= 1
+    l2 = n >= 2
+    l3 = n == 3
+    qsq = q * q
+    xsq = x * x
+    u = (1.0 - x) * (1.0 + x)
+
+    if not lm1:
+        # Needed if series, and otherwise useful when z = 0.
+        dt, d2t, d3t = 0.0, 0.0, 0.0
+
+    if lm1 or m > 0 or x < 0.0 or np.abs(u) > sw:
+        # Direct computation, series is not required.
+        y = np.sqrt(np.abs(u))
+        z = np.sqrt(qsqfm1 + qsq * xsq)
+        qx = q * x
+
+        if qx <= 0.0:
+            a = z - qx
+            b = q * z - x
+
+        if qx <= 0.0 and lm1:
+            aa = qsqfm1 / a
+            bb = qsqfm1 * (qsq * u - xsq) / b
+
+        if (qx == 0.0 and lm1) or (qx > 0.0):
+            aa = z + qx
+            bb = q * z + x
+
+        if qx > 0.0:
+            a = qsqfm1 / aa
+            b = qsqfm1 * (qsq * u - xsq) / bb
+
+        if lm1:
+            t, dt, d2t, d3t = 0, b, bb, aa
+        else:
+
+            if qx * u >= 0.0:
+                g = x * z + q * u
+            else:
+                g = (xsq - qsq * u) / (x * z - q * u)
+
+            f = a * y
+
+            if x <= 1.0:
+                t = m * np.pi + np.arctan2(f, g)
+            else:
+                if f > sw:
+                    t = np.log(f + g)
+                else:
+                    fg1 = f / (g + 1.0)
+                    term = 2.0 * fg1
+                    fg1sq = fg1 * fg1
+                    t = term
+                    twoi1 = 1.0
+
+                    # Emulate FORTRAN-77 do-loop by carrying out first iteration
+                    # outside of it.
+                    twoi1 = twoi1 + 2.0
+                    term = term * fg1sq
+                    told = t
+                    t = t + term / twoi1
+
+                    while t != told:
+                        twoi1 = twoi1 + 2.0
+                        term = term * fg1sq
+                        told = t
+                        t = t + term / twoi1
+
+            t = 2.0 * (t / y + b) / u
+
+            if l1 and z != 0.0:
+                qz = q / z
+                qz2 = qz * qz
+                qz = qz * qz2
+                dt = (3.0 * x * t - 4.0 * (a + qx * qsqfm1) / z) / u
+                if l2:
+                    d2t = (3.0 * t + 5.0 * x * dt + 4.0 * qz * qsqfm1) / u
+                if l3:
+                    d3t = (8.0 * dt + 7.0 * x * d2t - 12.0 * qz * qz2 * x * qsqfm1) / u
+
+    else:
+        # Compute by series.
+        u0i = 1.0
+
+        if l1:
+            u1i = 1.0
+        if l2:
+            u2i = 1.0
+        if l3:
+            u3i = 1.0
+
+        term = 4.0
+        tq = q * qsqfm1
+        i = 0
+
+        if q < 0.5:
+            tqsum = 1.0 - q * qsq
+
+        if q >= 0.5:
+            tqsum = (1.0 / (1.0 + q) + q) * qsqfm1
+
+        ttmold = term / 3.0
+        t = ttmold * tqsum
+
+        # Again, emulate do-loop from FORTRAN-77 by carrying out the first
+        # iteration out of it.
+
+        i = i + 1
+        p = i
+        u0i = u0i * u
+
+        if l1 and i > 1:
+            u1i = u1i * u
+        if l2 and i > 2:
+            u2i = u2i * u
+        if l3 and i > 3:
+            u3i = u3i * u
+
+        term = term * (p - 0.5) / p
+        tq = tq * qsq
+        tqsum = tqsum + tq
+        told = t
+        tterm = term / (2.0 * p + 3.0)
+        tqterm = tterm * tqsum
+        t = t - u0i * ((1.5 * p + 0.25) * tqterm / (p * p - 0.25) - ttmold * tq)
+        ttmold = tterm
+        tqterm = tqterm * p
+
+        if l1:
+            dt = dt + tqterm * u1i
+        if l2:
+            d2t = d2t + tqterm * u2i * (p - 1.0)
+        if l3:
+            d3t = d3t + tqterm * u3i * (p - 1.0) * (p - 2.0)
+
+        # The loop starts below this line in its second iteration.
+        while i < n or t != told:
+
+            i = i + 1
+            p = i
+            u0i = u0i * u
+
+            if l1 and i > 1:
+                u1i = u1i * u
+            if l2 and i > 2:
+                u2i = u2i * u
+            if l3 and i > 3:
+                u3i = u3i * u
+
+            term = term * (p - 0.5) / p
+            tq = tq * qsq
+            tqsum = tqsum + tq
+            told = t
+            tterm = term / (2.0 * p + 3.0)
+            tqterm = tterm * tqsum
+            t = t - u0i * ((1.5 * p + 0.25) * tqterm / (p * p - 0.25) - ttmold * tq)
+            ttmold = tterm
+            tqterm = tqterm * p
+
+            if l1:
+                dt = dt + tqterm * u1i
+            if l2:
+                d2t = d2t + tqterm * u2i * (p - 1.0)
+            if l3:
+                d3t = d3t + tqterm * u3i * (p - 1.0) * (p - 2.0)
+
+        if l3:
+            d3t = 8.0 * x * (1.5 * d2t - xsq * d3t)
+        if l2:
+            d2t = 2.0 * (2.0 * xsq * d2t - dt)
+        if l1:
+            dt = -2.0 * x * dt
+
+        t = t / xsq
+
+    return t, dt, d2t, d3t
+
+
+def xlamb(m, q, qsqfm1, tin, maxiter, atol, rtol):
+    r"""
+    Auxiliary routine for finding the independent variable as function of the
+    number of revolutions, the transfer angle parameter and the non-dimensional
+    time of flight.
+
+    Parameters
+    ----------
+    m: float
+        Number of revolutions.
+    q: float
+        The transfer angle parameter.
+    qsqfm1: float
+        Equivalent to :math:`1-q^2`.
+    tin: float
+        The actual non-dimensional time of flight.
+    maxiter: int
+        Maximum number of iterations.
+    atol: float
+        Desired absolute tolerance.
+    rtol: float
+        Desired relative tolerance.
+
+    Returns
+    -------
+    n_sol: int
+        Number of solutions.
+    x: float
+        First solution.
+    xpl: float
+        Second solution, if available.
+    numiter: int
+        Number of iterations.
+
+    """
+
+    # Declare auxiliary parameters.
+    xpl = 0
+    c0, c1, c2, c3, c41, c42 = 1.7, 0.5, 0.03, 0.15, 1.0, 0.24
+    thr2 = np.arctan2(qsqfm1, 2.0 * q) / np.pi
+
+    # Boolean variables to emulate original code as max as possible.
+    goto3 = False
+
+    # Auxiliary function for 8th root.
+    def d8rt(x):
+        return np.sqrt(np.sqrt(np.sqrt(x)))
+
+    # Start computing the initial guess. The process is different depending
+    # on the number of revolutions.
+    if m == 0:
+
+        # Single-rev starter from T (at x = 0) and bilinear usually.
+        n = 1
+
+        # Call TLAMB routine.
+        t0, dt, d2t, d3t = tlamb(m, q, qsqfm1, 0.0, 0)
+        tdiff = tin - t0
+
+        if tdiff <= 0.0:
+            x = t0 * tdiff / (-4.0 * tin)
+            # -4.0 is the value of dt for x = 0.
+        else:
+            x = -tdiff / (tdiff + 4.0)
+            w = x + c0 * np.sqrt(2.0 * (1.0 - thr2))
+
+            if w < 0.0:
+                x = x - np.sqrt(d8rt(-w)) * (x + np.sqrt(tdiff / (tdiff + 1.5 * t0)))
+
+            w = 4.0 / (4.0 + tdiff)
+            x = x * (1.0 + x * (c1 * w - c2 * x * np.sqrt(w)))
+
+    else:
+
+        # With multi-revolutions first get T(min) as basis for starter.
+        xm = 1.0 / (1.5 * (m + 0.5) * np.pi)
+
+        if thr2 < 0.5:
+            xm = d8rt(2.0 * thr2) * xm
+        if thr2 > 0.5:
+            xm = (2.0 - d8rt(2.0 - 2.0 * thr2)) * xm
+
+        # For locating Tmin, an iterative process is required. Original
+        # implementation imposed 12 iterations but they were not considered to
+        # be part of numerical routine as they belong to the initial guess.
+        # Here, we impose it not to exceeded the number of iterations
+        for numiter in range(1, maxiter + 1):
+
+            # Call TLAMB routine
+            tmin, dt, d2t, d3t = tlamb(m, q, qsqfm1, xm, 3)
+
+            if d2t == 0.0:
+                break
+
+            xmold = xm
+            xm = xm - dt * d2t / (d2t * d2t - dt * d3t / 2.0)
+
+            # Compute the absolute and relative tolerances and check if within
+            # desired range
+            if np.abs(xmold - xm) < rtol * np.abs(xmold) + atol:
+                break
+
+        # Check last
+        if (numiter + 1) == maxiter:
+            raise ValueError("Exceeded max iterations!")
+
+        # Now proceed from T(min) to full starter
+        tdiffm = tin - tmin
+
+        if tdiffm < 0.0:
+            # Exit if no solution with this m
+            raise ValueError("No feasible solution, try lower M!")
+
+        elif tdiffm == 0.0:
+            # Found minimum time condition
+            x = xm
+            n = 1
+            return n, x, None, numiter
+
+        else:
+
+            n = 3
+
+            if d2t == 0.0:
+                d2t = 6.0 * m * np.pi
+
+            x = np.sqrt(tdiffm / (d2t / 2.0 + tdiffm / (1.0 - xm) ** 2))
+            w = xm + x
+            w = w * 4.0 / (4.0 + tdiffm) + (1.0 - w) ** 2
+            x = (
+                x
+                * (
+                    1.0
+                    - (1.0 + m + c41 * (thr2 - 0.5))
+                    / (1.0 + c3 * m)
+                    * x
+                    * (c1 * w + c2 * x * np.sqrt(w))
+                )
+                + xm
+            )
+            d2t2 = d2t / 2.0
+
+            if x >= 1.0:
+                n = 1
+                # GOTO 3
+                goto3 = True
+
+    # --- THE ITERATION METHOD STARTS NOW ---
+    # The iterative Halley method starts now
+    # TODO: THIS IS THE HALLEY METHOD
+
+    while True:
+
+        # 5: LINE OF STATEMENT
+        if goto3 is False:
+
+            # Enable desired number of iterations
+            for numiter in range(1, maxiter + 1):
+
+                t, dt, d2t, d3t = tlamb(m, q, qsqfm1, x, 2)
+                t = tin - t
+
+                if dt != 0.0:
+                    xold = x
+                    x = x + t * dt / (dt * dt + t * d2t / 2.0)
+
+                # This line was included so the code can work with tolerances
+                x_atol, x_rtol = np.abs(x - xold), np.abs(x / xold - 1)
+                if x_atol <= atol and x_rtol <= rtol:
+                    break
+
+            if n != 3:
+                # Exit if only one solution normally when m = 0
+                return n, x, xpl, numiter
+
+            n = 2
+            xpl = x
+        else:
+            # Update the goto condition
+            goto3 = False
+
+        # 3: LINE OF STATEMENT
+        # Second multi-rev starter
+        t0, dt, d2t, d3t = tlamb(m, q, qsqfm1, 0.0, 0)
+        tdiff0 = t0 - tmin
+        tdiff = tin - t0
+
+        if tdiff <= 0.0:
+            x = xm - np.sqrt(tdiffm / (d2t2 - tdiffm * (d2t2 / tdiff0 - 1.0 / xm ** 2)))
+        else:
+            x = -tdiff / (tdiff + 4.0)
+            w = x + c0 * np.sqrt(2.0 * (1.0 - thr2))
+
+            if w < 0.0:
+                x = x - np.sqrt(d8rt(-w)) * (x + np.sqrt(tdiff / (tdiff + 1.5 * t0)))
+
+            w = 4.0 / (4.0 + tdiff)
+            x = x * (
+                1.0
+                + (1.0 + m + c42 * (thr2 - 0.5))
+                / (1.0 + c3 * m)
+                * x
+                * (c1 * w - c2 * x * np.sqrt(w))
+            )
+
+            if x <= -1.0:
+                n = n - 1
+
+                # No finite solution with x < xm
+                if n == 1:
+                    x = xpl
+
+    return n, x, xpl
+
+
+def vlamb(mu, r1_norm, r2_norm, dtheta, tof, low_path, maxiter, atol, rtol):
+    r"""
+    Auxiliary routine for computing the velocity vector components, both
+    radian and tangential ones.
+
+    Parameters
+    ----------
+    mu: float
+        Gravitational parameter, equivalent to :math:`GM` of attractor body.
+    r1_norm: float
+        Norm of the initial position vector.
+    r2_norm: float
+        Norm of the final position vector.
+    dtheta: float
+        Transfer angle betwen initial and final vectors.
+    tof: float
+        Time of flight betwen initial and final position vectors.
+    low_path: bool
+        If two solutions are available, it selects between high or low path.
+    maxiter: int
+        Maximum number of iterations.
+    atol: float
+        Absolute tolerance :math:`abs(x_{i+1} - x_{i})`
+    rtol: float
+        Relative tolerance :math:`abs(\frac{x_{i+1}}{x_{i}} - 1)`
+
+    Returns
+    -------
+    n_sol: int
+        Number of solutions
+    vri: float
+        Radial velocity component at the initial position vector.
+    vti: float
+        Tangential velocity component at the initial position vector.
+    vrf: float
+        Radial velocity component at the final position vector.
+    vtf: float
+        Tangential velocity component at the final position vector.
+    numiter: int
+        Number of iterations required to compute solution.
+
+    """
+
+    # The following yields m = 0 when th = 2pi exactly
+    # Neither this nor the original code works for
+    # th < 0.0
+    thr2 = dtheta
+    m = 0
+
+    while thr2 > 2 * np.pi:
+        thr2 = thr2 - 2 * np.pi
+        m = m + 1
+    thr2 = thr2 / 2.0
+
+    # Compute auxliary parameters
+    dr = r1_norm - r2_norm
+    r1r2 = r1_norm * r2_norm
+    r1r2th = 4.0 * r1r2 * np.sin(thr2) ** 2
+    csq = dr * dr + r1r2th
+    c = np.sqrt(csq)
+    s = (r1_norm + r2_norm + c) / 2.0
+    mus = np.sqrt(mu * s / 2.0)
+    qsqfm1 = c / s
+    q = np.sqrt(r1r2) * np.cos(thr2) / s
+
+    if c != 0.0:
+        rho = dr / c
+        sig = r1r2th / csq
+    else:
+        rho = 0.0
+        sig = 1.0
+
+    t = 4.0 * mus * tof / s ** 2
+
+    # Compute the number of solutions and
+    n_sol, x1_sol, x2_sol, numiter = xlamb(m, q, qsqfm1, t, maxiter, atol, rtol)
+
+    # Filter the solution
+    if n_sol > 1:
+        if low_path is True:
+            x_sol = np.max([x1_sol, x2_sol])
+        else:
+            x_sol = np.min([x1_sol, x2_sol])
+    else:
+        x_sol = x1_sol
+
+    # Compute radial and tangential velocity components
+    _, qzminx, qzplx, zplqx = tlamb(m, q, qsqfm1, x_sol, -1)
+    vt2 = mus * zplqx * np.sqrt(sig)
+    vr1 = mus * (qzminx - qzplx * rho) / r1_norm
+    vt1 = vt2 / r1_norm
+    vr2 = -mus * (qzminx + qzplx * rho) / r2_norm
+    vt2 = vt2 / r2_norm
+
+    return vr1, vt1, vr2, vt2, numiter
diff --git a/src/lamberthub/utils/__init__.py b/src/lamberthub/utils/__init__.py
new file mode 100644
index 0000000..cb9d78e
--- /dev/null
+++ b/src/lamberthub/utils/__init__.py
@@ -0,0 +1,1 @@
+""" A sub-package holding different utilities and auxiliary functions """
diff --git a/src/lamberthub/utils/angles.py b/src/lamberthub/utils/angles.py
new file mode 100644
index 0000000..a35226e
--- /dev/null
+++ b/src/lamberthub/utils/angles.py
@@ -0,0 +1,88 @@
+""" Utilities related to angles computations """
+
+import numpy as np
+from numpy import cross, dot
+from numpy.linalg import norm
+
+
+def get_transfer_angle(r1, r2, prograde):
+    """
+    Solves for the transfer angle being known the sense of rotation.
+
+    Parameters
+    ----------
+    r1: np.array
+        Initial position vector.
+    r2: np.array
+        Final position vector.
+    prograde: bool
+        If True, it assumes prograde motion, otherwise assumes retrograde.
+
+    Returns
+    -------
+    dtheta: float
+        Transfer angle in radians.
+
+    """
+
+    # Check if both position vectors are collinear. If so, check if the transfer
+    # angle is 0 or pi.
+    if np.all(np.cross(r1, r2) == 0):
+        return 0 if np.all(np.sign(r1) == np.sign(r2)) else np.pi
+
+    # Solve for a unitary vector normal to the vector plane. Its direction and
+    # sense the one given by the cross product (right-hand) from r1 to r2.
+    h = cross(r1, r2) / norm(np.cross(r1, r2))
+
+    # Compute the projection of the normal vector onto the reference plane.
+    alpha = dot(np.array([0, 0, 1]), h)
+
+    # Get the minimum angle (0 <= dtheta <= pi) between r1 and r2.
+    r1_norm, r2_norm = [norm(vec) for vec in [r1, r2]]
+    theta0 = np.arccos(dot(r1, r2) / (r1_norm * r2_norm))
+
+    # Fix the value of theta if necessary
+    if prograde is True:
+        dtheta = theta0 if alpha > 0 else 2 * np.pi - theta0
+    else:
+        dtheta = theta0 if alpha < 0 else 2 * np.pi - theta0
+
+    return dtheta
+
+
+def get_orbit_normal_vector(r1, r2, prograde):
+    """
+    Computes a unitary normal vector aligned with the specific angular momentum
+    one of the orbit.
+
+    Parameters
+    ----------
+    r1: np.array
+        Initial position vector.
+    r2: np.array
+        Final position vector.
+    prograde: bool
+        If True, it assumes prograde motion, otherwise assumes retrograde.
+
+    Returns
+    -------
+    i_h: np.array
+        Unitary vector aligned with orbit specific angular momentum.
+
+    """
+
+    # Compute the normal vector and its projection onto the vertical axis
+    i_h = np.cross(r1, r2) / norm(np.cross(r1, r2))
+
+    # Solve the projection onto the positive vertical direction of the
+    # fundamental plane.
+    alpha = dot(np.array([0, 0, 1]), i_h)
+
+    # An prograde orbit always has a positive vertical component of its specific
+    # angular momentum. Therefore, we just need to check for this condition
+    if prograde is True:
+        i_h = i_h if alpha > 0 else -i_h
+    else:
+        i_h = i_h if alpha < 0 else -i_h
+
+    return i_h
diff --git a/src/lamberthub/utils/assertions.py b/src/lamberthub/utils/assertions.py
new file mode 100644
index 0000000..645aff8
--- /dev/null
+++ b/src/lamberthub/utils/assertions.py
@@ -0,0 +1,133 @@
+""" A collection of checkers for raising custom exceptions if required """
+
+import numpy as np
+
+
+def assert_parameters_are_valid(mu, r1, r2, tof, M):
+    """
+    Checks if solver input parameters are safe.
+
+    """
+
+    # Run all required checks
+    assert_gravitational_parameter_is_positive(mu)
+    assert_position_vectors_are_valid(r1, r2)
+    assert_time_of_flight_is_positive(tof)
+    assert_number_of_revolutions_not_negative(M)
+
+    return True
+
+
+def assert_gravitational_parameter_is_positive(mu):
+    """
+    Checks if the gravitational parameter is positive.
+
+    Parameters
+    ----------
+    mu: float
+        Gravitational parameter
+
+    """
+
+    # Check positive gravitational parameter
+    if mu <= 0:
+        raise ValueError("Gravitational parameter must be positive!")
+    else:
+        return True
+
+
+def assert_position_vector_is_valid(r):
+    """
+    Checks if position vector has proper dimensions and is not the null one.
+
+    Parameters
+    ----------
+    r: np.array
+        Initial position vector.
+
+    """
+
+    # Check that vector belongs to three-dimensional space
+    if r.shape != (3,):
+        raise ValueError("Vector must be three-dimensional!")
+
+    if np.all(r == 0):
+        raise ValueError("Position vector cannot be the null vector [0,0,0]!")
+
+    return True
+
+
+def assert_position_vectors_are_valid(r1, r2):
+    """
+    Checks if position vectors are safe in dimension and values.
+
+    Parameters
+    ----------
+    r1: np.array
+        Initial position vector.
+    r2: np.array
+        Final position vector.
+
+    """
+
+    # Check if position vectors have proper dimensions
+    for r in [r1, r2]:
+        assert_position_vector_is_valid(r1)
+
+    # Check both vectors are different
+    if np.all(np.equal(r1, r2)):
+        raise ValueError("Initial and final position vectors cannot be equal!")
+
+    return True
+
+
+def assert_time_of_flight_is_positive(tof):
+    """
+    Checks if time of flight is positive.
+
+    Parameters
+    ----------
+    tof: float
+        Time of flight.
+
+    """
+
+    if tof <= 0:
+        raise ValueError("Time of flight must be positive!")
+    else:
+        return True
+
+
+def assert_number_of_revolutions_not_negative(M):
+    """
+    Checks if the number of revolutions is zero or positive, that is, it does
+    not have a negative value.
+
+    Parameters
+    ----------
+    M: int
+        Number of revolutions
+
+    """
+
+    if M < 0:
+        raise ValueError("Number of revolutions must be equal or greater than zero!")
+    else:
+        return True
+
+
+def assert_transfer_angle_not_zero(dtheta):
+    """
+    Checks if the transfer angle is the null value, if so, raises an exception.
+
+    Parameters
+    ----------
+    dtheta: float
+        Transfer angle value.
+
+    """
+
+    if dtheta == 0:
+        raise ValueError("Transfer angle was found to be zero!")
+    else:
+        return True


Test patch:
diff --git a/tests/test_all_solvers.py b/tests/test_all_solvers.py
new file mode 100644
index 0000000..7c7e258
--- /dev/null
+++ b/tests/test_all_solvers.py
@@ -0,0 +1,177 @@
+""" Basic unitary tests for checking all available solvers """
+
+import numpy as np
+import pytest
+from numpy.testing import assert_allclose
+
+from lamberthub import ALL_SOLVERS
+
+
+@pytest.mark.parametrize("solver", ALL_SOLVERS)
+def test_transfer_angle_is_zero_raises_exception(solver):
+    with pytest.raises(ValueError) as excinfo:
+        r1, r2 = [i * np.ones(3) for i in range(1, 3)]
+        solver(1.00, r1, r2, 1.00)
+    assert "Transfer angle was found to be zero!" in excinfo.exconly()
+
+
+@pytest.mark.parametrize("solver", ALL_SOLVERS)
+def test_case_from_vallado_book(solver):
+    """
+    Example 5.7 from Fundamentals of Astrodynamics and Applications (4th
+    Edition), by David A. Vallado
+    """
+
+    # Initial conditions
+    mu_earth = 3.986004418e5  # [km ** 3 / s ** 2]
+    r1 = np.array([15945.34, 0.0, 0.0])  # [km]
+    r2 = np.array([12214.83899, 10249.46731, 0.0])  # [km]
+    tof = 76.0 * 60  # [s]
+
+    # Solving the problem
+    v1, v2 = solver(mu_earth, r1, r2, tof, prograde=True, low_path=True)
+
+    # Expected final results
+    expected_v1 = np.array([2.058913, 2.915965, 0.0])  # [km / s]
+    expected_v2 = np.array([-3.451565, 0.910315, 0.0])  # [km / s]
+
+    # Assert the results
+    assert_allclose(v1, expected_v1, rtol=1e-5)
+    assert_allclose(v2, expected_v2, rtol=1e-4)
+
+
+@pytest.mark.parametrize("solver", ALL_SOLVERS)
+def test_case_from_curtiss_book(solver):
+    """
+    Example 5.2 from Orbital Mechanics for Engineering Students (3rd
+    Edition), by Howard D. Curtiss
+    """
+
+    # Initial conditions
+    mu_earth = 3.986004418e5  # [km ** 3 / s ** 2]
+    r1 = np.array([5000.0, 10000.0, 2100.0])  # [km]
+    r2 = np.array([-14600.0, 2500.0, 7000.0])  # [km]
+    tof = 3600  # [s]
+
+    # Solve the problem
+    v1, v2 = solver(mu_earth, r1, r2, tof, prograde=True)
+
+    # Expected final results
+    expected_v1 = np.array([-5.9925, 1.9254, 3.2456])  # [km / s]
+    expected_v2 = np.array([-3.3125, -4.1966, -0.38529])  # [ km / s ]
+
+    # Assert the results
+    assert_allclose(v1, expected_v1, rtol=1e-4)
+    assert_allclose(v2, expected_v2, rtol=1e-4)
+
+
+@pytest.mark.parametrize("solver", ALL_SOLVERS)
+def test_case_from_battin_book(solver):
+    """
+    Example 7.12 from An Introduction to the Mathematics and Methods of
+    Astrodynamics (Revised Edition), by Richard H. Battin
+    """
+    # Initial conditions
+    mu_sun = 39.47692641  # [AU ** 3 / year ** 2]
+    r1 = np.array([0.159321004, 0.579266185, 0.052359607])  # [AU]
+    r2 = np.array([0.057594337, 0.605750797, 0.068345246])  # [AU]
+    tof = 0.010794065  # [year]
+
+    # Solving the problem
+    v1, v2 = solver(mu_sun, r1, r2, tof)
+
+    # Expected final results
+    expected_v1 = np.array([-9.303603251, 3.018641330, 1.536362143])  # [AU / year]
+
+    # Assert the results
+    assert_allclose(v1, expected_v1, rtol=1e-5)
+
+
+TABLE_OF_TRANSFERS_I = {
+    "M0_prograde_low": [
+        np.array([-2.09572809, 3.92602196, -4.94516810]),  # [km / s]
+        np.array([2.46309613, 0.84490197, 6.10890863]),  # [km / s]
+    ],
+    "M0_retrograde_high": [
+        np.array([1.94312182, -4.35300015, 4.54630439]),  # [km / s]
+        np.array([-2.38885563, -1.42519647, -5.95772225]),  # [km / s]
+    ],
+}
+
+
+@pytest.mark.parametrize("solver", ALL_SOLVERS)
+@pytest.mark.parametrize("case", TABLE_OF_TRANSFERS_I)
+def test_case_from_der_article_I(solver, case):
+    """
+    Example from Astrodynamics 102, by Gim J. Der. see:
+    http://derastrodynamics.com/docs/astrodynamics_102_v2.pdf
+    """
+
+    # Initial conditions
+    mu_earth = 3.986004418e5  # [km ** 3 / s ** 2]
+    r1 = np.array([2249.171260, 1898.007100, 5639.599193])  # [km]
+    r2 = np.array([1744.495443, -4601.556054, 4043.864391])  # [km]
+    tof = 1618.50  # [s]
+
+    # Unpack problem conditions
+    M, sense, path = case.split("_")
+
+    # Convert proper type
+    M = int(M[-1])
+    sense = True if sense == "prograde" else False
+    path = True if path == "low" else False
+
+    # Solve the problem
+    v1, v2 = solver(mu_earth, r1, r2, tof, M=M, prograde=sense, low_path=path)
+
+    # Expected final results
+    expected_v1, expected_v2 = TABLE_OF_TRANSFERS_I[case]
+
+    # Assert the results
+    assert_allclose(v1, expected_v1, rtol=5e-6)
+    assert_allclose(v2, expected_v2, rtol=5e-6)
+
+
+TABLE_OF_TRANSFERS_II = {
+    "M0_prograde_high": [
+        np.array([2.000652697, 0.387688615, -2.666947760]),  # [km / s]
+        np.array([-3.79246619, -1.77707641, 6.856814395]),  # [km / s]
+    ],
+    "M0_retrograde_high": [
+        np.array([2.96616042, -1.27577231, -0.75545632]),  # [km / s]
+        np.array([5.8437455, -0.20047673, -5.48615883]),  # [km / s]
+    ],
+}
+
+
+@pytest.mark.parametrize("solver", ALL_SOLVERS)
+@pytest.mark.parametrize("case", TABLE_OF_TRANSFERS_II)
+def test_case_from_der_article_II(solver, case):
+    """
+    Example from Astrodynamics 102, by Gim J. Der. see:
+    http://derastrodynamics.com/docs/astrodynamics_102_v2.pdf
+    """
+
+    # Initial conditions
+    mu_earth = 3.986004418e5  # [km ** 3 / s ** 2]
+    r1 = np.array([22592.145603, -1599.915239, -19783.950506])  # [km]
+    r2 = np.array([1922.067697, 4054.157051, -8925.727465])  # [km]
+    tof = 36000  # [s]
+
+    # Unpack problem conditions
+    M, sense, path = case.split("_")
+
+    # Convert proper type
+    M = int(M[-1])
+    sense = True if sense == "prograde" else False
+    path = True if path == "low" else False
+
+    # Solve the problem
+    v1, v2 = solver(mu_earth, r1, r2, tof, M=M, prograde=sense, low_path=path)
+
+    # Expected final results
+    expected_v1, expected_v2 = TABLE_OF_TRANSFERS_II[case]
+
+    # Assert the results
+    assert_allclose(v1, expected_v1, rtol=5e-6)
+    assert_allclose(v2, expected_v2, rtol=5e-6)
diff --git a/tests/test_multirev_solvers.py b/tests/test_multirev_solvers.py
new file mode 100644
index 0000000..c23cdce
--- /dev/null
+++ b/tests/test_multirev_solvers.py
@@ -0,0 +1,74 @@
+""" A collection of  tests only for multi-revolution solvers """
+
+import numpy as np
+import pytest
+from numpy.testing import assert_allclose
+
+from lamberthub import MULTI_REV_SOLVERS
+
+TABLE_OF_TRANSFERS = {
+    "M1_prograde_high": [
+        np.array([0.50335770, 0.61869408, -1.57176904]),  # [km / s]
+        np.array([-4.18334626, -1.13262727, 6.13307091]),  # [km / s]
+    ],
+    "M1_prograde_low": [
+        np.array([-2.45759553, 1.16945801, 0.43161258]),  # [km / s]
+        np.array([-5.53841370, 0.01822220, 5.49641054]),  # [km / s]
+    ],
+    "M1_retrograde_high": [
+        np.array([1.33645655, -0.94654565, 0.30211211]),  # [km / s]
+        np.array([4.93628678, 0.39863416, -5.61593092]),  # [km / s]
+    ],
+    "M1_retrograde_low": [
+        np.array([-1.38861608, -0.47836611, 2.21280154]),  # [km / s]
+        np.array([3.92901545, 1.50871943, -6.52926969]),  # [km / s]
+    ],
+}
+"""
+Directly taken from example 1 from The Superior Lambert Algorithm (Der
+Astrodynamics), by Gim J. Der, see https://amostech.com/TechnicalPapers/2011/Poster/DER.pdf
+"""
+
+
+@pytest.mark.parametrize("solver", MULTI_REV_SOLVERS)
+@pytest.mark.parametrize("case", TABLE_OF_TRANSFERS)
+def test_multirev_case(solver, case):
+
+    # Initial conditions
+    mu_earth = 3.986004418e5  # [km ** 3 / s ** 2]
+    r1 = np.array([22592.145603, -1599.915239, -19783.950506])  # [km]
+    r2 = np.array([1922.067697, 4054.157051, -8925.727465])  # [km]
+    tof = 36000  # [s]
+
+    # Unpack problem conditions
+    M, sense, path = case.split("_")
+
+    # Convert proper type
+    M = int(M[-1])
+    sense = True if sense == "prograde" else False
+    path = True if path == "low" else False
+
+    # Solve the problem
+    v1, v2 = solver(mu_earth, r1, r2, tof, M=M, prograde=sense, low_path=path)
+
+    # Expected final results
+    expected_v1, expected_v2 = TABLE_OF_TRANSFERS[case]
+
+    # Assert the results
+    assert_allclose(v1, expected_v1, rtol=5e-6)
+    assert_allclose(v2, expected_v2, rtol=5e-6)
+
+
+@pytest.mark.parametrize("solver", MULTI_REV_SOLVERS)
+def test_exception_try_lower_M(solver):
+    """ Test that solver does not find any solution for particular input """
+
+    # Initial conditions
+    mu_earth = 3.986004418e5  # [km ** 3 / s ** 2]
+    r1 = np.array([22592.145603, -1599.915239, -19783.950506])  # [km]
+    r2 = np.array([1922.067697, 4054.157051, -8925.727465])  # [km]
+    tof = 5 * 3600  # [s]
+
+    with pytest.raises(ValueError) as excinfo:
+        solver(mu_earth, r1, r2, tof, M=1)
+    assert "ValueError: No feasible solution, try lower M!" in excinfo.exconly()
diff --git a/tests/test_robust_solvers.py b/tests/test_robust_solvers.py
new file mode 100644
index 0000000..7573a73
--- /dev/null
+++ b/tests/test_robust_solvers.py
@@ -0,0 +1,67 @@
+""" A collection of hard tests only for robust solvers """
+
+import numpy as np
+import pytest
+from numpy.testing import assert_allclose
+
+from lamberthub import ROBUST_SOLVERS
+
+TABLE_OF_TRANSFERS = {
+    "M0_prograde_high": [
+        np.array([8.79257809, 0.27867677, 0.02581527]),  # [km / s]
+        np.array([-8.68383320, -0.28592643, -0.03453010]),  # [km / s]
+    ],
+    "M1_prograde_high": [
+        np.array([7.63353091, 0.24582764, 0.02569470]),  # [km / s]
+        np.array([-7.50840227, -0.24335652, -0.02658981]),  # [km / s]
+    ],
+    "M1_prograde_low": [
+        np.array([8.19519089, 2.30595215, 1.75229388]),  # [km / s]
+        np.array([8.07984345, 2.30222567, 1.75189559]),  # [km / s]
+    ],
+    "M2_prograde_high": [
+        np.array([6.51890385, 0.21496104, 0.02618989]),  # [km / s]
+        np.array([-6.37230007, -0.20150975, -0.01832295]),  # [km / s]
+    ],
+    "M2_prograde_low": [
+        np.array([7.00660748, 1.96687296, 1.49423471]),  # [km / s]
+        np.array([6.87133644, 1.96250281, 1.49376762]),  # [km / s]
+    ],
+}
+"""
+Directly taken from example 2 from The Superior Lambert Algorithm (Der
+Astrodynamics), by Gim J. Der, see https://amostech.com/TechnicalPapers/2011/Poster/DER.pdf
+"""
+
+
+@pytest.mark.parametrize("solver", ROBUST_SOLVERS)
+@pytest.mark.parametrize("case", TABLE_OF_TRANSFERS)
+def test_hard_case(solver, case):
+    """
+    Example 2a from The Superior Lambert Algorithm (Der Astrodynamics), by Gim
+    J. Der, see: https://amostech.com/TechnicalPapers/2011/Poster/DER.pdf
+    """
+
+    # Initial conditions
+    mu_earth = 3.986004418e5  # [km ** 3 / s ** 2]
+    r1 = np.array([7231.58074563487, 218.02523761425, 11.79251215952])  # [km]
+    r2 = np.array([7357.06485698842, 253.55724281562, 38.81222241557])  # [km]
+    tof = 12300  # [s]
+
+    # Unpack problem conditions
+    M, sense, path = case.split("_")
+
+    # Convert proper type
+    M = int(M[-1])
+    sense = True if sense == "prograde" else False
+    path = True if path == "low" else False
+
+    # Solve the problem
+    v1, v2 = solver(mu_earth, r1, r2, tof, M=M, prograde=sense, low_path=path)
+
+    # Expected final results
+    expected_v1, expected_v2 = TABLE_OF_TRANSFERS[case]
+
+    # Assert the results
+    assert_allclose(v1, expected_v1, rtol=1e-6)
+    assert_allclose(v2, expected_v2, rtol=1e-6)
diff --git a/tests/tests_utils/test_angles.py b/tests/tests_utils/test_angles.py
new file mode 100644
index 0000000..448df47
--- /dev/null
+++ b/tests/tests_utils/test_angles.py
@@ -0,0 +1,21 @@
+import numpy as np
+import pytest
+from numpy.testing import assert_allclose
+
+from lamberthub.utils.angles import get_transfer_angle
+
+
+@pytest.mark.parametrize("sense", [True, False])
+def test_get_transfer_angle_collinear_vectors_same_sense(sense):
+
+    # Build two vectors with the same direction and sense using scalar
+    # proportion
+    r1 = np.array([1, 0, 0])
+    r2 = r1 if sense is True else -r1
+
+    # Compute the transfer angle
+    dtheta = get_transfer_angle(r1, r2, True)
+    expected_dtheta = 0 if sense is True else np.pi
+
+    # Check that transfer angle is zero
+    assert_allclose(dtheta, expected_dtheta)
diff --git a/tests/tests_utils/test_assertions.py b/tests/tests_utils/test_assertions.py
new file mode 100644
index 0000000..3b50a96
--- /dev/null
+++ b/tests/tests_utils/test_assertions.py
@@ -0,0 +1,48 @@
+import numpy as np
+import pytest
+
+from lamberthub.utils.assertions import (
+    assert_gravitational_parameter_is_positive,
+    assert_number_of_revolutions_not_negative,
+    assert_position_vector_is_valid,
+    assert_position_vectors_are_valid,
+    assert_time_of_flight_is_positive,
+)
+
+
+def test_exception_if_mu_is_negative():
+    with pytest.raises(ValueError) as excinfo:
+        assert_gravitational_parameter_is_positive(-1.00)
+    assert "Gravitational parameter must be positive!" in excinfo.exconly()
+
+
+def test_exception_if_vec_not_three_dimensional():
+    with pytest.raises(ValueError) as excinfo:
+        assert_position_vector_is_valid(np.zeros(4))
+    assert "Vector must be three-dimensional!" in excinfo.exconly()
+
+
+def test_exception_if_vec_is_the_null_vector():
+    with pytest.raises(ValueError) as excinfo:
+        assert_position_vector_is_valid(np.zeros(3))
+    assert "Position vector cannot be the null vector [0,0,0]!" in excinfo.exconly()
+
+
+def test_exception_position_vectors_have_same_components():
+    with pytest.raises(ValueError) as excinfo:
+        assert_position_vectors_are_valid(np.ones(3), np.ones(3))
+    assert "Initial and final position vectors cannot be equal!" in excinfo.exconly()
+
+
+def test_exception_if_time_of_flight_is_negative():
+    with pytest.raises(ValueError) as excinfo:
+        assert_time_of_flight_is_positive(-1.00)
+    assert "Time of flight must be positive!" in excinfo.exconly()
+
+
+def test_exception_if_number_of_revolutions_is_negative():
+    with pytest.raises(ValueError) as excinfo:
+        assert_number_of_revolutions_not_negative(-1)
+    assert (
+        "Number of revolutions must be equal or greater than zero!" in excinfo.exconly()
+    )