diff --git a/include/kep3/core_astro/constants.hpp b/include/kep3/core_astro/constants.hpp
index a5d8cae6..f400df77 100644
--- a/include/kep3/core_astro/constants.hpp
+++ b/include/kep3/core_astro/constants.hpp
@@ -29,7 +29,7 @@ inline constexpr double CAVENDISH = 6.67430e-11;             // Cavendish consta
 inline constexpr double MU_SUN = 1.32712440041279419e20;     // Sun's gravitational parameter (m^3/s^2 kg) - DE440
 inline constexpr double MU_EARTH = 398600435507000.0;        // Earth's gravitational parameter (m^3/s^2 kg) - DE440
 inline constexpr double MU_MOON =  4902800118000.0;          // Moon's gravitational parameter (m^3/s^2 kg) - DE440
-inline constexpr double EARTH_VELOCITY = 29784.691831696804; // Earth's velocity. (m/s)
+inline constexpr double EARTH_VELOCITY = 29784.69183430911;  // Earth's velocity. (m/s)
 inline constexpr double EARTH_J2 = 1.08262668E-03;
 inline constexpr double EARTH_RADIUS = 6378137; // Earth's radius (m)
 inline constexpr double DEG2RAD = (pi / 180.0);
diff --git a/test/leg_sims_flanagan_test.cpp b/test/leg_sims_flanagan_test.cpp
index 04718309..298fd869 100644
--- a/test/leg_sims_flanagan_test.cpp
+++ b/test/leg_sims_flanagan_test.cpp
@@ -246,20 +246,23 @@ TEST_CASE("compute_mismatch_constraints_test")
 
 TEST_CASE("mismatch_constraints_test2")
 {
-    // We test the correctness of the compute_mismatch_constraints computations against a ground truth (computed with a different program)
+    // The ground truth were computed with these values of the astro constants, hence we cannot use here the pykep ones.
+    double AU_OLD = 149597870700.0;
+    double EV_OLD = 29784.691831696804;
+    double MU_OLD = 1.32712440018e20;
+    // We test the correctness of the compute_mismatch_constraints computations against a ground truth (computed with a
+    // different program)
     std::array<std::array<double, 3>, 2> rvs{
-        {{1 * kep3::AU, 0.1 * kep3::AU, -0.1 * kep3::AU},
-         {0.2 * kep3::EARTH_VELOCITY, 1 * kep3::EARTH_VELOCITY, -0.2 * kep3::EARTH_VELOCITY}}};
+        {{1 * AU_OLD, 0.1 * AU_OLD, -0.1 * AU_OLD}, {0.2 * EV_OLD, 1 * EV_OLD, -0.2 * EV_OLD}}};
 
     std::array<std::array<double, 3>, 2> rvf{
-        {{1.2 * kep3::AU, -0.1 * kep3::AU, 0.1 * kep3::AU},
-         {-0.2 * kep3::EARTH_VELOCITY, 1.023 * kep3::EARTH_VELOCITY, -0.44 * kep3::EARTH_VELOCITY}}};
+        {{1.2 * AU_OLD, -0.1 * AU_OLD, 0.1 * AU_OLD}, {-0.2 * EV_OLD, 1.023 * EV_OLD, -0.44 * EV_OLD}}};
 
     double ms = 1500.;
     double mf = 1300.;
     std::vector<double> throttles
         = {0.10, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2, 0.21, 0.22, 0.23, 0.24};
-    kep3::leg::sims_flanagan sf(rvs, ms, throttles, rvf, mf, 324.0 * kep3::DAY2SEC, 0.12, 100, kep3::MU_SUN, 0.6);
+    kep3::leg::sims_flanagan sf(rvs, ms, throttles, rvf, mf, 324.0 * kep3::DAY2SEC, 0.12, 100, MU_OLD, 0.6);
     auto retval = sf.compute_mismatch_constraints();
     std::vector<double> ground_truth
         = {-1.9701274809621304e+11, 4.6965044246848071e+11, -1.5007523306033661e+11, -2.9975151466948650e+04,
@@ -292,7 +295,8 @@ TEST_CASE("grad_test")
     auto grad_a = udp.gradient(x);
     auto xgrad = xt::adapt(grad, {1u + 7u + 5u, 17u});
     auto xgrad_a = xt::adapt(grad_a, {1u + 7u + 5u, 17u});
-    REQUIRE(xt::linalg::norm(xgrad - xgrad_a) < 1e-8); // With the high fidelity gradient this is still the best we can achieve
+    REQUIRE(xt::linalg::norm(xgrad - xgrad_a)
+            < 1e-8); // With the high fidelity gradient this is still the best we can achieve
 }
 
 TEST_CASE("serialization_test")
