lb_to_radec.py

def lb_to_radec(l, b):
    """
    Convert galactic coordinates to RA, Dec

    Formulas from 'An Introduction to Modern Astrophysics (2nd Edition)' by
    Bradley W. Carroll, Dale A. Ostlie.

    Args:
        l (float): Galactic longitude [fractional degrees]
        b (float): Galactic latitude [fractional degrees]

    Returns:
        ra, dec (float): Right ascension and declination [fractional degrees]
    """

    gl = math.radians(l)
    gb = math.radians(b)

    # Coordinates of the galactic north pole (J2000)
    a_ngp = math.radians(12.9406333 * 15. + 180.)
    d_ngp = math.radians(27.1282500)
    l_ngp = math.radians(123.9320000)

    sd_ngp = math.sin(d_ngp)
    cd_ngp = math.cos(d_ngp)
    sb = math.sin(gb)
    cb = math.cos(gb)

    # Calculate right ascension
    y = cb*math.sin(l_ngp - gl)
    x = cd_ngp*sb - sd_ngp*cb*math.cos(l_ngp - gl)
    ra = math.atan(y/x) + a_ngp
    ra = math.degrees(ra)

    # Ensure value is in right quadrant
    if y > 0 and x > 0:
        ra = 180 + ra % 90
    elif y > 0 and x < 0:
        ra = 270 + ra % 90
    elif y < 0 and x < 0:
        ra = ra % 90
    elif y < 0 and x > 0:
        ra = 90 + (ra % 90)

    # Calculate declination
    dec = math.asin(sd_ngp*sb + cd_ngp*cb*math.cos(l_ngp - gl))
    dec = math.degrees(dec) % 360.
    if dec > 270:
        dec = - (360 - dec)
    return ra, dec

def lb_to_radec(l, b):
    """
    Convert galactic coordinates to RA, Dec using astropy. More accurate than above,
    but slower.

    Args:
        l (float): Galactic longitude [fractional degrees]
        b (float): Galactic latitude [fractional degrees]

    Returns:
        ra, dec (float): Right ascension and declination [fractional degrees]
    """

    from astropy.coordinates import SkyCoord
    from astropy import units as u
    
    c = SkyCoord(l=gl*u.degree, b=gb*u.degree, frame='galactic')
    ra = c.fk5.ra.degree
    dec = c.fk5.dec.degree
        
    return ra, dec