diff --git a/src/metrics/kerr_schild.h b/src/metrics/kerr_schild.h
index 14441f99..4df88cf4 100644
--- a/src/metrics/kerr_schild.h
+++ b/src/metrics/kerr_schild.h
@@ -40,6 +40,7 @@ namespace metric {
 
     const real_t dr, dtheta, dphi;
     const real_t dr_inv, dtheta_inv, dphi_inv;
+    const bool   small_angle;
 
     Inline auto Delta(real_t r) const -> real_t {
       return SQR(r) - TWO * r + SQR(a);
@@ -85,7 +86,8 @@ namespace metric {
       , dphi { (x3_max - x3_min) / nx3 }
       , dr_inv { ONE / dr }
       , dtheta_inv { ONE / dtheta }
-      , dphi_inv { ONE / dphi } {
+      , dphi_inv { ONE / dphi }
+      , small_angle { HALF * dtheta < constant::SMALL_ANGLE } {
       set_dxMin(find_dxMin());
     }
 
@@ -443,12 +445,21 @@ namespace metric {
     /**
      * differential area at the pole (used in axisymmetric solvers)
      * @param x1 radial coordinate along the axis (code units)
+     * @note uses small-angle approximation when the resolution is too high
      */
     Inline auto polar_area(real_t x1) const -> real_t {
-      return dr * (SQR(x1 * dr + x1_min) + SQR(a)) *
-             math::sqrt(ONE + TWO * (x1 * dr + x1_min) /
-                                (SQR(x1 * dr + x1_min) + SQR(a))) *
-             (ONE - math::cos(HALF * dtheta));
+      if (small_angle) {
+        return dr * (SQR(x1 * dr + x1_min) + SQR(a)) *
+              math::sqrt(ONE + TWO * (x1 * dr + x1_min) /
+                                  (SQR(x1 * dr + x1_min) + SQR(a))) *
+               (static_cast<real_t>(48) - SQR(dtheta)) * SQR(dtheta) /
+               static_cast<real_t>(384);
+      } else {
+        return dr * (SQR(x1 * dr + x1_min) + SQR(a)) *
+              math::sqrt(ONE + TWO * (x1 * dr + x1_min) /
+                                  (SQR(x1 * dr + x1_min) + SQR(a))) *
+              (ONE - math::cos(HALF * dtheta));
+      }
     }
 
     /**
diff --git a/src/metrics/qkerr_schild.h b/src/metrics/qkerr_schild.h
index 9b608f35..90f366ef 100644
--- a/src/metrics/qkerr_schild.h
+++ b/src/metrics/qkerr_schild.h
@@ -39,6 +39,7 @@ namespace metric {
     const real_t chi_min, eta_min, phi_min;
     const real_t dchi, deta, dphi;
     const real_t dchi_inv, deta_inv, dphi_inv;
+    const bool   small_angle;
 
     Inline auto Delta(real_t r) const -> real_t {
       return SQR(r) - TWO * r + SQR(a);
@@ -89,7 +90,8 @@ namespace metric {
       , dphi { (x3_max - phi_min) / nx3 }
       , dchi_inv { ONE / dchi }
       , deta_inv { ONE / deta }
-      , dphi_inv { ONE / dphi } {
+      , dphi_inv { ONE / dphi } 
+      , small_angle { eta2theta(HALF * deta) < constant::SMALL_ANGLE } {
       set_dxMin(find_dxMin());
     }
 
@@ -503,15 +505,27 @@ namespace metric {
      * differential area at the pole (used in axisymmetric solvers)
      * @note approximate solution for the polar area
      * @param x1 radial coordinate along the axis (code units)
+     * @note uses small-angle approximation when the resolution is too high
      */
     Inline auto polar_area(real_t x1) const -> real_t {
       if constexpr (D != Dim::_1D) {
-        return dchi * math::exp(x1 * dchi + chi_min) *
-               (SQR(r0 + math::exp(x1 * dchi + chi_min)) + SQR(a)) *
-               math::sqrt(
-                 ONE + TWO * (r0 + math::exp(x1 * dchi + chi_min)) /
-                         (SQR(r0 + math::exp(x1 * dchi + chi_min)) + SQR(a))) *
-               (ONE - math::cos(eta2theta(HALF * deta)));
+        if (small_angle) { 
+          const real_t dtheta = eta2theta(HALF * deta);
+          return dchi * math::exp(x1 * dchi + chi_min) *
+                (SQR(r0 + math::exp(x1 * dchi + chi_min)) + SQR(a)) *
+                math::sqrt(
+                  ONE + TWO * (r0 + math::exp(x1 * dchi + chi_min)) /
+                          (SQR(r0 + math::exp(x1 * dchi + chi_min)) + SQR(a))) *
+                (static_cast<real_t>(48) - SQR(dtheta)) * SQR(dtheta) /
+                 static_cast<real_t>(384);
+        } else {
+          return dchi * math::exp(x1 * dchi + chi_min) *
+                (SQR(r0 + math::exp(x1 * dchi + chi_min)) + SQR(a)) *
+                math::sqrt(
+                  ONE + TWO * (r0 + math::exp(x1 * dchi + chi_min)) /
+                          (SQR(r0 + math::exp(x1 * dchi + chi_min)) + SQR(a))) *
+                (ONE - math::cos(eta2theta(HALF * deta)));
+        }
       }
     }