This method selects a subset of particles that are located in a solid angle defined by a center and a set of points. The solid angle is constructed as a union of small cones whose axis goes through the center and a point. So there's such a cone for each position, and a point is in the solid angle if it lies within any of these cones. The width of the cones can be adjusted.

The method is implemented by partitioning the surface of the unit sphere into bins using the polar coordinates $(\theta, \phi)$ . The partitioning is always uniform in the zenith angle $\theta$ , while the partitioning in the azimuthal angle $\phi$ varies. For each reference point, the unit vector from the center to the point is constructed, and it is stored in all the bins that overlap with the cone defined by the point. Bins that are completely covered by a single cone are marked as such. Checking whether a point is in the solid angle is then straightforward with this data structure: one finds the bin that corresponds to the point, and checks whether the bin is completely covered. If it is not, one additionally needs to check whether it is within the specified cutoff of any of the stored points.

The above construction gives quite a lot of flexibility for constructing the bins without modifying the rest of the code. The current (quite inefficient) implementation is discussed below, but it should be optimized to get the most out of the code.

The current way of constructing the bins constructs the boundaries statically: the bin size in the zenith direction is set to approximately half the angle cutoff, and the bins in the azimuthal direction have sizes such that the shortest edge of the bin is approximately equal to half the angle cutoff (for the regions close to the poles, a single bin is used). Each reference point is then added to the bins as follows:

Find the zenith angle range that is spanned by the cone centered at the point (this is simple addition/subtraction).
Calculate the maximal span of the cone in the azimuthal direction using the formula
$\sin \Delta \phi_{max} = \frac{\sin \alpha}{\sin \theta}$
(a sine formula in spherical coordinates), where $\alpha$ is the width of the cone and $\theta$ is the zenith angle of the cone center. Similarly, the zenith angle at which this extent is achieved is calculated using
$\cos \theta_{max} = \frac{\cos \theta}{\cos \alpha}$
(Pythagoras's theorem in spherical coordinates).
For each zenith angle bin that is at least partially covered by the cone, calculate the span of the cone at the edges using
$\sin^2 \frac{\Delta \phi}{2} = \frac{\sin^2 \frac{\alpha}{2} - \sin^2 \frac{\theta - \theta'}{2}}{\sin \theta \sin \theta'}$
(distance in spherical geometry), where $\theta'$ is the zenith angle of the bin edge. Treat zenith angle bins that are completely covered by the cone (in the case that the cone is centered close to the pole) as a special case.
Using the values calculated above, loop through the azimuthal bins that are partially or completely covered by the cone and update them.

The total solid angle (for covered fraction calculations) is estimated by taking the total area of completely covered bins plus half the area of partially covered bins. The second one is an approximation, but should give reasonable estimates for the averages as well as in cases where the bin size is small.