diff --git a/Applications/Utils/GeoTools/addDataToRaster.cpp b/Applications/Utils/GeoTools/addDataToRaster.cpp
index 607403bc97233ce188ae4152b909e3c03ec6d7f3..425bfb5d68995ded4d0f5e2ed67ef65c81a03dcc 100644
--- a/Applications/Utils/GeoTools/addDataToRaster.cpp
+++ b/Applications/Utils/GeoTools/addDataToRaster.cpp
@@ -16,9 +16,9 @@
 #endif
 
 #include <algorithm>
-#include <boost/math/constants/constants.hpp>
 #include <cmath>
 #include <memory>
+#include <numbers>
 #include <numeric>
 
 #include "GeoLib/AABB.h"
@@ -45,10 +45,8 @@ double computeSinXSinY(GeoLib::Point const& point, GeoLib::AABB const& aabb)
     auto const aabb_size = aabb.getMaxPoint() - aabb.getMinPoint();
     auto const offset = aabb.getMinPoint();
 
-    return std::sin((point[0] - offset[0]) / aabb_size[0] *
-                    boost::math::double_constants::pi) *
-           std::sin((point[1] - offset[1]) / aabb_size[1] *
-                    boost::math::double_constants::pi);
+    return std::sin((point[0] - offset[0]) / aabb_size[0] * std::numbers::pi) *
+           std::sin((point[1] - offset[1]) / aabb_size[1] * std::numbers::pi);
 }
 
 double computeStepFunction(GeoLib::Point const& point, GeoLib::AABB const& aabb)
diff --git a/GeoLib/Polygon.cpp b/GeoLib/Polygon.cpp
index ab2616f1ab741acfab408b12733557044abbe4cf..b2a686412c240ae3928ef3ab898a8293b1f7845a 100644
--- a/GeoLib/Polygon.cpp
+++ b/GeoLib/Polygon.cpp
@@ -15,7 +15,6 @@
 #include "Polygon.h"
 
 #include <algorithm>
-#include <boost/math/constants/constants.hpp>
 
 #include "AnalyticalGeometry.h"
 #include "BaseLib/quicksort.h"
diff --git a/MaterialLib/MPL/Properties/Exponential.cpp b/MaterialLib/MPL/Properties/Exponential.cpp
index b5490c270c2adfdbb74b336b57cde1ec7406f9b1..5e2f6a10a02a6ebee4ec9e68c1a598f998c971b1 100644
--- a/MaterialLib/MPL/Properties/Exponential.cpp
+++ b/MaterialLib/MPL/Properties/Exponential.cpp
@@ -10,7 +10,6 @@
 
 #include "MaterialLib/MPL/Properties/Exponential.h"
 
-#include <boost/math/special_functions/pow.hpp>
 #include <cmath>
 
 namespace MaterialPropertyLib
diff --git a/MeshToolsLib/MeshQuality/AngleSkewMetric.cpp b/MeshToolsLib/MeshQuality/AngleSkewMetric.cpp
index 29ae02268ee7cc119a02f066211f6e49a6134add..a639aa9e80e4dd012c30499547ad51c76cb88086 100644
--- a/MeshToolsLib/MeshQuality/AngleSkewMetric.cpp
+++ b/MeshToolsLib/MeshQuality/AngleSkewMetric.cpp
@@ -14,14 +14,12 @@
 
 #include "AngleSkewMetric.h"
 
-#include <boost/math/constants/constants.hpp>
 #include <cmath>
+#include <numbers>
 
 #include "MathLib/MathTools.h"
 #include "MeshLib/Node.h"
 
