diff --git a/src/distribution_multi.cpp b/src/distribution_multi.cpp
index ea24e3ea34..61d1570245 100644
--- a/src/distribution_multi.cpp
+++ b/src/distribution_multi.cpp
@@ -62,11 +62,6 @@ Direction PolarAzimuthal::sample(uint64_t* seed) const
   // Sample azimuthal angle
   double phi = phi_->sample(seed);
 
-  // If the reference direction is along the z-axis, rotate the aziumthal angle
-  // to match spherical coordinate conventions.
-  // TODO: apply this change directly to rotate_angle
-  if (u_ref_.x == 0 && u_ref_.y == 0) phi += 0.5*PI;
-
   return rotate_angle(u_ref_, mu, &phi, seed);
 }
 
diff --git a/src/math_functions.cpp b/src/math_functions.cpp
index ac60c6ce0a..449c50467b 100644
--- a/src/math_functions.cpp
+++ b/src/math_functions.cpp
@@ -662,15 +662,9 @@ Direction rotate_angle(Direction u, double mu, const double* phi, uint64_t* seed
             mu*u.z - a*b*cosphi};
   } else {
     b = std::sqrt(1. - u.y*u.y);
-    return {mu*u.x + a*(u.x*u.y*cosphi + u.z*sinphi) / b,
-            mu*u.y - a*b*cosphi,
-            mu*u.z + a*(u.y*u.z*cosphi - u.x*sinphi) / b};
-    // TODO: use the following code to make PolarAzimuthal distributions match
-    // spherical coordinate conventions. Remove the related fixup code in
-    // PolarAzimuthal::sample.
-    //return {mu*u.x + a*(-u.x*u.y*sinphi + u.z*cosphi) / b,
-    //        mu*u.y + a*b*sinphi,
-    //        mu*u.z - a*(u.y*u.z*sinphi + u.x*cosphi) / b};
+    return {mu*u.x + a*(-u.x*u.y*sinphi + u.z*cosphi) / b,
+            mu*u.y + a*b*sinphi,
+            mu*u.z - a*(u.y*u.z*sinphi + u.x*cosphi) / b};
   }
 }