brick_utils.py 1.33 KB
Newer Older
Wu You's avatar
Wu You committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
import numpy as np
import healpy as hp



def angular_distance(ra1, dec1, ra2, dec2):
    """
    Calculate the angular distance between two points on a sphere.

    Parameters
    ----------
    ra1 : float
        Right Ascension of the first point (radians).
    dec1 : float
        Declination of the first point (radians).
    ra2 : float
        Right Ascension of the second point (radians).
    dec2 : float
        Declination of the second point (radians).

    Returns
    -------
    float
        Angular distance between the two points (radians).
    """
    return np.arccos(
        np.sin(dec1) * np.sin(dec2) +
        np.cos(dec1) * np.cos(dec2) * np.cos(ra1 - ra2)
    )

def radius_exact(nside, pix_id):
    """
    Compute the exact minimum enclosing circle radius of a HEALPix pixel.

    Parameters
    ----------
    nside : int
        HEALPix resolution parameter.
    pix_id : int
        HEALPix pixel ID.

    Returns
    -------
    float
        Exact minimum enclosing circle radius of the pixel (degrees).
    """
    theta_c, phi_c = hp.pix2ang(nside, pix_id)
    verts = hp.boundaries(nside, pix_id, step=1)
    ra_v = np.arctan2(verts[1], verts[0])
    dec_v = np.arcsin(verts[2])
    rads = [angular_distance(phi_c, np.pi/2 - theta_c, rv, dv)
            for rv, dv in zip(ra_v, dec_v)]
    return np.degrees(max(rads))