-using namespace boost::math::double_constants;
-
 namespace MeshToolsLib
 {
 namespace
@@ -30,7 +28,7 @@ template <unsigned long N>
 std::tuple<double, double> getMinMaxAngle(
     std::array<MeshLib::Node, N> const& nodes)
 {
-    double min_angle(two_pi);
+    double min_angle(2 * std::numbers::pi);
     double max_angle(0.0);
 
     for (decltype(N) i = 0; i < N; ++i)
@@ -45,11 +43,11 @@ std::tuple<double, double> getMinMaxAngle(
 
 double checkTriangle(MeshLib::Element const& elem)
 {
+    using namespace std::numbers;
     std::array const nodes = {*elem.getNode(0), *elem.getNode(1),
                               *elem.getNode(2)};
     auto const& [min_angle, max_angle] = getMinMaxAngle(nodes);
-    return std::max((max_angle - third_pi) / two_thirds_pi,
-                    (third_pi - min_angle) / third_pi);
+    return std::max((max_angle - pi / 3) / 2, (pi / 3 - min_angle)) * 3 / pi;
 }
 
 double checkQuad(MeshLib::Element const& elem)
@@ -58,8 +56,8 @@ double checkQuad(MeshLib::Element const& elem)
                               *elem.getNode(2), *elem.getNode(3)};
     auto const& [min_angle, max_angle] = getMinMaxAngle(nodes);
 
-    return std::max((max_angle - half_pi) / half_pi,
-                    (half_pi - min_angle) / half_pi);
+    using namespace std::numbers;
+    return std::max((max_angle - pi / 2), (pi / 2 - min_angle)) * 2 / pi;
 }
 
 double checkTetrahedron(MeshLib::Element const& elem)
@@ -77,8 +75,8 @@ double checkTetrahedron(MeshLib::Element const& elem)
     double const min_angle = *std::min_element(min.begin(), min.end());
     double const max_angle = *std::max_element(max.begin(), max.end());
 
-    return std::max((max_angle - third_pi) / two_thirds_pi,
-                    (third_pi - min_angle) / third_pi);
+    using namespace std::numbers;
+    return std::max((max_angle - pi / 3) / 2, (pi / 3 - min_angle)) * 3 / pi;
 }
 
 double checkHexahedron(MeshLib::Element const& elem)
@@ -96,8 +94,8 @@ double checkHexahedron(MeshLib::Element const& elem)
     double const min_angle = *std::min_element(min.begin(), min.end());
     double const max_angle = *std::max_element(max.begin(), max.end());
 
-    return std::max((max_angle - half_pi) / half_pi,
-                    (half_pi - min_angle) / half_pi);
+    using namespace std::numbers;
+    return std::max((max_angle - pi / 2), (pi / 2 - min_angle)) * 2 / pi;
 }
 
 double checkPrism(MeshLib::Element const& elem)
@@ -117,9 +115,10 @@ double checkPrism(MeshLib::Element const& elem)
     auto const min_angle_tri = std::min(min_angle_tri0, min_angle_tri1);
     auto const max_angle_tri = std::max(max_angle_tri0, max_angle_tri1);
 
-    double const tri_criterion(
-        std::max((max_angle_tri - third_pi) / two_thirds_pi,
-                 (third_pi - min_angle_tri) / third_pi));
+    using namespace std::numbers;
+    double const tri_criterion =
+        std::max((max_angle_tri - pi / 3) / 2, (pi / 3 - min_angle_tri)) * 3 /
+        pi;
 
     std::array<double, 3> min;
     std::array<double, 3> max;
@@ -134,8 +133,9 @@ double checkPrism(MeshLib::Element const& elem)
     double const min_angle_quad = *std::min_element(min.begin(), min.end());
     double const max_angle_quad = *std::max_element(max.begin(), max.end());
 
-    double const quad_criterion(std::max((max_angle_quad - half_pi) / half_pi,
-                                         (half_pi - min_angle_quad) / half_pi));
+    using namespace std::numbers;
+    double const quad_criterion =
+        std::max((max_angle_quad - pi / 2), (pi / 2 - min_angle_quad)) * 2 / pi;
 
     return std::min(tri_criterion, quad_criterion);
 }
diff --git a/MeshToolsLib/MeshSurfaceExtraction.cpp b/MeshToolsLib/MeshSurfaceExtraction.cpp
index c606cc2d54f13a8fe715631dd7f215c6bf7d3a37..763835f7d987ce6f75913193acfd754dc929294f 100644
--- a/MeshToolsLib/MeshSurfaceExtraction.cpp
+++ b/MeshToolsLib/MeshSurfaceExtraction.cpp
@@ -14,8 +14,8 @@
 
 #include "MeshSurfaceExtraction.h"
 
