diff --git a/src/particle.cpp b/src/particle.cpp
index b6a3f5dee1..6e87691957 100644
--- a/src/particle.cpp
+++ b/src/particle.cpp
@@ -428,7 +428,8 @@ Particle::cross_surface()
     }
     return;
 
-  } else if (surf->bc_ == BC_REFLECT && (settings::run_mode != RUN_MODE_PLOTTING)) {
+  } else if ((surf->bc_ == BC_REFLECT || surf->bc_ == BC_WHITE) 
+                                    && (settings::run_mode != RUN_MODE_PLOTTING)) {
     // =======================================================================
     // PARTICLE REFLECTS FROM SURFACE
 
@@ -457,10 +458,15 @@ Particle::cross_surface()
       score_surface_tally(this, model::active_meshsurf_tallies);
       this->r() = r;
     }
-
-    // Reflect particle off surface
-    Direction u = surf->reflect(this->r(), this->u());
-
+   
+    Direction u;
+    if(surf->bc_ == BC_REFLECT) {
+      // Reflect particle off surface
+      u = surf->reflect(this->r(), this->u());
+    } else if(surf->bc_ == BC_WHITE) {
+      // Diffuse reflect particle off surface
+      u = surf->diffuse_reflect(this->r(), this->u());
+    }
     // Make sure new particle direction is normalized
     this->u() = u / u.norm();
 
@@ -490,56 +496,6 @@ Particle::cross_surface()
     }
     return;
 
-  } else if (surf->bc_ == BC_WHITE && settings::run_mode != RUN_MODE_PLOTTING) {
-    // =======================================================================
-    // WHITE BOUNDARY
-
-    // Do not handle white boundary conditions on lower universes
-    if (n_coord_ != 1) {
-      this->mark_as_lost("Cannot diffuse reflect particle " + std::to_string(id_) +
-        " off surface in a lower universe.");
-      return;
-    }
-
-    if (!model::active_surface_tallies.empty()) {
-      score_surface_tally(this, model::active_surface_tallies);
-    }
-
-    if (!model::active_meshsurf_tallies.empty()) {
-      Position r {this->r()};
-      this->r() -= TINY_BIT * this->u();
-      score_surface_tally(this, model::active_meshsurf_tallies);
-      this->r() = r;
-    }
-
-    // Diffuse Reflect particle off surface
-    Direction u = surf->diffuse_reflect(this->r(), this->u()); 
-
-    // Make sure new particle direction is normalized
-    this->u() = u / u.norm();
-
-    // Reassign particle's cell and surface
-    coord_[0].cell = cell_last_[n_coord_last_ - 1];
-    surface_ = -surface_;
-
-    if (!settings::dagmc) {
-      n_coord_ = 1;
-      if (!find_cell(this, true)) {
-        this->mark_as_lost("Couldn't find particle after diffuse reflecting from surface "
-                           + std::to_string(surf->id_) + ".");
-        return;
-      }
-    }
-
-    // Set previous coordinate going slightly past surface crossing
-    // r_last_current_ = this->r() + TINY_BIT*this->u();
-
-    // Diagnostic message
-    if (settings::verbosity >= 10 || simulation::trace) {
-      write_message("    Diffuse reflected from surface " + std::to_string(surf->id_));
-    }
-    return;
-
   } else if (surf->bc_ == BC_PERIODIC && settings::run_mode != RUN_MODE_PLOTTING) {
     // =======================================================================
     // PERIODIC BOUNDARY
diff --git a/src/surface.cpp b/src/surface.cpp
index f759e7752b..2354659520 100644
--- a/src/surface.cpp
+++ b/src/surface.cpp
@@ -12,6 +12,7 @@
 #include "openmc/string_utils.h"
 #include "openmc/xml_interface.h"
 #include "openmc/random_lcg.h"
+#include "openmc/math_functions.h"
 
 namespace openmc {
 
@@ -148,7 +149,7 @@ Surface::Surface(pugi::xml_node surf_node)
     } else if (surf_bc == "reflective" || surf_bc == "reflect"
                || surf_bc == "reflecting") {
       bc_ = BC_REFLECT;
-    } else if (surf_bc == "white" || surf_bc == "diffuse") {
+    } else if (surf_bc == "white") {
       bc_ = BC_WHITE;
     } else if (surf_bc == "periodic") {
       bc_ = BC_PERIODIC;
@@ -201,51 +202,19 @@ Surface::diffuse_reflect(Position r, Direction u) const
   // Diffuse reflect direction according to the normal.
   // cosine distribution 
   
-  Direction n = normal(r);
-  n = n/n.norm();
-  const double projection = std::max(-1.0, std::min(1.0, n.dot(u)));
+  Direction n = this->normal(r);
+  n /= n.norm();
+  const double projection = n.dot(u);
   
   // sample from inverse function, u=sqrt(rand) since p(u)=2u, so F(u)=u^2 
-  const double cosine = (projection>=0.0) ? 
+  const double mu = (projection>=0.0) ? 
                   -std::sqrt(prn()) : std::sqrt(prn());  
   
-  // handle special case for very large cosine  
-  if(std::fabs(cosine)>1.0001){
-    std::cout<<"cosine="<<"\n";
-    exit(-1);
-  }else if(std::fabs(cosine)>=1.000){
-    return u = cosine * n;
-  }
-
   // sample azimuthal distribution uniformly 
-  double t1, t2, rr, ss, tt;
-  for(;;){
-    t1 = 2.0 * prn() - 1.0;
-    t2 = 2.0 * prn() - 1.0;  
-    rr = t1 * t1 + t2 * t2;
-    if(rr<=1.0) break;
-  }
-  rr = std::sqrt((1.0 - cosine * cosine)/rr);
-  t1 = t1 * rr;
-  t2 = t2 * rr;
-  
-  // determine the direction off surface relative to norm 
-  if (std::fabs(n[2]<=0.9)) { // for the general normal
-    tt = (std::sqrt(n[0]*n[0]+n[1]*n[1]));
-    ss = 1.0/tt;
-    u[0] = n[0] * cosine + (t1*n[0]*n[2]-t2*n[1])*ss;
-    u[1] = n[1] * cosine + (t1*n[1]*n[2]+t2*n[0])*ss;
-    u[2] = n[2] * cosine - t1*tt;    
-  } else {  // for the normal is almost along with z-axis  
-    tt = (std::sqrt(n[0]*n[0]+n[2]*n[2]));
-    ss = 1.0/tt;   
-    u[0] = n[0] * cosine + (t1*n[0]*n[1]+t2*n[2])*ss;
-    u[1] = n[1] * cosine - t1*tt;
-    u[2] = n[2] * cosine + (t1*n[1]*n[2]-t2*n[0])*ss; 
-  }
+  u = rotate_angle(n, mu, nullptr);
   
   // normalize the direction 
-  return u = u/u.norm();    
+  return u/u.norm();    
 }
 
 CSGSurface::CSGSurface() : Surface{} {};