Repository: jorgepiloto/lamberthub (language: python)

Issue:
Implement battin1984
💻 **Solver request: battin1984**

This solver by Battin inherits from Gauss but avoiding its singularities for 180 degrees transfer angle. Only drawback with this solver is that it was only developed for the direct transfer scenario. However, [this paper](file:///tmp/mozilla_lobo0/USING_BATTINS_METHOD_TO_OBTAIN_MULTIPLE-REVOLUTIO.pdf) by Shen extends the solver to the multi-revolution case.

📝 **References**

```
@inproceedings{battin1983elegant,
  title={An elegant Lambert algorithm(for spacecraft orbit estimation)},
  author={BATTIN, RH and VAUGHAN, RM},
  booktitle={International Astronautical Federation, International Astronautical Congress, 34 th, Budapest, Hungary},
  volume={10},
  number={15},
  year={1983}
}
```



Discussion / hints:
The continuous fractions used in this algorithm are blocking me... Anyway, Euler to the rescue! The genius developed a formula which for computing these mathematical operations as a series. In fact, I learnt this in the [original work for the paper](https://dspace.mit.edu/bitstream/handle/1721.1/15587/11880750-MIT.pdf?sequence=2&isAllowed=y). Notice that the supervisor of this work is R.H. Battin!
It turns out that this algorithm is also known as Battin-Vaughan algorithm 👀 
It must be said that the effort made by Battin when developing the continued fraction was to get rid of a trigonometric function evaluation. However, I do not think a custom series for computing a particular value of `arctan` is faster than `np.arctan` neither it uses less memory. In addition, even after implementing the original report continued fractions, some errors appeared as these were not covered in the article.

The analysis made by Randal Allen about Battin's Lambert recipe was super useful, see "Tan Root for Battin’s Lambert Recipe" who's DOI is: https://doi.org/10.2514/1.G001013  

Fix patch:
diff --git a/README.md b/README.md
index afad862..d51cde0 100644
--- a/README.md
+++ b/README.md
@@ -36,7 +36,7 @@ print([solver.__name__ for solver in ALL_SOLVERS])
 At this moment, the following algorithms are available:
 
 ```bash
->>> ['gauss1809', 'gooding1990', 'avanzini2008', 'arora2013', 'izzo2015']
+>>> ['gauss1809', 'battin1984', 'gooding1990', 'avanzini2008', 'arora2013', 'izzo2015']
 ```
 
 ## How can I use a solver?
diff --git a/src/lamberthub/__init__.py b/src/lamberthub/__init__.py
index 08620f3..64690be 100644
--- a/src/lamberthub/__init__.py
+++ b/src/lamberthub/__init__.py
@@ -1,6 +1,7 @@
 """ A collection of Lambert's problem solvers """
 
 from lamberthub.ecc_solvers.avanzini import avanzini2008
+from lamberthub.p_solvers.battin import battin1984
 from lamberthub.p_solvers.gauss import gauss1809
 from lamberthub.universal_solvers.arora import arora2013
 from lamberthub.universal_solvers.gooding import gooding1990
@@ -8,10 +9,10 @@
 
 __version__ = "0.2.dev0"
 
-ALL_SOLVERS = [gauss1809, gooding1990, avanzini2008, arora2013, izzo2015]
+ALL_SOLVERS = [gauss1809, battin1984, gooding1990, avanzini2008, arora2013, izzo2015]
 """ A list holding all lamberthub available solvers """
 
-ZERO_REV_SOLVERS = [gooding1990, avanzini2008, arora2013, izzo2015]
+ZERO_REV_SOLVERS = [battin1984, gooding1990, avanzini2008, arora2013, izzo2015]
 """ A list holding all direct-revolution lamberthub solvers """
 
 MULTI_REV_SOLVERS = [gooding1990, izzo2015]
diff --git a/src/lamberthub/p_solvers/battin.py b/src/lamberthub/p_solvers/battin.py
new file mode 100644
index 0000000..10f2936
--- /dev/null
+++ b/src/lamberthub/p_solvers/battin.py
@@ -0,0 +1,281 @@
+"""This module holds all methods devised by R. H. Battin."""
+
+import time
+
+import numpy as np
+
+from lamberthub.utils.angles import get_transfer_angle
+from lamberthub.utils.assertions import assert_parameters_are_valid
+
+
+def battin1984(
+    mu,
+    r1,
+    r2,
+    tof,
+    M=0,
+    prograde=True,
+    low_path=True,
+    maxiter=35,
+    atol=1e-5,
+    rtol=1e-7,
+    full_output=False,
+):
+    r"""
+    Battin's elegant algortihm for solving the Lambert's problem. This algorithm
+    is known to improve Gauss original one by removing the singularity for 180
+    transfer angles and increasing its performance.
+
+    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 and time per iteration are also returned.
+
+    Returns
+    -------
+    v1: numpy.array
+        Initial velocity vector.
+    v2: numpy.array
+        Final velocity vector.
+    numiter: int
+        Number of iterations.
+    tpi: float
+        Time per iteration in seconds.
+
+    Notes
+    -----
+    The algorithm originally devised by Gauss exploits the so-called ratio of
+    sector to triangle area, which is a numerical value related with the orbital
+    parameter. This Algorithm was used to the discovery of the orbit of Ceres by
+    the genius and adopted by many other authors of his time due to its
+    simplicity. However, the Algorithm is found to be singular for transfer
+    angles of 180 degrees and shows a low performance for really small angles.
+
+    References
+    ----------
+    [1] Battin, R. H., & Vaughan, R. M. (1984). An elegant Lambert algorithm.
+    Journal of Guidance, Control, and Dynamics, 7(6), 662-670.
+
+    [2] Battin, R. H. (1999). An introduction to the mathematics and methods of
+    astrodynamics. Aiaa.
+
+    [3] Vaughan, R. M. (1983). An improvement of Gauss' method for solving
+    Lambert's problem (Doctoral dissertation, Massachusetts Institute of
+    Technology).
+
+    """
+
+    # Check that input parameters are safe
+    assert_parameters_are_valid(mu, r1, r2, tof, M)
+
+    # Retrieve the fundamental geometry of the problem
+    r1_norm, r2_norm, c_norm = [np.linalg.norm(vec) for vec in [r1, r2, r2 - r1]]
+    semiperimeter = (r1_norm + r2_norm + c_norm) / 2
+    dtheta = get_transfer_angle(r1, r2, prograde)
+
+    # Compute the auxiliary lambda parameter. The variable _lambda is used here
+    # as lambda is a reserved keyword in Python for declaring inline functions.
+    # This expression is equivalent to equation (1) from [1]. Equation (7.122a)
+    # applies here. The sign is selected according to the value of the transfer
+    # angle, being positive if lower than 180 degrees or negative if greater.
+    # Therefore, Battin's lambda variable is the so-called transfer angle
+    # parameter.
+    _lambda = np.sqrt(semiperimeter * (semiperimeter - c_norm)) / semiperimeter
+    _lambda = np.abs(_lambda) if dtheta < np.pi else -np.abs(_lambda)
+
+    # Solve for the new auxiliary variables developed by Battin. Equations (30)
+    # and (31) from the report [1] apply here.
+    ll = ((1 - _lambda) / (1 + _lambda)) ** 2
+    m = (8 * mu * tof ** 2) / (semiperimeter ** 3 * (1 + _lambda) ** 6)
+
+    # The non-dimensional transfer time is computed using equation (54) from
+    # Battin's report [1]. The non-dimensional time for the parabolic case is
+    # also computed here using equation appearing in page 340 of Battin's book
+    # [2]. The author did not provide this in the original paper [1], which is a
+    # critical errata, as it is used in the initial guess computation.
+    T = np.sqrt(8 * mu / semiperimeter ** 3) * tof
+    T_p = (4 / 3) * (1 - _lambda ** 3)
+
+    # The initial guess procedure is set according to piece-wise no-numbered
+    # equation appearing in the last page of the report [1]. This relation is
+    # also shown in page 340 from book [2], where it is claimed that if T is
+    # greater than T_p, the transfer orbit is expected to be elliptic type,
+    # otherwise parabolic or hyperbolic are assumed.
+    x0 = ll if T > T_p else 0
+
+    # The iterative procedure starts
+    tic = time.perf_counter()
+    for numiter in range(1, maxiter + 1):
+
+        # Evaluate the xi function at particular value of x. This is equation
+        # (53) from original report [1] or expression (7.121) from the book [2].
+        # The function follows the series technique for an easy implementation
+        # from the computational point of view.
+        xi = _xi_at_x(x0)
+
+        # Evaluate the h coefficients given by equations (47) and (48) from the
+        # report [1] or the relations (7.111) and (7.112) from the book [2]. The
+        # denominator is common for both h1 and h2 coefficients, so it is
+        # to avoid computing it twice
+        h_denominator = (1 + 2 * x0 + ll) * (4 * x0 + xi * (3 + x0))
+        h1 = ((ll + x0) ** 2 * (1 + 3 * x0 + xi)) / h_denominator
+        h2 = (m * (x0 - ll + xi)) / (h_denominator)
+
+        # For a better convergence of the solution, Battin introduced in [1] the
+        # equations (55), (56) and (58), being this last one a continued
+        # fraction again.
+        B = (27 * h2) / (4 * (1 + h1) ** 3)
+        u = -B / (2 * np.sqrt(1 + B) + 1)
+        K = _K_at_u(u)
+
+        # Compute the desired positive root via an explicit equation defined in
+        # [1] under the number (57).
+        y = ((1 + h1) / 3) * (2 + np.sqrt(B + 1) / (1 - 2 * u * K))
+
+        # Compute the new value for x using equation (49) from original report
+        # [1], who's alternative from the book [2] is expression (7.113)
+        x = np.sqrt(((1 - ll) / 2) ** 2 + m / y ** 2) - (1 + ll) / 2
+
+        # Check if the new computed value lies within desired tolerance. If so,
+        # stop the iteration procedure, otherwise update the initial guess and
+        # keep performing the computation.
+        if np.abs(x - x0) <= atol:
+            tac = time.perf_counter()
+            tpi = (tac - tic) / numiter
+            break
+        else:
+            x0 = x
+    else:
+        raise ValueError("Exceeded maximum number of iterations!")
+
+    r11 = (1 + _lambda) ** 2 / (4 * tof * _lambda)
+    s11 = y * (1 + x)
+    t11 = (m * semiperimeter * (1 + _lambda) ** 2) / s11
+
+    # Compute the initial and final velocity vectors
+    v1 = -r11 * (s11 * (r1 - r2) - t11 * r1 / r1_norm)
+    v2 = -r11 * (s11 * (r1 - r2) + t11 * r2 / r2_norm)
+
+    return (v1, v2, numiter, tpi) if full_output is True else (v1, v2)
+
+
+def _xi_at_x(x, levels=125):
+    """Evaluates the xi function at a particular value of x.
+
+    Parameters
+    ----------
+    x: float
+        The independent variable.
+
+    Returns
+    -------
+    xi: float
+        The value of the xi function.
+
+    Notes
+    -----
+    This is equation (53) from original report [1]. However, the method
+    presented in [3] which makes use of a series is the one used here as it is
+    much simpler to implement.
+
+    """
+
+    # Compute the value of eta, given by equation (52) in [1].
+    eta = x / (np.sqrt(1 + x) + 1) ** 2
+
+    # Initial values
+    delta, u, sigma, m1, = (
+        0,
+        1,
+        1,
+        1,
+    )
+
+    while np.abs(u) > 1e-18 and m1 <= levels:
+        m1 += 1
+        gamma = (m1 + 3) ** 2 / (4 * (m1 + 3) ** 2 - 1)
+        delta = 1 / (1 + gamma * eta * delta)
+        u = u * (delta - 1)
+        sigma = sigma + u
+
+    xi = 8 * (np.sqrt(1 + x) + 1) / (3 + 1 / (5 + eta + (9 * eta / 7) * sigma))
+    return xi
+
+
+def _K_at_u(u, levels=1000):
+    """Evaluates the K function at a particular value of u.
+
+    Parameters
+    ----------
+    u: float
+        Battin's auxiliary variable.
+
+    Returns
+    -------
+    K: float
+        The value of the K function.
+
+    Notes
+    -----
+    This is equation (58) from original report [1]. However, the method
+    presented in [3] which makes use of a series is the one used here as it is
+    much simpler to implement.
+
+    """
+
+    # Initial values
+    delta, u0, sigma, n1, = (
+        1,
+        1,
+        1,
+        0,
+    )
+
+    while np.abs(u0) > 1e-18 and n1 <= levels:
+        if n1 == 0:
+            gamma = 4 / 27
+            delta = 1 / (1 - gamma * u * delta)
+            u0 = u0 * (delta - 1)
+            sigma = sigma + u0
+        else:
+            for val in [1, 2]:
+                if val == 1:
+                    gamma = (
+                        2
+                        * (3 * n1 + 1)
+                        * (6 * n1 - 1)
+                        / (9 * (4 * n1 - 1) * (4 * n1 + 1))
+                    )
+                else:
+                    gamma = (
+                        2
+                        * (3 * n1 + 2)
+                        * (6 * n1 + 1)
+                        / (9 * (4 * n1 + 1) * (4 * n1 + 3))
+                    )
+                delta = 1 / (1 - gamma * u * delta)
+                u0 = u0 * (delta - 1)
+                sigma = sigma + u0
+
+        n1 = n1 + 1
+
+    K = (sigma / 3) ** 2
+    return K


Test patch:
diff --git a/tests/test_all_solvers.py b/tests/test_all_solvers.py
index 5aadef5..4b18492 100644
--- a/tests/test_all_solvers.py
+++ b/tests/test_all_solvers.py
@@ -13,6 +13,9 @@
 """A list holding all solvers which present a decent accuracy. Most of the
 solvers from the NON_ROBUST_SOLVERS macro were developed in the old days."""
 
+(ATOL, RTOL) = (0.02, 0.001)
+"""The absolute and relative tolerances for the tests."""
+
 
 @pytest.mark.parametrize("solver", ALL_SOLVERS)
 def test_case_from_vallado_book(solver):
@@ -35,8 +38,8 @@ def test_case_from_vallado_book(solver):
     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)
+    assert_allclose(v1, expected_v1, atol=ATOL, rtol=RTOL)
+    assert_allclose(v2, expected_v2, atol=ATOL, rtol=ATOL)
 
 
 @pytest.mark.parametrize("solver", ALL_SOLVERS)
@@ -60,8 +63,8 @@ def test_case_from_curtiss_book(solver):
     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)