-#include <boost/math/constants/constants.hpp>
 #include <memory>
+#include <numbers>
 
 #include "BaseLib/Logging.h"
 #include "MeshLib/Elements/Line.h"
@@ -311,8 +311,7 @@ void MeshSurfaceExtraction::get2DSurfaceElements(
 
     bool const complete_surface = (dir.dot(dir) == 0);
 
-    double const pi(boost::math::constants::pi<double>());
-    double const cos_theta(std::cos(angle * pi / 180.0));
+    double const cos_theta(std::cos(angle * std::numbers::pi / 180.0));
     Eigen::Vector3d const norm_dir(dir.normalized());
 
     for (auto const* elem : all_elements)
diff --git a/NumLib/Fem/FiniteElement/TemplateIsoparametric.h b/NumLib/Fem/FiniteElement/TemplateIsoparametric.h
index da04325acb2196a6888ccc32f169297c3d1754e9..07dc174fab954492dd6c490389832f7b8ac80980 100644
--- a/NumLib/Fem/FiniteElement/TemplateIsoparametric.h
+++ b/NumLib/Fem/FiniteElement/TemplateIsoparametric.h
@@ -12,8 +12,8 @@
 
 #pragma once
 
-#include <boost/math/constants/constants.hpp>
 #include <cassert>
+#include <numbers>
 
 #include "MeshLib/Elements/Element.h"
 #include "MeshLib/Node.h"
@@ -148,7 +148,7 @@ private:
         // located at edge of the unit triangle, it is possible that
         // r becomes zero.
         auto const r = interpolateZerothCoordinate(shape.N);
-        shape.integralMeasure = boost::math::constants::two_pi<double>() * r;
+        shape.integralMeasure = 2. * std::numbers::pi * r;
     }
 
     const MeshLib::Element* _ele;
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1P.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_1P.cpp
index 089623cf2c428c773c9631051bbd9339bf2b209f..bf075fd94fe48c2c4e39ba4863258b7cc88cb5df 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1P.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1P.cpp
@@ -10,7 +10,7 @@
 
 #include "BHE_1P.h"
 
-#include <boost/math/constants/constants.hpp>
+#include <numbers>
 
 #include "FlowAndTemperatureControl.h"
 #include "Physics.h"
