From 243c2495333c4bc15aaa345a2ba33330a66f2caf Mon Sep 17 00:00:00 2001
From: Lewis Gross <43077972+lewisgross1296@users.noreply.github.com>
Date: Wed, 15 Jul 2026 08:57:53 -0500
Subject: [PATCH] normalize new_u after periodic crossing to prevent
 unnormalized direction (#4015)

---
 src/boundary_condition.cpp | 10 ++++++++++
 1 file changed, 10 insertions(+)

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);