+    assert_allclose(v1, expected_v1, atol=ATOL, rtol=RTOL)
+    assert_allclose(v2, expected_v2, atol=ATOL, rtol=RTOL)
 
 
 @pytest.mark.parametrize("solver", ALL_SOLVERS_LAMBERTHUB)
@@ -83,7 +86,7 @@ def test_case_from_battin_book(solver):
     expected_v1 = np.array([-9.303603251, 3.018641330, 1.536362143])  # [AU / year]
 
     # Assert the results
-    assert_allclose(v1, expected_v1, rtol=1e-5)
+    assert_allclose(v1, expected_v1, atol=ATOL, rtol=RTOL)
 
 
 @pytest.mark.parametrize("solver", ALL_SOLVERS)
@@ -110,8 +113,8 @@ def test_case_from_GMAT_hyperbolic_prograde(solver):
     expected_v2 = np.array([-3.6379, 4.4932, 1.7735])  # [km / s]
 
     # Assert the results
-    assert_allclose(v1, expected_v1, atol=5e-3, rtol=1e-4)
-    assert_allclose(v2, expected_v2, atol=5e-3, rtol=1e-4)
+    assert_allclose(v1, expected_v1, atol=ATOL, rtol=RTOL)
+    assert_allclose(v2, expected_v2, atol=ATOL, rtol=RTOL)
 
 
 @pytest.mark.parametrize("solver", ALL_SOLVERS)
