diff --git a/src/boundary_condition.cpp b/src/boundary_condition.cpp
index 5bbda4830..d166f1247 100644
--- a/src/boundary_condition.cpp
+++ b/src/boundary_condition.cpp
@@ -38,6 +38,9 @@ void VacuumBC::handle_particle(Particle& p, const Surface& surf) const
 void ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const
 {
   Direction u = surf.reflect(p.r(), p.u(), &p);
+
+  // normalize reflected u to ensure no floating point error leads to
+  // unnormalized directions
   u /= u.norm();
 
   // Handle the effects of the surface albedo on the particle's weight.
@@ -53,6 +56,9 @@ void ReflectiveBC::handle_particle(Particle& p, const Surface& surf) const
 void WhiteBC::handle_particle(Particle& p, const Surface& surf) const
 {
   Direction u = surf.diffuse_reflect(p.r(), p.u(), p.current_seed());
+
+  // normalize outgoing u to ensure no floating point error leads to
+  // unnormalized directions
   u /= u.norm();
 
   // Handle the effects of the surface albedo on the particle's weight.
@@ -241,6 +247,10 @@ void RotationalPeriodicBC::handle_particle(
   new_u[axis_1_idx_] = cos_theta * u[axis_1_idx_] - sin_theta * u[axis_2_idx_];
   new_u[axis_2_idx_] = sin_theta * u[axis_1_idx_] + cos_theta * u[axis_2_idx_];
 
+  // normalize new_u to ensure no floating point error leads to unnormalized
+  // directions
+  new_u /= new_u.norm();
+
   // Handle the effects of the surface albedo on the particle's weight.
   BoundaryCondition::handle_albedo(p, surf);