Add TODOs to rotate_angle for a change in azimuth

src/distribution_multi.cpp

@@ -59,8 +59,14 @@ Direction PolarAzimuthal::sample() const
   double mu = mu_->sample();
   if (mu == 1.0) return u_ref_;
 
-  // Sample the azimuthal angle with an offset to match azimuth conventions
-  double phi = phi_->sample() + 0.5*PI;
+  // Sample azimuthal angle
+  double phi = phi_->sample();
+
+  // 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);
 }
 

src/math_functions.cpp

@@ -665,6 +665,12 @@ Direction rotate_angle(Direction u, double mu, const double* phi)
     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};
   }
 }