Function rotate_sc 

pub fn rotate_sc(
    u: [f64; 3],
    polar: (f64, f64),
    azimuth: (f64, f64),
) -> [f64; 3]

rotate with the polar and azimuthal angles given as (sin, cos) pairs, so a caller that already has them (the collision code gets the lab deflection from the half-angle tangent, without an inverse trigonometric call) pays for no sin_cos.

u must be a unit vector to rounding and polar must have sin^2 + cos^2 = 1 to rounding. The result is renormalised to first order, v (3 - |v|^2) / 2, which is exact to second order in |v|^2 - 1 (about 1e-16 here) and so costs a multiply instead of a square root and a division; the rounding drift of a long chain of rotations therefore does not accumulate.