diff --git a/benchmark/CMakeLists.txt b/benchmark/CMakeLists.txt
index 1a00de3b..1530f755 100644
--- a/benchmark/CMakeLists.txt
+++ b/benchmark/CMakeLists.txt
@@ -17,3 +17,4 @@ function(ADD_kep3_BENCHMARK arg1)
 endfunction()
 
 ADD_kep3_BENCHMARK(convert_anomalies_benchmark)
+ADD_kep3_BENCHMARK(propagate_lagrangian_benchmark)
diff --git a/benchmark/convert_anomalies_benchmark.cpp b/benchmark/convert_anomalies_benchmark.cpp
index 11488d1a..d242302a 100644
--- a/benchmark/convert_anomalies_benchmark.cpp
+++ b/benchmark/convert_anomalies_benchmark.cpp
@@ -15,8 +15,11 @@
 
 #include <kep3/core_astro/convert_anomalies.hpp>
 
-using namespace kep3;
-using namespace std::chrono;
+using kep3::e2m;
+using kep3::m2e;
+using std::chrono::duration_cast;
+using std::chrono::high_resolution_clock;
+using std::chrono::microseconds;
 
 // In this benchmark we test the speed and accuracy of the Kepler's equation
 // solvers
@@ -25,9 +28,10 @@ void perform_test_speed(double min_ecc, double max_ecc, unsigned N) {
   //
   // Engines
   //
+  // NOLINTNEXTLINE(cert-msc32-c, cert-msc51-cpp)
   std::mt19937 rng_engine(122012203u);
   //
-  // Distribtuions
+  // Distributions
   //
   std::uniform_real_distribution<double> ecc_d(min_ecc, max_ecc);
   std::uniform_real_distribution<double> M_d(-100, 100.);
@@ -51,16 +55,17 @@ void perform_test_speed(double min_ecc, double max_ecc, unsigned N) {
   }
   auto stop = high_resolution_clock::now();
   auto duration = duration_cast<microseconds>(stop - start);
-  fmt::print("{:.3f}s\n", duration.count() / 1e6);
+  fmt::print("{:.3f}s\n", (static_cast<double>(duration.count()) / 1e6));
 }
 
 void perform_test_accuracy(double min_ecc, double max_ecc, unsigned N) {
   //
   // Engines
   //
+  // NOLINTNEXTLINE(cert-msc32-c, cert-msc51-cpp)
   std::mt19937 rng_engine(122012203u);
   //
-  // Distribtuions
+  // Distributions
   //
   std::uniform_real_distribution<double> ecc_d(min_ecc, max_ecc);
   std::uniform_real_distribution<double> M_d(-100, 100.);
@@ -84,13 +89,14 @@ void perform_test_accuracy(double min_ecc, double max_ecc, unsigned N) {
   }
   auto max_it = max_element(std::begin(err), std::end(err));
   auto min_it = min_element(std::begin(err), std::end(err));
-  auto avg = std::accumulate(err.begin(), err.end(), 0.0) / err.size();
+  auto avg = std::accumulate(err.begin(), err.end(), 0.0) /
+             static_cast<double>(err.size());
   fmt::print("{:.3e} avg, {:.3e} min, {:.3e} max\n", avg, *min_it, *max_it);
 }
 
 int main() {
   unsigned seed = 7898935u;
-  fmt::print("\nComputes speed different eccentricity ranges:\n");
+  fmt::print("\nComputes speed at different eccentricity ranges:\n");
   perform_test_speed(0, 0.5, 1000000);
   perform_test_speed(0.5, 0.9, 1000000);
   perform_test_speed(0.9, 0.99, 1000000);
diff --git a/benchmark/propagate_lagrangian_benchmark.cpp b/benchmark/propagate_lagrangian_benchmark.cpp
new file mode 100644
index 00000000..8a9d7812
--- /dev/null
+++ b/benchmark/propagate_lagrangian_benchmark.cpp
@@ -0,0 +1,166 @@
+// Copyright 2023, 2024 Dario Izzo (dario.izzo@gmail.com), Francesco Biscani
+// (bluescarni@gmail.com)
+//
+// This file is part of the kep3 library.
+//
+// This Source Code Form is subject to the terms of the Mozilla
+// Public License v. 2.0. If a copy of the MPL was not distributed
+// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+#include "kep3/core_astro/par2ic.hpp"
+#include "kep3/core_astro/propagate_lagrangian.hpp"
+#include <chrono>
+#include <iostream>
+#include <random>
+
+#include <fmt/core.h>
+#include <fmt/ranges.h>
+
+#include <kep3/core_astro/convert_anomalies.hpp>
+#include <xtensor/xtensor.hpp>
+
+using std::chrono::duration_cast;
+using std::chrono::high_resolution_clock;
+using std::chrono::microseconds;
+
+constexpr double pi{boost::math::constants::pi<double>()};
+
+// In this benchmark we test the speed and accuracy of the Lagrangian
+// propagation solvers
+
+void perform_test_speed(
+    double min_ecc, double max_ecc, unsigned N,
+    const std::function<void(std::array<std::array<double, 3>, 2> &, double,
+                             double)> &propagate) {
+  //
+  // Engines
+  //
+  // NOLINTNEXTLINE(cert-msc32-c, cert-msc51-cpp)
+  std::mt19937 rng_engine(122012203u);
+  //
+  // Distributions
+  //
+  std::uniform_real_distribution<double> sma_d(0.5, 20.);
+  std::uniform_real_distribution<double> ecc_d(min_ecc, max_ecc);
+  std::uniform_real_distribution<double> incl_d(0., pi);
+  std::uniform_real_distribution<double> Omega_d(0, 2 * pi);
+  std::uniform_real_distribution<double> omega_d(0., pi);
+  std::uniform_real_distribution<double> f_d(0, 2 * pi);
+  std::uniform_real_distribution<double> tof_d(0.1, 20.);
+
+  // We generate the random dataset
+  std::vector<std::array<std::array<double, 3>, 2>> pos_vels(N);
+  std::vector<double> tofs(N);
+  for (auto i = 0u; i < N; ++i) {
+    auto ecc = ecc_d(rng_engine);
+    auto sma = sma_d(rng_engine);
+    ecc > 1. ? sma = -sma : sma = sma;
+    double f = pi;
+    while (std::cos(f) < -1. / ecc && sma < 0.) {
+      f = f_d(rng_engine);
+    }
+    pos_vels[i] = kep3::par2ic({sma, ecc, incl_d(rng_engine),
+                                Omega_d(rng_engine), omega_d(rng_engine), f},
+                               1.);
+    tofs[i] = tof_d(rng_engine);
+  }
+
+  // We log progress
+  fmt::print("{:.2f} min_ecc, {:.2f} max_ecc, on {} data points: ", min_ecc,
+             max_ecc, N);
+
+  auto start = high_resolution_clock::now();
+  for (auto i = 0u; i < N; ++i) {
+    propagate(pos_vels[i], tofs[i], 1.);
+  }
+  auto stop = high_resolution_clock::now();
+  auto duration = duration_cast<microseconds>(stop - start);
+  fmt::print("{:.3f}s\n", (static_cast<double>(duration.count()) / 1e6));
+}
+
+void perform_test_accuracy(
+    double min_ecc, double max_ecc, unsigned N,
+    const std::function<void(std::array<std::array<double, 3>, 2> &, double,
+                             double)> &propagate) {
+  //
+  // Engines
+  //
+  // NOLINTNEXTLINE(cert-msc32-c, cert-msc51-cpp)
+  std::mt19937 rng_engine(122012203u);
+  //
+  // Distributions
+  //
+  std::uniform_real_distribution<double> sma_d(0.5, 20.);
+  std::uniform_real_distribution<double> ecc_d(min_ecc, max_ecc);
+  std::uniform_real_distribution<double> incl_d(0., pi);
+  std::uniform_real_distribution<double> Omega_d(0, 2 * pi);
+  std::uniform_real_distribution<double> omega_d(0., pi);
+  std::uniform_real_distribution<double> f_d(0, 2 * pi);
+  std::uniform_real_distribution<double> tof_d(0.1, 20.);
+
+  // We generate the random dataset
+  std::vector<std::array<std::array<double, 3>, 2>> pos_vels(N);
+  std::vector<double> tofs(N);
+  for (auto i = 0u; i < N; ++i) {
+    double f = pi;
+    double ecc = 10.;
+    double sma = -1.;
+    while (std::cos(f) < -1. / ecc && sma < 0.) {
+      ecc = ecc_d(rng_engine);
+      sma = sma_d(rng_engine);
+      ecc > 1. ? sma = -sma : sma = sma;
+      f = f_d(rng_engine);
+    }
+
+    pos_vels[i] = {
+        sma, ecc, incl_d(rng_engine), Omega_d(rng_engine), omega_d(rng_engine),
+        f};
+    tofs[i] = tof_d(rng_engine);
+  }
+  // We log progress
+  fmt::print("{:.2f} min_ecc, {:.2f} max_ecc, on {} data points: ", min_ecc,
+             max_ecc, N);
+  std::vector<double> err(N);
+  auto pos_vels_old(pos_vels);
+  for (auto i = 0u; i < N; ++i) {
+    propagate(pos_vels[i], tofs[i], 1.);
+    propagate(pos_vels[i], -tofs[i], 1.);
+    err[i] = std::abs(pos_vels[i][0][0] - pos_vels_old[i][0][0]) +
+             std::abs(pos_vels[i][0][1] - pos_vels_old[i][0][1]) +
+             std::abs(pos_vels[i][0][2] - pos_vels_old[i][0][2]) +
+             std::abs(pos_vels[i][1][0] - pos_vels_old[i][1][0]) +
+             std::abs(pos_vels[i][1][1] - pos_vels_old[i][1][1]) +
+             std::abs(pos_vels[i][1][2] - pos_vels_old[i][1][2]);
+  }
+  auto max_it = max_element(std::begin(err), std::end(err));
+  auto min_it = min_element(std::begin(err), std::end(err));
+  auto avg = std::accumulate(err.begin(), err.end(), 0.0) /
+             static_cast<double>(err.size()) / 6.;
+  fmt::print("{:.3e} avg, {:.3e} min, {:.3e} max\n", avg, *min_it, *max_it);
+}
+
+int main() {
+  fmt::print("\nComputes speed at different eccentricity ranges:\n");
+  perform_test_speed(0, 0.5, 1000000, &kep3::propagate_lagrangian);
+  perform_test_speed(0.5, 0.9, 1000000, &kep3::propagate_lagrangian);
+  perform_test_speed(0.9, 0.99, 1000000, &kep3::propagate_lagrangian);
+  perform_test_speed(1.1, 10., 1000000, &kep3::propagate_lagrangian);
+
+  fmt::print("\nComputes error at different eccentricity ranges:\n");
+  perform_test_accuracy(0, 0.5, 100000, &kep3::propagate_lagrangian);
+  perform_test_accuracy(0.5, 0.9, 100000, &kep3::propagate_lagrangian);
+  perform_test_accuracy(0.9, 0.99, 100000, &kep3::propagate_lagrangian);
+
+  fmt::print("\nComputes speed at different eccentricity ranges [Universal "
+             "Variable]:\n");
+  perform_test_speed(0, 0.5, 1000000, &kep3::propagate_lagrangian_u);
+  perform_test_speed(0.5, 0.9, 1000000, &kep3::propagate_lagrangian_u);
+  perform_test_speed(0.9, 0.99, 1000000, &kep3::propagate_lagrangian_u);
+  perform_test_speed(1.1, 10., 1000000, &kep3::propagate_lagrangian_u);
+
+  fmt::print("\nComputes error at different eccentricity ranges [Universal "
+             "Variable]:\n");
+  perform_test_accuracy(0, 0.5, 100000, &kep3::propagate_lagrangian_u);
+  perform_test_accuracy(0.5, 0.9, 100000, &kep3::propagate_lagrangian_u);
+  perform_test_accuracy(0.9, 0.99, 100000, &kep3::propagate_lagrangian_u);
+}
\ No newline at end of file
diff --git a/kep3_devel.yml b/kep3_devel.yml
index 1d41011e..c7871bf5 100644
--- a/kep3_devel.yml
+++ b/kep3_devel.yml
@@ -11,4 +11,5 @@ dependencies:
   - heyoka >=0.21.0
   - spdlog
   - xtensor
+  - xtensor-blas
 
diff --git a/src/core_astro/propagate_lagrangian.cpp b/src/core_astro/propagate_lagrangian.cpp
index b14f0ae5..e96ed538 100644
--- a/src/core_astro/propagate_lagrangian.cpp
+++ b/src/core_astro/propagate_lagrangian.cpp
@@ -16,6 +16,7 @@
 #include <kep3/core_astro/kepler_equations.hpp>
 #include <kep3/core_astro/propagate_lagrangian.hpp>
 #include <kep3/core_astro/special_functions.hpp>
+#include <stdexcept>
 
 namespace kep3 {
 
@@ -28,7 +29,6 @@ constexpr double pi{boost::math::constants::pi<double>()};
  * numerical technique. All units systems can be used, as long
  * as the input parameters are all expressed in the same system.
  */
-
 void propagate_lagrangian(std::array<std::array<double, 3>, 2> &pos_vel_0,
                           const double dt, const double mu) {
   auto &[r0, v0] = pos_vel_0;
@@ -44,50 +44,64 @@ void propagate_lagrangian(std::array<std::array<double, 3>, 2> &pos_vel_0,
 
   if (a > 0) { // Solve Kepler's equation, elliptical case
     sqrta = std::sqrt(a);
-    double DM = std::sqrt(mu / pow(a, 3)) * dt;
+    double DM = std::sqrt(mu / std::pow(a, 3)) * dt;
 
     // Solve Kepler Equation for ellipses in DE (eccentric anomaly difference)
     const int digits = std::numeric_limits<double>::digits;
+    std::uintmax_t max_iter = 100u;
     double DE = boost::math::tools::halley_iterate(
         [DM, sigma0, sqrta, a, R](double DE) {
           return std::make_tuple(kepDE(DE, DM, sigma0, sqrta, a, R),
                                  d_kepDE(DE, sigma0, sqrta, a, R),
                                  dd_kepDE(DE, sigma0, sqrta, a, R));
         },
-        DM, DM - pi, DM + pi, digits);
-
-    double r = a + (R - a) * cos(DE) + sigma0 * sqrta * sin(DE);
+        DM, DM - pi, DM + pi, digits, max_iter);
+    if (max_iter == 100u) {
+      throw std::domain_error(
+          "Maximum number of iterations exceeded when solving Kepler's "
+          "equation for the eccentric anomaly in propagate_lagrangian.");
+    }
+    double r = a + (R - a) * std::cos(DE) + sigma0 * sqrta * std::sin(DE);
 
     // Lagrange coefficients
-    F = 1 - a / R * (1 - cos(DE));
-    G = a * sigma0 / sqrt(mu) * (1 - cos(DE)) + R * sqrt(a / mu) * sin(DE);
-    Ft = -sqrt(mu * a) / (r * R) * sin(DE);
-    Gt = 1 - a / r * (1 - cos(DE));
+    F = 1 - a / R * (1 - std::cos(DE));
+    G = a * sigma0 / std::sqrt(mu) * (1 - std::cos(DE)) +
+        R * std::sqrt(a / mu) * std::sin(DE);
+    Ft = -std::sqrt(mu * a) / (r * R) * std::sin(DE);
+    Gt = 1 - a / r * (1 - std::cos(DE));
   } else { // Solve Kepler's equation, hyperbolic case
-    sqrta = sqrt(-a);
-    double DN = sqrt(-mu / pow(a, 3)) * dt;
+    sqrta = std::sqrt(-a);
+    double DN = std::sqrt(-mu / a / a / a) * dt;
     double IG = 0.;
-    dt > 0 ? IG = 1 : IG = -1; // TODO(darioizzo): find a better initial guess.
-                               // I tried with 0 and D (both have numercial
-                               // problems and result in exceptions)
+    dt > 0. ? IG = 3.
+            : IG = -3.; // TODO(darioizzo): find a better initial guess.
+                        // I tried with 0 and DN (both have numercial
+                        // problems and result in exceptions)
 
     // Solve Kepler Equation for ellipses in DH (hyperbolic anomaly difference)
     const int digits = std::numeric_limits<double>::digits;
+    std::uintmax_t max_iter = 100u;
     double DH = boost::math::tools::halley_iterate(
         [DN, sigma0, sqrta, a, R](double DH) {
           return std::make_tuple(kepDH(DH, DN, sigma0, sqrta, a, R),
                                  d_kepDH(DH, sigma0, sqrta, a, R),
                                  dd_kepDH(DH, sigma0, sqrta, a, R));
         },
-        IG, IG - pi, IG + pi, digits);
+        IG, IG - pi, IG + pi, digits, max_iter);
+    if (max_iter == 100u) {
+      throw std::domain_error(
+          "Maximum number of iterations exceeded when solving Kepler's "
+          "equation for the hyperbolic anomaly in propagate_lagrangian.");
+    }
 
-    double r = a + (R - a) * cosh(DH) + sigma0 * sqrta * sinh(DH);
+    double r = a + (R - a) * std::cosh(DH) + sigma0 * sqrta * std::sinh(DH);
 
     // Lagrange coefficients
-    F = 1 - a / R * (1 - cosh(DH));
-    G = a * sigma0 / sqrt(mu) * (1 - cosh(DH)) + R * sqrt(-a / mu) * sinh(DH);
-    Ft = -sqrt(-mu * a) / (r * R) * sinh(DH);
-    Gt = 1 - a / r * (1 - cosh(DH));
+    F = 1 - a / R * (1 - std::cosh(DH));
+    G = a * sigma0 / std::sqrt(mu) * (1 - std::cosh(DH)) +
+        R * std::sqrt(-a / mu) * std::sinh(DH);
+    Ft = -std::sqrt(-mu * a) / (r * R) * std::sinh(DH);
+    Gt = 1 - a / r * (1 - std::cosh(DH));
   }
 
   double temp[3] = {r0[0], r0[1], r0[2]};
@@ -97,14 +111,14 @@ void propagate_lagrangian(std::array<std::array<double, 3>, 2> &pos_vel_0,
   }
 }
 
-/// Universial Variales version
+/// Universial Variables version
 /**
  * This function has the same prototype as kep3::propagate_lgrangian, but
  * internally makes use of universal variables formulation for the Lagrange
  * Coefficients.
  */
 void propagate_lagrangian_u(std::array<std::array<double, 3>, 2> &pos_vel_0,
-                            const double dt, const double mu) {
+                            const double dt, const double mu) { // NOLINT
   // If time is negative we need to invert time and velocities. Unlike the other
   // formulation of the propagate lagrangian we cannot rely on negative times to
   // automatically mean back-propagation
@@ -119,46 +133,54 @@ void propagate_lagrangian_u(std::array<std::array<double, 3>, 2> &pos_vel_0,
   }
 
   double F = 0., G = 0., Ft = 0., Gt = 0.;
-  double R0 = sqrt(r0[0] * r0[0] + r0[1] * r0[1] + r0[2] * r0[2]);
-  double V0 = sqrt(v0[0] * v0[0] + v0[1] * v0[1] + v0[2] * v0[2]);
+  double R0 = std::sqrt(r0[0] * r0[0] + r0[1] * r0[1] + r0[2] * r0[2]);
+  double V0 = std::sqrt(v0[0] * v0[0] + v0[1] * v0[1] + v0[2] * v0[2]);
   // the reciprocal of the semi-major axis
   double alpha = 2 / R0 - V0 * V0 / mu;
   // initial radial velocity
   double VR0 = (r0[0] * v0[0] + r0[1] * v0[1] + r0[2] * v0[2]) / R0;
 
   // solve kepler's equation in universal variables
-  double IG = 1;
-  alpha > 0 ? IG = sqrt(mu) * dt_copy * std::abs(alpha)
-            : IG = 1; // initial guess for the universal anomaly. For
-                      // hyperbolas it is 1.... can be better?
+  double IG = 0;
+  alpha > 0. ? IG = std::sqrt(mu) * dt_copy * std::abs(alpha)
+             : IG = 3.; // TODO(darioizzo): initial guess for the universal anomaly. For
+                        // hyperbolas is 3 .... can be better?
 
   // Solve Kepler Equation in DS (univrsal anomaly difference)
   const int digits = std::numeric_limits<double>::digits;
-  double DS = boost::math::tools::halley_iterate(
+  std::uintmax_t max_iter = 100u;
+  // NOTE: Halley iterate here seems to introduce issues in corner cases
+  // and be slower anyway.
+  double DS = boost::math::tools::newton_raphson_iterate(
       [dt_copy, R0, VR0, alpha, mu](double DS) {
         return std::make_tuple(kepDS(DS, dt_copy, R0, VR0, alpha, mu),
-                               d_kepDS(DS, R0, VR0, alpha, mu),
-                               dd_kepDS(DS, R0, VR0, alpha, mu));
+                               d_kepDS(DS, R0, VR0, alpha, mu));
       },
-      IG, IG - pi, IG + pi, digits);
-
+      IG, IG - 2 * pi, IG + 2 * pi, digits,
+      max_iter); // limiting the IG error within
+                 // only pi will not work.
+  if (max_iter == 100u) {
+    throw std::domain_error(
+        "Maximum number of iterations exceeded when solving Kepler's "
+        "equation for the universal anomaly in propagate_lagrangian_u.");
+  }
   // evaluate the lagrangian coefficients F and G
   double S = stumpff_s(alpha * DS * DS);
   double C = stumpff_c(alpha * DS * DS);
   //
   double z = alpha * DS * DS;
   F = 1 - DS * DS / R0 * C;
-  G = dt_copy - 1 / sqrt(mu) * DS * DS * DS * S;
+  G = dt_copy - 1 / std::sqrt(mu) * DS * DS * DS * S;
 
   double r0_copy[3] = {r0[0], r0[1], r0[2]};
   // compute the final position
   r0[0] = F * r0[0] + G * v0[0];
   r0[1] = F * r0[1] + G * v0[1];
   r0[2] = F * r0[2] + G * v0[2];
-  double RF = sqrt(r0[0] * r0[0] + r0[1] * r0[1] + r0[2] * r0[2]);
+  double RF = std::sqrt(r0[0] * r0[0] + r0[1] * r0[1] + r0[2] * r0[2]);
 
   // compute the lagrangian coefficients Ft, Gt
-  Ft = sqrt(mu) / RF / R0 * (z * S - 1) * DS;
+  Ft = std::sqrt(mu) / RF / R0 * (z * S - 1) * DS;
   Gt = 1 - DS * DS / RF * C;
 
   // compute the final velocity
diff --git a/test/convert_anomalies_test.cpp b/test/convert_anomalies_test.cpp
index 0af282a4..ce210f83 100644
--- a/test/convert_anomalies_test.cpp
+++ b/test/convert_anomalies_test.cpp
@@ -30,7 +30,7 @@ TEST_CASE("m2e") {
   //
   std::mt19937 rng_engine(rd());
   //
-  // Distribtuions
+  // Distributions
   //
   std::uniform_real_distribution<double> ecc_difficult_d(0.9, 0.99);
   std::uniform_real_distribution<double> ecc_easy_d(0., 0.9);
@@ -107,7 +107,7 @@ TEST_CASE("zeta2e") {
   //
   std::mt19937 rng_engine(rd());
   //
-  // Distribtuions
+  // Distributions
   //
   std::uniform_real_distribution<double> ecc_difficult_d(1.01, 1.1);
   std::uniform_real_distribution<double> ecc_easy_d(2., 100.);
diff --git a/test/propagate_lagrangian_test.cpp b/test/propagate_lagrangian_test.cpp
index 90c9c2a6..c74501eb 100644
--- a/test/propagate_lagrangian_test.cpp
+++ b/test/propagate_lagrangian_test.cpp
@@ -17,12 +17,12 @@
 
 #include <kep3/core_astro/ic2par.hpp>
 #include <kep3/core_astro/kepler_equations.hpp>
+#include <kep3/core_astro/par2ic.hpp>
 #include <kep3/core_astro/propagate_lagrangian.hpp>
 
 #include "catch.hpp"
 
 using Catch::Matchers::WithinAbs;
-using Catch::Matchers::WithinRel;
 using kep3::propagate_lagrangian;
 using kep3::propagate_lagrangian_u;
 
@@ -84,36 +84,115 @@ void test_propagate_lagrangian(
   // Engines
   // NOLINTNEXTLINE(cert-msc32-c, cert-msc51-cpp)
   std::mt19937 rng_engine(12201203u);
-  // Distribution for the initial Cartesian state
-  std::uniform_real_distribution<double> r_d(-2., 2.);
-  std::uniform_real_distribution<double> v_d(-2., 2.);
-  std::uniform_real_distribution<double> time(0.1, 20.);
-  // Testing on N random calls
-  for (auto i = 0u; i < N; ++i) {
-    std::array<double, 3> pos = {r_d(rng_engine), r_d(rng_engine),
-                                 r_d(rng_engine)};
-    std::array<double, 3> vel = {v_d(rng_engine), v_d(rng_engine),
-                                 v_d(rng_engine)};
-    std::array<std::array<double, 3>, 2> pos_vel = {pos, vel};
-    auto par_before = kep3::ic2par(pos_vel, 1.0);
-    // We filter out cases of little significance (too close to singularity)
-    if (std::abs(par_before[0]) > 0.5 && std::abs(par_before[0]) < 10. &&
-        std::abs(1 - par_before[1]) > 1e-1) {
-      propagate(pos_vel, time(rng_engine), 1.);
+  { // Random istribution of the initial Cartesian state (will mainly produce
+    // hyperbolas)
+    std::uniform_real_distribution<double> r_d(-2., 2.);
+    std::uniform_real_distribution<double> v_d(-2., 2.);
+    std::uniform_real_distribution<double> time(0.1, 20.);
+    // Testing on N random calls
+    for (auto i = 0u; i < N; ++i) {
+      std::array<double, 3> pos = {r_d(rng_engine), r_d(rng_engine),
+                                   r_d(rng_engine)};
+      std::array<double, 3> vel = {v_d(rng_engine), v_d(rng_engine),
+                                   v_d(rng_engine)};
+      std::array<std::array<double, 3>, 2> pos_vel = {pos, vel};
+      auto par_before = kep3::ic2par(pos_vel, 1.0);
+      // We filter out cases of little significance (too close to singularity)
+      if (std::abs(par_before[0]) > 0.5 && std::abs(par_before[0]) < 10. &&
+          std::abs(1 - par_before[1]) > 1e-1) {
+        propagate(pos_vel, time(rng_engine), 1.);
+        auto par_after = kep3::ic2par(pos_vel, 1.0);
+        if (std::isfinite(par_before[0]) && std::isfinite(par_after[0])) {
+          REQUIRE_THAT(par_before[0], WithinAbs(par_after[0], 1e-8));
+          REQUIRE_THAT(par_before[1], WithinAbs(par_after[1], 1e-8));
+          REQUIRE_THAT(par_before[2], WithinAbs(par_after[2], 1e-8));
+          REQUIRE_THAT(par_before[3], WithinAbs(par_after[3], 1e-8));
+          REQUIRE_THAT(par_before[4], WithinAbs(par_after[4], 1e-8));
+        }
+      }
+    }
+  }
+  { // Targeting Ellipses
+    std::uniform_real_distribution<double> sma_d(1.1, 100.);
+    std::uniform_real_distribution<double> ecc_d(0, 0.99);
+    std::uniform_real_distribution<double> incl_d(0., pi);
+    std::uniform_real_distribution<double> Omega_d(0, 2 * pi);
+    std::uniform_real_distribution<double> omega_d(0., pi);
+    std::uniform_real_distribution<double> f_d(0, 2 * pi);
+    std::uniform_real_distribution<double> time_d(0.1, 20.);
+
+    // Testing on N random calls on ellipses
+    for (auto i = 0u; i < N; ++i) {
+      double sma = sma_d(rng_engine);
+      double ecc = ecc_d(rng_engine);
+      double incl = incl_d(rng_engine);
+      double Omega = Omega_d(rng_engine);
+      double omega = omega_d(rng_engine);
+      double f = f_d(rng_engine);
+
+      std::array<double, 6> par_before = {sma, ecc, incl, Omega, omega, f};
+      auto pos_vel = kep3::par2ic(par_before, 1.);
+
+      propagate(pos_vel, time_d(rng_engine), 1.);
       auto par_after = kep3::ic2par(pos_vel, 1.0);
-      fmt::print("\n\n{}", par_before);
-      fmt::print("\n{}", par_after);
-
-      REQUIRE_THAT(par_before[0], WithinAbs(par_after[0], 1e-11));
-      REQUIRE_THAT(par_before[1], WithinAbs(par_after[1], 1e-11));
-      REQUIRE_THAT(par_before[2], WithinAbs(par_after[2], 1e-11));
-      REQUIRE_THAT(par_before[3], WithinAbs(par_after[3], 1e-11));
-      REQUIRE_THAT(par_before[4], WithinAbs(par_after[4], 1e-11));
+      REQUIRE_THAT(par_before[0], WithinAbs(par_after[0], 1e-8));
+      REQUIRE_THAT(par_before[1], WithinAbs(par_after[1], 1e-8));
+      REQUIRE_THAT(par_before[2], WithinAbs(par_after[2], 1e-8));
+      REQUIRE_THAT(par_before[3], WithinAbs(par_after[3], 1e-8));
+      REQUIRE_THAT(par_before[4], WithinAbs(par_after[4], 1e-8));
+    }
+  }
+
+  { // Targeting Hyperbolas
+    std::uniform_real_distribution<double> sma_d(1.1, 100.);
+    std::uniform_real_distribution<double> ecc_d(2., 20.);
+    std::uniform_real_distribution<double> incl_d(0., pi);
+    std::uniform_real_distribution<double> Omega_d(0, 2 * pi);
+    std::uniform_real_distribution<double> omega_d(0., pi);
+    std::uniform_real_distribution<double> f_d(0, 2 * pi);
+    std::uniform_real_distribution<double> time_d(0.1, 20.);
+    // Testing on N random calls on hyperbolas
+    for (auto i = 0u; i < N; ++i) {
+      double sma = sma_d(rng_engine);
+      double ecc = ecc_d(rng_engine);
+      double incl = incl_d(rng_engine);
+      double Omega = Omega_d(rng_engine);
+      double omega = omega_d(rng_engine);
+      double f = f_d(rng_engine);
+      if (std::cos(f) > -1 / ecc) {
+        std::array<double, 6> par_before = {-sma, ecc, incl, Omega, omega, f};
+        auto pos_vel = kep3::par2ic(par_before, 1.);
+
+        propagate(pos_vel, time_d(rng_engine), 1.);
+        auto par_after = kep3::ic2par(pos_vel, 1.0);
+        REQUIRE_THAT(par_before[0], WithinAbs(par_after[0], 1e-8));
+        REQUIRE_THAT(par_before[1], WithinAbs(par_after[1], 1e-8));
+        REQUIRE_THAT(par_before[2], WithinAbs(par_after[2], 1e-8));
+        REQUIRE_THAT(par_before[3], WithinAbs(par_after[3], 1e-8));
+        REQUIRE_THAT(par_before[4], WithinAbs(par_after[4], 1e-8));
+      }
     }
   }
 }
 
 TEST_CASE("propagate_lagrangian") {
-  //test_propagate_lagrangian(&propagate_lagrangian, 10000u);
-  // test_propagate_lagrangian(&propagate_lagrangian_u, 10000u);
+  // We test both Normal and Universal variables version with the same data.
+  test_propagate_lagrangian(&propagate_lagrangian, 10000u);
+  test_propagate_lagrangian(&propagate_lagrangian_u, 10000u);
+}
+
+TEST_CASE("infinite loop") {
+  std::array<std::array<double, 3>, 2> pos_vel = {
+      3.2479950146598147, 4.866100102242875,  0.8564969484971678,
+      0.3140399734911721, 0.5042257639093218, 0.09475180867356801};
+  double tof = 6.023574175415248;
+  kep3::propagate_lagrangian(pos_vel, -tof, 1.);
+}
+
+TEST_CASE("extreme_orbit") {
+  std::array<std::array<double, 3>, 2> pos_vel = {
+      0.8086322075411211,   -1.3297145067523164, -2.4969299661382585,
+      -0.02869546877795607, 0.05765808202641542, -0.03999826575867087};
+  double tof = 4.454030166101634;
+  kep3::propagate_lagrangian_u(pos_vel, tof, 1.);
 }
\ No newline at end of file