diff --git a/test/stark_problem_test.cpp b/test/stark_problem_test.cpp
index 6e026a75..52cec9ca 100644
--- a/test/stark_problem_test.cpp
+++ b/test/stark_problem_test.cpp
@@ -34,29 +34,39 @@ TEST_CASE("construct")
 
 TEST_CASE("propagate")
 {
-    auto s = kep3::stark_problem(kep3::MU_SUN, 3000 * kep3::G0, 1e-16);
+    // The ground truth were computed with these values of the astro constants, hence we cannot use here the pykep ones.
+    double AU_OLD = 149597870700.0;
+    double EV_OLD = 29784.691831696804;
+    double MU_OLD = 1.32712440018e20;
+
+    auto s = kep3::stark_problem(MU_OLD, 3000 * kep3::G0, 1e-16);
     // Mass is zero
-    REQUIRE_THROWS_AS(s.propagate({kep3::AU, 1000, 01000, 0.123, kep3::EARTH_VELOCITY, 10.1234, 0.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC), std::domain_error);
+    REQUIRE_THROWS_AS(s.propagate({AU_OLD, 1000, 01000, 0.123, EV_OLD, 10.1234, 0.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC), std::domain_error);
     // Radius is zero
-    REQUIRE_THROWS_AS(s.propagate({0., 0., 0., 0.123, kep3::EARTH_VELOCITY, 10.1234, 0.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC), std::domain_error);
+    REQUIRE_THROWS_AS(s.propagate({0., 0., 0., 0.123, EV_OLD, 10.1234, 0.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC), std::domain_error);
     // Thrust is zero
-    REQUIRE_NOTHROW(s.propagate({kep3::AU, 1000, 01000, 0.123, kep3::EARTH_VELOCITY, 10.1234, 1000.}, {0., 0., 0.}, 300. * kep3::DAY2SEC));
+    REQUIRE_NOTHROW(s.propagate({AU_OLD, 1000, 01000, 0.123, EV_OLD, 10.1234, 1000.}, {0., 0., 0.}, 300. * kep3::DAY2SEC));
 
-    auto res = s.propagate({kep3::AU, 1000, 1000, 0.123, kep3::EARTH_VELOCITY, 10.1234, 1000.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC);
+    auto res = s.propagate({AU_OLD, 1000, 1000, 0.123, EV_OLD, 10.1234, 1000.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC);
     std::array<double, 7> const ground_truth{86447318382.73862, -109196478678.49327, 9774837.693084948, 25535.89383875655, 18933.302734531168, -47.00078590023982, 954.2200884785473};
     REQUIRE(L_infinity_norm_rel(res, ground_truth) <=1e-13);
 }
 TEST_CASE("propagate_var")
 {
-    auto s = kep3::stark_problem(kep3::MU_SUN, 3000 * kep3::G0, 1e-16);
+    // The ground truth were computed with these values of the astro constants, hence we cannot use here the pykep ones.
+    double AU_OLD = 149597870700.0;
+    double EV_OLD = 29784.691831696804;
+    double MU_OLD = 1.32712440018e20;
+
+    auto s = kep3::stark_problem(MU_OLD, 3000 * kep3::G0, 1e-16);
     // Mass is zero
-    REQUIRE_THROWS_AS(s.propagate_var({kep3::AU, 1000, 01000, 0.123, kep3::EARTH_VELOCITY, 10.1234, 0.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC), std::domain_error);
+    REQUIRE_THROWS_AS(s.propagate_var({AU_OLD, 1000, 01000, 0.123, EV_OLD, 10.1234, 0.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC), std::domain_error);
     // Radius is zero
-    REQUIRE_THROWS_AS(s.propagate_var({0., 0., 0., 0.123, kep3::EARTH_VELOCITY, 10.1234, 0.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC), std::domain_error);
+    REQUIRE_THROWS_AS(s.propagate_var({0., 0., 0., 0.123, EV_OLD, 10.1234, 0.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC), std::domain_error);
     // Thrust is zero
-    REQUIRE_THROWS_AS(s.propagate_var({kep3::AU, 1000, 01000, 0.123, kep3::EARTH_VELOCITY, 10.1234, 1000.}, {0., 0., 0.}, 300. * kep3::DAY2SEC), std::domain_error);
+    REQUIRE_THROWS_AS(s.propagate_var({AU_OLD, 1000, 01000, 0.123, EV_OLD, 10.1234, 1000.}, {0., 0., 0.}, 300. * kep3::DAY2SEC), std::domain_error);
 
-    auto res = s.propagate_var({kep3::AU, 1000, 1000, 0.123, kep3::EARTH_VELOCITY, 10.1234, 1000.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC);
+    auto res = s.propagate_var({AU_OLD, 1000, 1000, 0.123, EV_OLD, 10.1234, 1000.}, {0.05, 0.01, 0.01}, 300. * kep3::DAY2SEC);
     std::array<double, 7> const ground_truth{86447318382.73862, -109196478678.49327, 9774837.693084948, 25535.89383875655, 18933.302734531168, -47.00078590023982, 954.2200884785473};
     REQUIRE(L_infinity_norm_rel(std::get<0>(res), ground_truth) <=1e-13);
 }
diff --git a/test/ta_stark_test.cpp b/test/ta_stark_test.cpp
index a7d09dc6..003aa046 100644
--- a/test/ta_stark_test.cpp
+++ b/test/ta_stark_test.cpp
@@ -23,8 +23,8 @@ using heyoka::taylor_adaptive;
 using heyoka::taylor_outcome;
 
 using kep3::ta::get_ta_stark;
-using kep3::ta::get_ta_stark_var;
 using kep3::ta::get_ta_stark_cache_dim;
+using kep3::ta::get_ta_stark_var;
 using kep3::ta::get_ta_stark_var_cache_dim;
 
 using kep3_tests::L_infinity_norm_rel;
@@ -37,9 +37,9 @@ TEST_CASE("caches")
     auto ta_cached = get_ta_stark(1e-16);
     REQUIRE(get_ta_stark_cache_dim() == 1u);
     ta_cached = get_ta_stark(1e-16);
-    REQUIRE(get_ta_stark_cache_dim() ==1u);
+    REQUIRE(get_ta_stark_cache_dim() == 1u);
     ta_cached = get_ta_stark(1e-8);
-    REQUIRE(get_ta_stark_cache_dim() ==2u);
+    REQUIRE(get_ta_stark_cache_dim() == 2u);
 
     // The variational integrator.
     REQUIRE(get_ta_stark_var_cache_dim() == 0u);
@@ -48,7 +48,7 @@ TEST_CASE("caches")
     ta_cached = get_ta_stark_var(1e-16);
     REQUIRE(get_ta_stark_var_cache_dim() == 1u);
     ta_cached = get_ta_stark_var(1e-8);
-    REQUIRE(get_ta_stark_var_cache_dim() == 2u);    
+    REQUIRE(get_ta_stark_var_cache_dim() == 2u);
 }
 
 TEST_CASE("dynamics")
@@ -72,12 +72,15 @@ TEST_CASE("dynamics")
     }
     {
         // We test a generic case.
+        // The ground truth were computed with these values of the astro constants, hence we cannot use here the pykep
+        // ones.
+        double MU_OLD = 1.32712440018e20;
         taylor_adaptive<double> ta(ta_cached); // making a copy as to be able to modify the object.
         ta.set_time(0.);
         std::vector<double> ic{164557657760.1, 179517444829.19998, 47871318621.12, 32763.159550000004,
                                23827.7524,     9531.10096,         1200.};
         std::copy(ic.begin(), ic.end(), ta.get_state_data());
-        std::vector<double> pars{kep3::MU_SUN, 3000.*kep3::G0, 0.05, 0., 0.032};
+        std::vector<double> pars{MU_OLD, 3000. * kep3::G0, 0.05, 0., 0.032};
         std::copy(pars.begin(), pars.end(), ta.get_pars_data());
         auto out = ta.propagate_until(3888000.0);
         std::vector<double> const ground_truth
@@ -107,24 +110,28 @@ TEST_CASE("variational_dynamics")
         std::copy(pars.begin(), pars.end(), ta.get_pars_data());
         auto out = ta.propagate_until(2. * kep3::pi);
         REQUIRE(std::get<0>(out) == taylor_outcome::time_limit);
-        REQUIRE(L_infinity_norm_rel(std::vector<double>(ta.get_state().begin(), ta.get_state().begin() + 7), ic) <= 1e-13);
+        REQUIRE(L_infinity_norm_rel(std::vector<double>(ta.get_state().begin(), ta.get_state().begin() + 7), ic)
+                <= 1e-13);
     }
     {
         // We test a generic case.
+        // The ground truth were computed with these values of the astro constants, hence we cannot use here the pykep
+        // ones.
+        double MU_OLD = 1.32712440018e20;
         taylor_adaptive<double> ta(ta_cached); // making a copy as to be able to modify the object.
         ta.set_time(0.);
         std::vector<double> ic{164557657760.1, 179517444829.19998, 47871318621.12, 32763.159550000004,
                                23827.7524,     9531.10096,         1200.};
         std::copy(ic.begin(), ic.end(), ta.get_state_data());
-        std::vector<double> pars{kep3::MU_SUN, 3000.*kep3::G0, 0.05, 0., 0.032};
+        std::vector<double> pars{MU_OLD, 3000. * kep3::G0, 0.05, 0., 0.032};
         std::copy(pars.begin(), pars.end(), ta.get_pars_data());
         auto out = ta.propagate_until(3888000.0);
         std::vector<double> const ground_truth
             = {284296823432.0578,  263961690798.0665, 82814214381.94377, 29341.50292902231,
                20219.294034700008, 8592.028822618351, 1192.1548315009777};
         REQUIRE(std::get<0>(out) == taylor_outcome::time_limit);
-        REQUIRE(kep3_tests::L_infinity_norm_rel(std::vector<double>(ta.get_state().begin(), ta.get_state().begin() + 7), ground_truth)
+        REQUIRE(kep3_tests::L_infinity_norm_rel(std::vector<double>(ta.get_state().begin(), ta.get_state().begin() + 7),
+                                                ground_truth)
                 <= 1e-13);
     }
 }
-