@@ -138,8 +141,8 @@ def test_case_from_GMAT_hyperbolic_retrograde(solver):
     expected_v2 = np.array([-4.3016, -3.4314, -2.5467])  # [km / s]
 
     # Assert the results
-    assert_allclose(v1, expected_v1, atol=5e-3, rtol=1e-4)
-    assert_allclose(v2, expected_v2, atol=5e-3, rtol=1e-4)
+    assert_allclose(v1, expected_v1, atol=ATOL, rtol=RTOL)
+    assert_allclose(v2, expected_v2, atol=ATOL, rtol=RTOL)
 
 
 TABLE_OF_TRANSFERS_I = {
@@ -183,8 +186,8 @@ def test_case_from_der_article_I(solver, case):
     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)
+    assert_allclose(v1, expected_v1, atol=ATOL, rtol=RTOL)
+    assert_allclose(v2, expected_v2, atol=ATOL, rtol=RTOL)
 
 
 TABLE_OF_TRANSFERS_II = {
@@ -228,5 +231,5 @@ def test_case_from_der_article_II(solver, case):
     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)
+    assert_allclose(v1, expected_v1, atol=ATOL, rtol=RTOL)
+    assert_allclose(v2, expected_v2, atol=ATOL, rtol=RTOL)
diff --git a/tests/tests_p_solvers/test_battin.py b/tests/tests_p_solvers/test_battin.py
new file mode 100644
index 0000000..7b5a059
--- /dev/null
+++ b/tests/tests_p_solvers/test_battin.py
@@ -0,0 +1,19 @@
+import numpy as np
+import pytest
+
+from lamberthub import battin1984
+
+
+def test_exceeded_maximum_number_of_iterations():
+    # Initial conditions
+    mu_earth = 3.986004418e5  # [km ** 3 / s ** 2]
+    r1 = np.array([1.0, 0.0, 0.0])  # [km]
+    r2 = np.array([0.0, 1.0, 0.0])  # [km]
+    tof = 1000  # [s]
+
+    # Solving the problem with only two iteration so an error is raised
+    with pytest.raises(ValueError) as excinfo:
+        v1, v2 = battin1984(
+            mu_earth, r1, r2, tof, maxiter=1, prograde=True, low_path=True
+        )
+    assert "Exceeded maximum number of iterations!" in excinfo.exconly()