healpix.py 文件源码

def ang2pix_ring_single(nside, theta, phi):
    """Calculate the pixel indexes in RING ordering scheme for one single
    pair of angular coordinate on the sphere.

    Parameters
    ----------
    theta : float
        The polar angle (i.e., latitude), ? ? [0, ?]. (unit: rad)
    phi : float
        The azimuthal angle (i.e., longitude), ? ? [0, 2?). (unit: rad)

    Returns
    -------
    ipix : int
        The index of the pixel corresponding to the input coordinate.

    NOTE
    ----
    * Only support the *RING* ordering scheme
    * This is the JIT-optimized version that partially replaces the
      ``healpy.ang2pix``

    References
    ----------
    - HEALPix software: ``src/C/subs/chealpix.c``: ``ang2pix_ring_z_phi()``
      http://healpix.sourceforge.net/
    """
    z = np.cos(theta)  # colatitude
    za = np.absolute(z)
    tt = (2.0 / np.pi) * np.remainder(phi, 2*np.pi)  # range: [0, 4)
    if za <= 2.0/3.0:
        # Equatorial region
        temp1 = nside * (tt + 0.5)
        temp2 = nside * z * 0.75
        jp = int(temp1 - temp2)  # Index of ascending edge line
        jm = int(temp1 + temp2)  # Index of descending edge line
        # Ring number counted from z=2/3
        iring = nside + 1 + jp - jm  # range: [1, 2n+1]
        kshift = 1 - (iring & 1)  # kshift=1 if ir even, 0 otherwise
        ip = int((jp + jm - nside + kshift + 1) / 2)
        ip = np.remainder(ip, 4*nside)
        ipix = nside * (nside-1) * 2 + (iring-1) * 4 * nside + ip
    else:
        # North & South polar caps
        tp = tt - int(tt)
        tmp = nside * np.sqrt(3 * (1-za))
        jp = int(tp * tmp)
        jm = int((1.0-tp) * tmp)
        # Ring number counted from the closest pole
        iring = jp + jm + 1
        ip = int(tt * iring)
        ip = np.remainder(ip, 4*iring)
        #
        if z > 0:
            ipix = 2 * iring * (iring-1) + ip
        else:
            ipix = 12 * nside * nside - 2 * iring * (iring+1) + ip
    #
    return ipix