@@ -112,7 +112,7 @@ void BHE_1P::updateHeatTransferCoefficients(double const flow_rate)
 std::array<double, BHE_1P::number_of_unknowns> BHE_1P::calcThermalResistances(
     double const Nu)
 {
-    constexpr double pi = boost::math::constants::pi<double>();
+    constexpr double pi = std::numbers::pi;
 
     double const lambda_r = refrigerant.thermal_conductivity;
     double const lambda_g = grout.lambda_g;
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
index e2d10c6847a201010d220bbac9ece63311a6cc77..a6df8711e4539c918d4759f2646a35bde1f35ab6 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_1U.cpp
@@ -10,7 +10,7 @@
 
 #include "BHE_1U.h"
 
-#include <boost/math/constants/constants.hpp>
+#include <numbers>
 
 #include "FlowAndTemperatureControl.h"
 #include "Physics.h"
@@ -173,7 +173,7 @@ void BHE_1U::updateHeatTransferCoefficients(double const flow_rate)
 std::array<double, BHE_1U::number_of_unknowns> BHE_1U::calcThermalResistances(
     double const Nu)
 {
-    constexpr double pi = boost::math::constants::pi<double>();
+    constexpr double pi = std::numbers::pi;
 
     double const lambda_r = refrigerant.thermal_conductivity;
     double const lambda_g = grout.lambda_g;
diff --git a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
index 9034e31d6805c1715dc48f8686615ebb26fd054d..4c269046e0b3b94fa34b3c26a7dbe40017e88cba 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
+++ b/ProcessLib/HeatTransportBHE/BHE/BHE_2U.cpp
@@ -10,7 +10,7 @@
 
 #include "BHE_2U.h"
 
-#include <boost/math/constants/constants.hpp>
+#include <numbers>
 
 #include "FlowAndTemperatureControl.h"
 #include "Physics.h"
@@ -193,7 +193,7 @@ void BHE_2U::updateHeatTransferCoefficients(double const flow_rate)
 std::array<double, BHE_2U::number_of_unknowns> BHE_2U::calcThermalResistances(
     double const Nu)
 {
-    constexpr double pi = boost::math::constants::pi<double>();
+    constexpr double pi = std::numbers::pi;
 
     double const lambda_r = refrigerant.thermal_conductivity;
     double const lambda_g = grout.lambda_g;
diff --git a/ProcessLib/HeatTransportBHE/BHE/BoreholeGeometry.h b/ProcessLib/HeatTransportBHE/BHE/BoreholeGeometry.h
index 0dd6b5d25a253b92151fec2b01f6db5791836610..0b42c799d109c47b1bcb8a89e8b1ec81b85796a6 100644
--- a/ProcessLib/HeatTransportBHE/BHE/BoreholeGeometry.h
+++ b/ProcessLib/HeatTransportBHE/BHE/BoreholeGeometry.h
@@ -10,7 +10,7 @@
 
 #pragma once
 
-#include <boost/math/constants/constants.hpp>
+#include <numbers>
 
 namespace BaseLib
 {
@@ -36,11 +36,7 @@ struct BoreholeGeometry
      */
     double const diameter;
 
-    double area() const
-    {
-        constexpr double pi = boost::math::constants::pi<double>();
-        return pi * diameter * diameter / 4;
-    }
+    double area() const { return std::numbers::pi * diameter * diameter / 4; }
 };
 
 BoreholeGeometry createBoreholeGeometry(BaseLib::ConfigTree const& config);
diff --git a/ProcessLib/HeatTransportBHE/BHE/Pipe.h b/ProcessLib/HeatTransportBHE/BHE/Pipe.h
index 54959fb151ca179863106288c3f7866e8b4f9fa5..9b0c9824d06984bbf4f3396b60ce0cfa983b380d 100644
--- a/ProcessLib/HeatTransportBHE/BHE/Pipe.h
+++ b/ProcessLib/HeatTransportBHE/BHE/Pipe.h
@@ -10,7 +10,8 @@
 
 #pragma once
 
-#include <boost/math/constants/constants.hpp>
+#include <cmath>
+#include <numbers>
 
 namespace BaseLib
 {
@@ -38,19 +39,16 @@ struct Pipe
 
     double wallThermalResistance() const
     {
-        constexpr double pi = boost::math::constants::pi<double>();
-
         double const outside_diameter = outsideDiameter();
 
         return std::log(outside_diameter / diameter) /
-               (2.0 * pi * wall_thermal_conductivity);
+               (2.0 * std::numbers::pi * wall_thermal_conductivity);
     }
 
 private:
     double circleArea(double const diameter) const
     {
-        constexpr double pi = boost::math::constants::pi<double>();
-        return pi * diameter * diameter / 4;
+        return std::numbers::pi * diameter * diameter / 4;
     }
 };
 
diff --git a/ProcessLib/HeatTransportBHE/BHE/ThermalResistancesCoaxial.h b/ProcessLib/HeatTransportBHE/BHE/ThermalResistancesCoaxial.h
index e5e6d3f15a3daddc9ed70711fdeac56bed41af0d..bb755156934772ba6c98568460cdb51bce42d703 100644
--- a/ProcessLib/HeatTransportBHE/BHE/ThermalResistancesCoaxial.h
+++ b/ProcessLib/HeatTransportBHE/BHE/ThermalResistancesCoaxial.h
@@ -10,6 +10,8 @@
 
 #pragma once
 
+#include <numbers>
+
 #include "GroutParameters.h"
 #include "Pipe.h"
 #include "RefrigerantProperties.h"
@@ -47,9 +49,10 @@ inline AdvectiveThermalResistanceCoaxial calculateAdvectiveThermalResistance(
     double const hydraulic_diameter =
         coaxialPipesAnnulusDiameter(inner_pipe, outer_pipe);
 
-    auto advective_thermal_resistance = [&](double Nu, double diameter_ratio) {
-        constexpr double pi = boost::math::constants::pi<double>();
-        return 1.0 / (Nu * fluid.thermal_conductivity * pi) * diameter_ratio;
+    auto advective_thermal_resistance = [&](double Nu, double diameter_ratio)
+    {
+        return 1.0 / (Nu * fluid.thermal_conductivity * std::numbers::pi) *
+               diameter_ratio;
     };
     return {advective_thermal_resistance(Nu_inner_pipe, 1.),
             advective_thermal_resistance(
@@ -70,8 +73,6 @@ calculateGroutAndGroutSoilExchangeThermalResistance(
     Pipe const& outer_pipe, GroutParameters const& grout_parameters,
     double const borehole_diameter)
 {
-    constexpr double pi = boost::math::constants::pi<double>();
-
     double const outer_pipe_outside_diameter = outer_pipe.outsideDiameter();
     double const chi =
         std::log(std::sqrt(borehole_diameter * borehole_diameter +
@@ -81,7 +82,7 @@ calculateGroutAndGroutSoilExchangeThermalResistance(
         std::log(borehole_diameter / outer_pipe_outside_diameter);
     double const R_g =
         std::log(borehole_diameter / outer_pipe_outside_diameter) / 2 /
-        (pi * grout_parameters.lambda_g);
+        (std::numbers::pi * grout_parameters.lambda_g);
     double const conductive_b = chi * R_g;
     double const grout_soil = (1 - chi) * R_g;
     return {conductive_b, grout_soil};
diff --git a/ProcessLib/PhaseField/PhaseFieldFEM-impl.h b/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
index 56b62e26e9688b6e1f59fed60fd4af319dd7a11b..79e6b0a4aa499e58d1c0d9626d070bab3a90cfee 100644
--- a/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
+++ b/ProcessLib/PhaseField/PhaseFieldFEM-impl.h
@@ -10,6 +10,8 @@
  */
 #pragma once
 
+#include <numbers>
+
 #include "NumLib/DOF/DOFTableUtil.h"
 #include "ProcessLib/CoupledSolutionsForStaggeredScheme.h"
 
@@ -238,7 +240,7 @@ void PhaseFieldLocalAssembler<ShapeFunction, DisplacementDim>::
             case PhaseFieldModel::COHESIVE:
             {
                 auto const local_Jac_COHESIVE =
-                    (2.0 / boost::math::double_constants::pi * gc *
+                    (2.0 / std::numbers::pi * gc *
                      (-N.transpose() * N / ls + dNdx.transpose() * dNdx * ls) *
                      w)
                         .eval();
@@ -381,7 +383,7 @@ void PhaseFieldLocalAssembler<ShapeFunction, DisplacementDim>::computeEnergy(
             case PhaseFieldModel::COHESIVE:
             {
                 element_surface_energy +=
-                    gc / boost::math::double_constants::pi *
+                    gc / std::numbers::pi *
                     ((1 - d_ip * d_ip) / ls + (dNdx * d).dot((dNdx * d)) * ls) *
                     w;
                 break;
diff --git a/Tests/GeoLib/TestEquidistantPointGeneration.cpp b/Tests/GeoLib/TestEquidistantPointGeneration.cpp
index 20fabb3e6e0b28c9362089e2fc8cbfbfb1a33c02..49202b9514e5a2df421dcb32e9cc709bca14ea61 100644
--- a/Tests/GeoLib/TestEquidistantPointGeneration.cpp
+++ b/Tests/GeoLib/TestEquidistantPointGeneration.cpp
@@ -10,7 +10,7 @@
 
 #include <gtest/gtest.h>
 
-#include <boost/math/constants/constants.hpp>
+#include <numbers>
 #include <random>
 
 #include "GeoLib/Utils.h"
@@ -22,7 +22,7 @@ struct GeoLibGenerateEquidistantPoints : public testing::Test
         // enable random engine
         std::random_device rd;
         std::mt19937 random_engine_mt19937(rd());
-        constexpr auto pi = boost::math::double_constants::pi;
+        constexpr auto pi = std::numbers::pi;
         std::normal_distribution<> normal_dist_phi(-pi, pi);  // azimuthal angle
         std::normal_distribution<> normal_dist_theta(0, pi);  // polar angle
         // generate point on unit sphere
diff --git a/Tests/GeoLib/TestSurfaceIsPointInSurface.cpp b/Tests/GeoLib/TestSurfaceIsPointInSurface.cpp
index d1b81f509c640cca40da2cf3d75cde6a3ce95611..f617335d878781821dc383fd82a5699fd538df59 100644
--- a/Tests/GeoLib/TestSurfaceIsPointInSurface.cpp
+++ b/Tests/GeoLib/TestSurfaceIsPointInSurface.cpp
@@ -12,8 +12,8 @@
 #include <gtest/gtest.h>
 
 #include <array>
-#include <boost/math/constants/constants.hpp>
 #include <memory>
+#include <numbers>
 #include <random>
 #include <vector>
 
@@ -114,8 +114,7 @@ TEST(GeoLib, SurfaceIsPointInSurface)
                                                            n_steps, f));
 
         // random rotation angles
-        std::normal_distribution<> normal_dist_angles(
-            0, boost::math::double_constants::two_pi);
+        std::normal_distribution<> normal_dist_angles(0, 2. * std::numbers::pi);
         std::array<double, 3> euler_angles = {
             {normal_dist_angles(random_engine_mt19937),
              normal_dist_angles(random_engine_mt19937),
diff --git a/Tests/MaterialLib/TestGasPressureDependentPermeability.cpp b/Tests/MaterialLib/TestGasPressureDependentPermeability.cpp
index 8ea5b381a7f30df66b841dd3af787d3bf6bd6b88..a35c3cc7d7f6934e2c9e5b9d137e7c074615ea5d 100644
--- a/Tests/MaterialLib/TestGasPressureDependentPermeability.cpp
+++ b/Tests/MaterialLib/TestGasPressureDependentPermeability.cpp
@@ -12,7 +12,6 @@
 #include <gtest/gtest.h>
 
 #include <Eigen/Core>
-#include <boost/math/constants/constants.hpp>
 
 #include "MaterialLib/MPL/Medium.h"
 #include "MaterialLib/MPL/Properties/CreateGasPressureDependentPermeability.h"
diff --git a/Tests/MaterialLib/TestStrainDependentPermeability.cpp b/Tests/MaterialLib/TestStrainDependentPermeability.cpp
index afacff3478220b0bb8edc2a43fd8dbdfdaa80ad5..529ec131c67765099cb8097b26d59dc3d8b110e9 100644
--- a/Tests/MaterialLib/TestStrainDependentPermeability.cpp
+++ b/Tests/MaterialLib/TestStrainDependentPermeability.cpp
@@ -12,7 +12,6 @@
 #include <gtest/gtest.h>
 
 #include <Eigen/Core>
-#include <boost/math/constants/constants.hpp>
 
 #include "MaterialLib/MPL/Medium.h"
 #include "MaterialLib/MPL/Properties/CreateStrainDependentPermeability.h"
diff --git a/Tests/MeshLib/TestTetQualityCriteria.cpp b/Tests/MeshLib/TestTetQualityCriteria.cpp
index 4af1e19668fa39ad137ee9bb2bd5bc9bc51b523a..a80bc0999cb732e89eb02af499bf48266ed4b4a2 100644
--- a/Tests/MeshLib/TestTetQualityCriteria.cpp
+++ b/Tests/MeshLib/TestTetQualityCriteria.cpp
@@ -10,9 +10,9 @@
 #include <gtest/gtest.h>
 
 #include <Eigen/Core>
-#include <boost/math/constants/constants.hpp>
 #include <cmath>
 #include <memory>
+#include <numbers>
 #include <numeric>
 #include <random>
 
@@ -148,7 +148,7 @@ std::pair<double, double> getMinMaxAngleFromTriangle(Eigen::Vector3d n0,
 double computeCriterionForTet(Eigen::Vector3d n0, Eigen::Vector3d n1,
                               Eigen::Vector3d n2, Eigen::Vector3d n3)
 {
-    using namespace boost::math::double_constants;
+    using namespace std::numbers;
     auto const alpha0 = getMinMaxAngleFromTriangle(n0, n2, n1);
     auto const alpha1 = getMinMaxAngleFromTriangle(n0, n1, n3);
     auto const alpha2 = getMinMaxAngleFromTriangle(n1, n2, n3);
@@ -157,8 +157,7 @@ double computeCriterionForTet(Eigen::Vector3d n0, Eigen::Vector3d n1,
         std::min({alpha0.first, alpha1.first, alpha2.first, alpha3.first});
     auto const max =
         std::max({alpha0.second, alpha1.second, alpha2.second, alpha3.second});
-    return std::max((max - third_pi) / two_thirds_pi,
-                    (third_pi - min) / third_pi);
+    return std::max((max - pi / 3) * 3 / (2 * pi), (pi / 3 - min) * 3 / pi);
 }
 
 TEST_F(TetElementQuality, EquiAngleSkew)
diff --git a/Tests/MeshLib/TestTriQualityCriteria.cpp b/Tests/MeshLib/TestTriQualityCriteria.cpp
index 245db1d0e5ad87e89de7ec257a76f8c8c3f33ca0..7f1da3c848afbd66b6df5536a31920ee7f80da74 100644
--- a/Tests/MeshLib/TestTriQualityCriteria.cpp
+++ b/Tests/MeshLib/TestTriQualityCriteria.cpp
@@ -9,9 +9,9 @@
 
 #include <gtest/gtest.h>
 
-#include <boost/math/constants/constants.hpp>
 #include <cmath>
 #include <memory>
+#include <numbers>
 #include <numeric>
 #include <random>
 
@@ -103,7 +103,7 @@ TEST_F(TriElementQuality, EdgeRatio)
 
 TEST_F(TriElementQuality, EquiAngleSkew)
 {
-    using namespace boost::math::double_constants;
+    using namespace std::numbers;
     auto const type = MeshLib::MeshQualityType::EQUIANGLESKEW;
     auto const element_quality_vector =
         getElementQualityVectorFromRegularTriMesh(type);
@@ -112,11 +112,11 @@ TEST_F(TriElementQuality, EquiAngleSkew)
                                        std::pow(1.0 / n_subdivisions[1], 2));
     std::array const angles = {
         std::asin((1.0 / n_subdivisions[0]) / hypothenuse),
-        std::asin((1.0 / n_subdivisions[1]) / hypothenuse), half_pi};
+        std::asin((1.0 / n_subdivisions[1]) / hypothenuse), pi / 2};
     auto const& min_max = std::minmax_element(angles.begin(), angles.end());
     auto const expected_value =
-        std::max((*min_max.second - third_pi) / two_thirds_pi,
-                 (third_pi - *min_max.first) / third_pi);
+        std::max((*min_max.second - pi / 3) * 3 / (2 * pi),
+                 (pi / 3 - *min_max.first) * 3 / pi);
 
     for (auto const element_quality : element_quality_vector)
     {
diff --git a/Tests/NumLib/ODEs.h b/Tests/NumLib/ODEs.h
index 2453a4a66a7729013637ded97197caf0b342757d..98d7d61b869e772dc5bf241b436ad5f028c462f5 100644
--- a/Tests/NumLib/ODEs.h
+++ b/Tests/NumLib/ODEs.h
@@ -10,7 +10,7 @@
 
 #pragma once
 
-#include <boost/math/constants/constants.hpp>
+#include <numbers>
 
 #include "MathLib/LinAlg/LinAlg.h"
 #include "MathLib/LinAlg/UnifiedMatrixSetters.h"
@@ -142,8 +142,7 @@ public:
 
 const double ODETraits<ODE1>::t0 = 0.0;
 
-const double ODETraits<ODE1>::t_end =
-    2.0 * 3.141592653589793238462643383279502884197169399375105820974944;
+const double ODETraits<ODE1>::t_end = 2. * std::numbers::pi;
 // ODE 1 end //////////////////////////////////////////////////////
 
 // ODE 2 //////////////////////////////////////////////////////////
@@ -364,6 +363,5 @@ public:
 
 const double ODETraits<ODE3>::t0 = 0.0;
 
-const double ODETraits<ODE3>::t_end =
-    0.5 * boost::math::constants::pi<double>();
+const double ODETraits<ODE3>::t_end = std::numbers::pi / 2.;
 // ODE 3 end //////////////////////////////////////////////////////