diff --git a/src/global/utils/numeric.h b/src/global/utils/numeric.h
index afd66233..a3385336 100644
--- a/src/global/utils/numeric.h
+++ b/src/global/utils/numeric.h
@@ -80,18 +80,19 @@ inline constexpr double INV_64 = 0.015625;
 #define CROSS_x3(ax1, ax2, ax3, bx1, bx2, bx3) ((ax1) * (bx2) - (ax2) * (bx1))
 
 namespace constant {
-  inline constexpr std::uint64_t RandomSeed  = 0x123456789abcdef0;
-  inline constexpr double        HALF_PI     = 1.57079632679489661923;
-  inline constexpr double        PI          = 3.14159265358979323846;
-  inline constexpr double        INV_PI      = 0.31830988618379067154;
-  inline constexpr double        PI_SQR      = 9.86960440108935861882;
-  inline constexpr double        INV_PI_SQR  = 0.10132118364233777144;
-  inline constexpr double        TWO_PI      = 6.28318530717958647692;
-  inline constexpr double        E           = 2.71828182845904523536;
-  inline constexpr double        SQRT2       = 1.41421356237309504880;
-  inline constexpr double        INV_SQRT2   = 0.70710678118654752440;
-  inline constexpr double        SQRT3       = 1.73205080756887729352;
-  inline constexpr double        SMALL_ANGLE = 1e-3;
+  inline constexpr std::uint64_t RandomSeed       = 0x123456789abcdef0;
+  inline constexpr double        HALF_PI          = 1.57079632679489661923;
+  inline constexpr double        PI               = 3.14159265358979323846;
+  inline constexpr double        INV_PI           = 0.31830988618379067154;
+  inline constexpr double        PI_SQR           = 9.86960440108935861882;
+  inline constexpr double        INV_PI_SQR       = 0.10132118364233777144;
+  inline constexpr double        TWO_PI           = 6.28318530717958647692;
+  inline constexpr double        E                = 2.71828182845904523536;
+  inline constexpr double        SQRT2            = 1.41421356237309504880;
+  inline constexpr double        INV_SQRT2        = 0.70710678118654752440;
+  inline constexpr double        SQRT3            = 1.73205080756887729352;
+  inline constexpr double        SMALL_ANGLE      = 1e-3;
+  inline constexpr double        SMALL_ANGLE_GR   = 1e-5;
 } // namespace constant
 
 namespace convert {
diff --git a/src/kernels/currents_deposit.hpp b/src/kernels/currents_deposit.hpp
index 44bc52fd..06eeadad 100644
--- a/src/kernels/currents_deposit.hpp
+++ b/src/kernels/currents_deposit.hpp
@@ -127,10 +127,22 @@ namespace kernel {
                                                           vp);
           inv_energy = ONE / math::sqrt(ONE + NORM_SQR(ux1(p), ux2(p), ux3(p)));
         } else {
-          metric.template transform<Idx::D, Idx::U>(xp,
+          coord_t<Dim::_2D> xp_ { ZERO };
+          xp_[0] = xp[0];
+          real_t theta_Cd { xp[1] };
+          const real_t theta_Ph { metric.template convert<2, Crd::Cd, Crd::Ph>(theta_Cd) };
+          const real_t small_angle { constant::SMALL_ANGLE_GR };
+          const auto large_angle { constant::PI - small_angle };
+          if (theta_Ph < small_angle) {
+            theta_Cd = metric.template convert<2, Crd::Ph, Crd::Cd>(small_angle);
+          } else if (theta_Ph >= large_angle) {
+            theta_Cd = metric.template convert<2, Crd::Ph, Crd::Cd>(large_angle);
+          }
+          xp_[1] = theta_Cd;
+          metric.template transform<Idx::D, Idx::U>(xp_,
                                                     { ux1(p), ux2(p), ux3(p) },
                                                     vp);
-          inv_energy = metric.alpha(xp) / math::sqrt(ONE + ux1(p) * vp[0] + ux2(p) * vp[1] +
+          inv_energy = metric.alpha(xp_) / math::sqrt(ONE + ux1(p) * vp[0] + ux2(p) * vp[1] +
                                         ux3(p) * vp[2]);
         }
         if (Kokkos::isnan(vp[2]) || Kokkos::isinf(vp[2])) {
diff --git a/src/kernels/particle_pusher_gr.hpp b/src/kernels/particle_pusher_gr.hpp
index 2318d543..1eb12494 100644
--- a/src/kernels/particle_pusher_gr.hpp
+++ b/src/kernels/particle_pusher_gr.hpp
@@ -678,6 +678,19 @@ namespace kernel::gr {
       xp[0] = i_di_to_Xi(i1(p), dx1(p));
       xp[1] = i_di_to_Xi(i2(p), dx2(p));
 
+      coord_t<Dim::_2D> xp_ { ZERO };
+      xp_[0] = xp[0];
+      real_t theta_Cd { xp[1] };
+      const real_t theta_Ph { metric.template convert<2, Crd::Cd, Crd::Ph>(theta_Cd) };
+      const real_t small_angle { constant::SMALL_ANGLE_GR };
+      const auto large_angle { constant::PI - small_angle };
+      if (theta_Ph < small_angle) {
+        theta_Cd = metric.template convert<2, Crd::Ph, Crd::Cd>(small_angle);
+      } else if (theta_Ph >= large_angle) {
+        theta_Cd = metric.template convert<2, Crd::Ph, Crd::Cd>(large_angle);
+      }
+      xp_[1] = theta_Cd;
+
       vec_t<Dim::_3D> Dp_cntrv { ZERO }, Bp_cntrv { ZERO }, Dp_hat { ZERO },
         Bp_hat { ZERO };
       interpolateFields(p, Dp_cntrv, Bp_cntrv);
@@ -694,7 +707,7 @@ namespace kernel::gr {
       vp[0] = vp_upd[0];
       vp[1] = vp_upd[1];
       vp[2] = vp_upd[2];
-      GeodesicMomentumPush<Massive_t>(Massive_t {}, xp, vp, vp_upd);
+      GeodesicMomentumPush<Massive_t>(Massive_t {}, xp_, vp, vp_upd);
       /* u**_i(n) -> u_i(n + 1/2) */
       vp[0] = vp_upd[0];
       vp[1] = vp_upd[1];