Sky temperature

t_sky.py

def load_T_sky():
    """
    Read the Haslam sky temperature map into a list from which temperatures can
    be retrieved. The temperature sky map is given in the weird units of
    HealPix, and despite looking up info on this coordinate system, I don't
    have the foggiest idea of how to transform these to galactic coordinates. I
    have therefore directly copied the following code from psrpoppy in the
    assumption Sam Bates managed to figure it out.
    Returns:
        t_sky_list (list): List of sky temperatures in HealPix? coordinates?
    """

    model = os.path.join(os.path.dirname(__file__), '../data/models/tsky/')
    path = os.path.join(model, 'haslam_2014.dat')

    t_sky_list = []
    with open(path) as f:
        for line in f:
            str_idx = 0
            while str_idx < len(line):
                # each temperature occupies space of 5 chars
                temp_string = line[str_idx:str_idx+5]
                try:
                    t_sky_list.append(float(temp_string))
                except:
                    pass
                str_idx += 5

    return t_sky_list
  
      def calc_T_sky(self, src):
        """
        Calculate the sky temperature from the Haslam table, before scaling to
        the survey frequency. The temperature sky map is given in the weird
        units of HealPix and despite looking up info on this coordinate system,
        I don't have the foggiest idea of how to transform these to galactic
        coordinates. I have therefore directly copied the following code from
        psrpoppy in the assumption Sam Bates managed to figure it out.
        Args:
            src (class): Needed for coordinates
        Returns:
            src.T_sky (float): Sky temperature [K]
        """

        # ensure l is in range 0 -> 360
        B = src.gb
        if src.gl < 0.:
            L = 360 + src.gl
        else:
            L = src.gl

        # convert from l and b to list indices
        j = B + 90.5
        if j > 179:
            j = 179

        nl = L - 0.5
        if L < 0.5:
            nl = 359
        i = float(nl) / 4.

        T_sky_haslam = self.T_sky_list[180*int(i) + int(j)]

        # scale temperature
        # Assuming dominated by syncrotron radiation
        T_sky = T_sky_haslam * (self.central_freq/408.0)**(-2.6)

        src.T_sky = T_sky