Fix a small error in the construction of the circular sky grid (#5175)

Author: titodalcanton

bin/pycbc_make_sky_grid

@@ -8,7 +8,7 @@ gravitational wave triggers and calculate the coherent SNRs and related
 statistics. 
 
 The grid is constructed following the method described in Section V of
-https://arxiv.org/abs/1410.6042.
+Williamson et al 2014, https://arxiv.org/abs/1410.6042.
 
 Please refer to help(pycbc.types.angle_as_radians) for the recommended
 configuration file syntax for angle arguments.
@@ -84,24 +84,31 @@ def make_multi_det_grid(args):
     ]
     angular_spacing = min(ang_spacings)
 
-    sky_points = np.zeros((1, 2))
+    # The first point is always exactly at the North pole
+    sky_points = [
+        [0, np.pi / 2]
+    ]
 
+    # Generate the sky grid centered at the North pole
     number_of_rings = int(args.sky_error / angular_spacing)
+    for i in range(1, number_of_rings + 1):
+        # Note: Williamson et al 2014 states that the number of points in the
+        # i-th ring is 2πi / x, with x the angular spacing. This is not correct.
+        # In the approximation where the error radius is tiny (i.e. the initial
+        # grid stays close to the North pole) the correct number is 2πi.
+        # However, here we implement the general case, correct over the whole
+        # extent of the polar angle, for which the number is instead
+        # 2π cos(π/2 - ix) / x.
+        polar_angle = np.pi / 2 - i * angular_spacing
+        number_of_points = round(
+            2 * np.pi * np.cos(polar_angle) / angular_spacing
+        )
+        sky_points += [
+            [2 * np.pi * j / number_of_points, polar_angle]
+            for j in range(number_of_points)
+        ]
 
-    # Generate the sky grid centered at the North pole
-    for i in range(number_of_rings + 1):
-        if i == 0:
-            sky_points[0][0] = 0
-            sky_points[0][1] = np.pi / 2
-        else:
-            number_of_points = int(2 * np.pi * i)
-            for j in range(number_of_points):
-                sky_points = np.row_stack(
-                    (
-                        sky_points,
-                        np.array([j / i, np.pi / 2 - i * angular_spacing]),
-                    )
-                )
+    sky_points = np.array(sky_points)
 
     # Convert spherical coordinates to cartesian coordinates
     cart = spher_to_cart(sky_points)