diff --git a/docs/users_guide/framework_basics/quantity_arithmetics.md b/docs/users_guide/framework_basics/quantity_arithmetics.md
index ceb35262..4795baba 100644
--- a/docs/users_guide/framework_basics/quantity_arithmetics.md
+++ b/docs/users_guide/framework_basics/quantity_arithmetics.md
@@ -343,7 +343,7 @@ Among others, we can find there the following:
 - `inverse()`,
 - `hypot()`,
 - `sin()`, `cos()`, `tan()`,
-- `asin()`, `acos()`, `atan()`.
+- `asin()`, `acos()`, `atan()`, `atan2()`.
 
 In the library, we can also find _mp-units/random.h_ header file with all the pseudo-random number
 generators working on quantity types.
diff --git a/src/systems/include/mp-units/systems/angular/math.h b/src/systems/include/mp-units/systems/angular/math.h
index f2954f05..bbfaa033 100644
--- a/src/systems/include/mp-units/systems/angular/math.h
+++ b/src/systems/include/mp-units/systems/angular/math.h
@@ -119,4 +119,24 @@ template<ReferenceOf<dimensionless> auto R, typename Rep>
     return quantity{atan(q.numerical_value_in(one)), radian};
 }
 
+template<auto R1, typename Rep1, auto R2, typename Rep2>
+  requires requires(Rep1 v1, Rep2 v2) {
+    common_reference(R1, R2);
+    requires requires { atan2(v1, v2); } || requires { std::atan2(v1, v2); };
+  }
+[[nodiscard]] inline QuantityOf<angle> auto atan2(const quantity<R1, Rep1>& y, const quantity<R2, Rep2>& x) noexcept
+{
+  constexpr auto ref = common_reference(R1, R2);
+  constexpr auto unit = get_unit(ref);
+  using std::atan2;
+  if constexpr (!treat_as_floating_point<Rep1> || !treat_as_floating_point<Rep2>) {
+    // check what is the return type when called with the integral value
+    using rep = decltype(atan2(y.force_numerical_value_in(unit), x.force_numerical_value_in(unit)));
+    // use this type ahead of calling the function to prevent narrowing if a unit conversion is needed
+    return quantity{atan2(value_cast<rep>(y).numerical_value_in(unit), value_cast<rep>(x).numerical_value_in(unit)),
+                    radian};
+  } else
+    return quantity{atan2(y.numerical_value_in(unit), x.numerical_value_in(unit)), radian};
+}
+
 }  // namespace mp_units::angular
diff --git a/src/systems/include/mp-units/systems/si/math.h b/src/systems/include/mp-units/systems/si/math.h
index a2f5971f..eafddcbf 100644
--- a/src/systems/include/mp-units/systems/si/math.h
+++ b/src/systems/include/mp-units/systems/si/math.h
@@ -120,4 +120,25 @@ template<ReferenceOf<dimensionless> auto R, typename Rep>
     return quantity{atan(q.numerical_value_in(one)), radian};
 }
 
+template<auto R1, typename Rep1, auto R2, typename Rep2>
+  requires requires(Rep1 v1, Rep2 v2) {
+    common_reference(R1, R2);
+    requires requires { atan2(v1, v2); } || requires { std::atan2(v1, v2); };
+  }
+[[nodiscard]] inline QuantityOf<isq::angular_measure> auto atan2(const quantity<R1, Rep1>& y,
+                                                                 const quantity<R2, Rep2>& x) noexcept
+{
+  constexpr auto ref = common_reference(R1, R2);
+  constexpr auto unit = get_unit(ref);
+  using std::atan2;
+  if constexpr (!treat_as_floating_point<Rep1> || !treat_as_floating_point<Rep2>) {
+    // check what is the return type when called with the integral value
+    using rep = decltype(atan2(y.force_numerical_value_in(unit), x.force_numerical_value_in(unit)));
+    // use this type ahead of calling the function to prevent narrowing if a unit conversion is needed
+    return quantity{atan2(value_cast<rep>(y).numerical_value_in(unit), value_cast<rep>(x).numerical_value_in(unit)),
+                    radian};
+  } else
+    return quantity{atan2(y.numerical_value_in(unit), x.numerical_value_in(unit)), radian};
+}
+
 }  // namespace mp_units::si
diff --git a/test/runtime/math_test.cpp b/test/runtime/math_test.cpp
index c4e6b759..643f8378 100644
--- a/test/runtime/math_test.cpp
+++ b/test/runtime/math_test.cpp
@@ -375,6 +375,22 @@ TEST_CASE("SI inverse trigonometric functions", "[inv trig][si]")
   }
 }
 
+TEST_CASE("SI atan2 functions", "[atan2][si]")
+{
+  SECTION("atan2 should work on the same quantities")
+  {
+    REQUIRE_THAT(si::atan2(-1. * isq::length[km], 1. * isq::length[km]), AlmostEquals(-45. * deg));
+    REQUIRE_THAT(si::atan2(0. * isq::length[km], 1. * isq::length[km]), AlmostEquals(0. * deg));
+    REQUIRE_THAT(si::atan2(1. * isq::length[km], 1. * isq::length[km]), AlmostEquals(45. * deg));
+  }
+  SECTION("atan2 should work with different units of the same dimension")
+  {
+    REQUIRE_THAT(si::atan2(-1. * isq::length[km], 1000. * isq::length[m]), AlmostEquals(-45. * deg));
+    REQUIRE_THAT(si::atan2(0. * isq::length[km], 1000. * isq::length[m]), AlmostEquals(0. * deg));
+    REQUIRE_THAT(si::atan2(1. * isq::length[km], 1000. * isq::length[m]), AlmostEquals(45. * deg));
+  }
+}
+
 
 TEST_CASE("Angle trigonometric functions", "[trig][angle]")
 {
@@ -449,3 +465,23 @@ TEST_CASE("Angle inverse trigonometric functions", "[inv trig][angle]")
     REQUIRE_THAT(atan(1 * one), AlmostEquals(45. * angle[deg]));
   }
 }
+
+TEST_CASE("Angle atan2 functions", "[atan2][angle]")
+{
+  using namespace mp_units::angular;
+  using namespace mp_units::angular::unit_symbols;
+  using mp_units::angular::unit_symbols::deg;
+
+  SECTION("atan2 should work on the same quantities")
+  {
+    REQUIRE_THAT(atan2(-1. * isq::length[km], 1. * isq::length[km]), AlmostEquals(-45. * angle[deg]));
+    REQUIRE_THAT(atan2(0. * isq::length[km], 1. * isq::length[km]), AlmostEquals(0. * angle[deg]));
+    REQUIRE_THAT(atan2(1. * isq::length[km], 1. * isq::length[km]), AlmostEquals(45. * angle[deg]));
+  }
+  SECTION("atan2 should work with different units of the same dimension")
+  {
+    REQUIRE_THAT(atan2(-1. * isq::length[km], 1000. * isq::length[m]), AlmostEquals(-45. * angle[deg]));
+    REQUIRE_THAT(atan2(0. * isq::length[km], 1000. * isq::length[m]), AlmostEquals(0. * angle[deg]));
+    REQUIRE_THAT(atan2(1. * isq::length[km], 1000. * isq::length[m]), AlmostEquals(45. * angle[deg]));
+  }
+}