<p>The aim of this work is to propose a new approach for computing the integral of a product of several real spherical harmonics over the full unit sphere as well as over arbitrary hemispheres. Such integrals arise in the computation of the entries of elementary matrices when solving the neutral particle transport equation using numerical schemes that combine the discontinuous Galerkin method for the spatial variable (with a polyhedral mesh) and the spherical harmonics method for handling the angular variable. While computing these integrals over the entire unit sphere remains relatively straightforward, their evaluation over spherical subdomains, such as hemispheres, presents additional challenges due to the complexity of the integration bounds. To overcome these difficulties, we propose a method based on an appropriate change of variables via rotation, taking advantage of the invariance of real spherical harmonics under such rotations.</p>