docs: def ang2pix(nside, theta, phi):
      r"""Convert angle :math:`\theta` :math:`\phi` to pixel.
  
          This is translated from chealpix.c; but refer to Section 4.1 of
          http://adsabs.harvard.edu/abs/2005ApJ...622..759G
      """
      nside, theta, phi = numpy.lib.stride_tricks.broadcast_arrays(nside, theta, phi)
  
      def equatorial(nside, tt, z):
          t1 = nside * (0.5 + tt)
          t2 = nside * z * 0.75
          jp = (t1 - t2).astype('i8')
          jm = (t1 + t2).astype('i8')
          ir = nside + 1 + jp - jm # in {1, 2n + 1}
          kshift = 1 - (ir & 1) # kshift=1 if ir even, 0 odd
  
          ip = (jp + jm - nside + kshift + 1) // 2 # in {0, 4n - 1}
  
          ip = ip % (4 * nside)
          return nside * (nside - 1) * 2 + (ir - 1) * 4 * nside + ip
  
      def polecaps(nside, tt, z, s):
          tp = tt - numpy.floor(tt)
          za = numpy.abs(z)
          tmp = nside * s / ((1 + za) / 3) ** 0.5
          mp = za > 0.99
          tmp[mp] = nside[mp] * (3 *(1-za[mp])) ** 0.5
          jp = (tp * tmp).astype('i8')
          jm = ((1 - tp) * tmp).astype('i8')
          ir = jp + jm + 1
          ip = (tt * ir).astype('i8')
          ip = ip % (4 * ir)
  
          r1 = 2 * ir * (ir - 1)
          r2 = 2 * ir * (ir + 1)
  
          r = numpy.empty_like(r1)
  
          r[z > 0] = r1[z > 0] + ip[z > 0]
          r[z < 0] = 12 * nside[z < 0] * nside[z < 0] - r2[z < 0] + ip[z < 0]
          return r
  
      z = numpy.cos(theta)
      s = numpy.sin(theta)
  
      tt = (phi / (0.5 * numpy.pi) ) % 4 # in [0, 4]
  
      result = numpy.zeros(z.shape, dtype='i8')
      mask = (z < 2. / 3) & (z > -2. / 3)
  
      result[mask] = equatorial(nside[mask], tt[mask], z[mask])
      result[~mask] = polecaps(nside[~mask], tt[~mask], z[~mask], s[~mask])
      return result