diff --git a/src/mesh.cpp b/src/mesh.cpp
index acd90033f32..789113ccc31 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -1893,18 +1893,19 @@ double CylindricalMesh::find_r_crossing(
   const double inv_denominator = 1.0 / denominator;
 
   const double p = (u.x * r.x + u.y * r.y) * inv_denominator;
-  double c = r.x * r.x + r.y * r.y - r0 * r0;
-  double D = p * p - c * inv_denominator;
+  double R = std::sqrt(r.x * r.x + r.y * r.y);
+  double D = p * p - (R - r0) * (R + r0) * inv_denominator;
 
   if (D < 0.0)
     return INFTY;
 
   D = std::sqrt(D);
 
-  // the solution -p - D is always smaller as -p + D : Check this one first
-  if (std::abs(c) <= RADIAL_MESH_TOL)
+  // Particle is already on the shell surface; avoid spurious crossing
+  if (std::abs(R - r0) <= RADIAL_MESH_TOL * (1.0 + std::abs(r0)))
     return INFTY;
 
+  // Check -p - D first because it is always smaller as -p + D
   if (-p - D > l)
     return -p - D;
   if (-p + D > l)
@@ -2178,15 +2179,16 @@ double SphericalMesh::find_r_crossing(
   if (r0 == 0.0)
     return INFTY;
   const double p = r.dot(u);
-  double c = r.dot(r) - r0 * r0;
-  double D = p * p - c;
+  double R = r.norm();
+  double D = p * p - (R - r0) * (R + r0);
 
-  if (std::abs(c) <= RADIAL_MESH_TOL)
+  // Particle is already on the shell surface; avoid spurious crossing
+  if (std::abs(R - r0) <= RADIAL_MESH_TOL * (1.0 + std::abs(r0)))
     return INFTY;
 
   if (D >= 0.0) {
     D = std::sqrt(D);
-    // the solution -p - D is always smaller as -p + D : Check this one first
+    // Check -p - D first because it is always smaller as -p + D
     if (-p - D > l)
       return -p - D;
     if (-p + D > l)