diff --git a/benchmark/propagate_lagrangian_benchmark.cpp b/benchmark/propagate_lagrangian_benchmark.cpp
index 9504499b..f2ef4bc6 100644
--- a/benchmark/propagate_lagrangian_benchmark.cpp
+++ b/benchmark/propagate_lagrangian_benchmark.cpp
@@ -44,7 +44,7 @@ void perform_test_speed(
   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> omega_d(0., 2 * pi);
   std::uniform_real_distribution<double> f_d(0, 2 * pi);
   std::uniform_real_distribution<double> tof_d(0.1, 20.);
 
@@ -94,7 +94,7 @@ void perform_test_accuracy(
   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> omega_d(0., 2 * pi);
   std::uniform_real_distribution<double> f_d(0, 2 * pi);
   std::uniform_real_distribution<double> tof_d(0.1, 20.);
 
diff --git a/src/core_astro/par2ic.cpp b/src/core_astro/par2ic.cpp
index 35351ba9..5c7b240b 100644
--- a/src/core_astro/par2ic.cpp
+++ b/src/core_astro/par2ic.cpp
@@ -23,8 +23,6 @@
 
 namespace kep3 {
 
-constexpr double pi4{boost::math::constants::quarter_pi<double>()};
-
 // keplerian osculating elements [a,e,i,W,w,f] -> r,v.
 // The last osculating elements is the true anomaly.
 // The semi-major axis a needs to be positive
diff --git a/src/core_astro/propagate_lagrangian.cpp b/src/core_astro/propagate_lagrangian.cpp
index a0857e4c..e578f3dc 100644
--- a/src/core_astro/propagate_lagrangian.cpp
+++ b/src/core_astro/propagate_lagrangian.cpp
@@ -41,7 +41,7 @@ void propagate_lagrangian(std::array<std::array<double, 3>, 2> &pos_vel_0,
   double sigma0 =
       (r0[0] * v0[0] + r0[1] * v0[1] + r0[2] * v0[2]) / std::sqrt(mu);
 
-  if (a > 0) { // Solve Kepler's equation, elliptical case
+  if (a > 0) { // Solve Kepler's equation in DE, elliptical case
     sqrta = std::sqrt(a);
     double DM = std::sqrt(mu / std::pow(a, 3)) * dt;
     double IG = DM;
@@ -69,7 +69,7 @@ void propagate_lagrangian(std::array<std::array<double, 3>, 2> &pos_vel_0,
         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
+  } else { // Solve Kepler's equation in DH, hyperbolic case
     sqrta = std::sqrt(-a);
     double DN = std::sqrt(-mu / a / a / a) * dt;
     double IG = 0.;
@@ -105,7 +105,7 @@ void propagate_lagrangian(std::array<std::array<double, 3>, 2> &pos_vel_0,
   }
 
   double temp[3] = {r0[0], r0[1], r0[2]};
-  for (int i = 0; i < 3; i++) {
+  for (auto i = 0u; i < 3; i++) {
     r0[i] = F * r0[i] + G * v0[i];
     v0[i] = Ft * temp[i] + Gt * v0[i];
   }
@@ -140,7 +140,7 @@ void propagate_lagrangian_u(std::array<std::array<double, 3>, 2> &pos_vel_0,
   // 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
+  // solve kepler's equation in the universal anomaly DS
   double IG = 0;
   alpha > 0. ? IG = std::sqrt(mu) * dt_copy * std::abs(alpha)
              : IG = 3.; // TODO(darioizzo): initial guess for the universal