Source code for pysm.utils

# Licensed under a 3-clause BSD style license - see LICENSE.rst

# This sub-module is destined for common non-package specific utility
# functions.

import numpy as np
from numba import njit

from .. import units as u

def has_polarization(m):
    """Checks if a map or a group of map is polarized

    Works with a map of shape (IQU, npix) or
    (channels, IQU, npix)"""
    if isinstance(m, np.ndarray):  # array
        if m.ndim == 1:
            return False
            return m.shape[-2] == 3
    elif isinstance(m[0], np.ndarray):  # IQU tuple
        return len(m) == 3
    elif isinstance(m[0][0], np.ndarray):  # multiple IQU tuples
        return len(m[0]) == 3
        raise TypeError("Map format not understood")

[docs]def normalize_weights(freqs, weights): if freqs.isscalar or len(freqs) == 1: return np.array([1.0]) else: if weights is None: weights = np.ones(len(freqs), dtype=np.float) weights = (weights * u.uK_RJ).to_value((u.Jy /, equivalencies=u.cmb_equivalencies(freqs)) return weights / np.trapz(weights, freqs.value)
[docs]def bandpass_unit_conversion(freqs, weights, output_unit): """Unit conversion from uK_RJ to output unit given a bandpass Parameters ---------- freqs : astropy.units.Quantity Frequency array in a unit compatible with GHz """ freqs = check_freq_input(freqs) factors = (np.ones(len(freqs), dtype=np.float) * u.uK_RJ).to_value(output_unit, equivalencies=u.cmb_equivalencies(freqs)) if len(freqs) > 1: w = normalize_weights(freqs, weights) factor = np.trapz(factors * w, freqs.value) else: factor = factors[0] return factor * u.Unit(u.Unit(output_unit) / u.uK_RJ)
@njit def trapz_step_inplace(freqs, weights, i, m, output): """Execute a step of the trapezoidal rule and accumulate into output freqs : ndarray Frequency axis, generally in GHz, but doesn't matter as long as weights were normalized accordingly weights : ndarray Frequency bandpass response, normalized to unit integral (with trapz) i : integer Index of the current step in the arrays m : ndarray Emission evaluated at the current frequency to be accumulated output : ndarray Array where the integrated emission is accumulated. """ # case for a single frequency, compensate for the .5 factor if i == 0 and len(freqs) == 1: delta_freq = 2 # first step of the integration elif i == 0: delta_freq = freqs[1] - freqs[0] # last step elif i == (len(freqs) - 1): delta_freq = freqs[-1] - freqs[-2] # middle steps else: delta_freq = freqs[i + 1] - freqs[i - 1] output += 0.5 * m * weights[i] * delta_freq
[docs]def check_freq_input(freqs): """ Function to check that the input to `Model.get_emission` is a np.ndarray. This function will convet input integers or arrays to a single element numpy array. Parameters ---------- freqs: int, float, list, ndarray Returns ------- ndarray Frequencies in numpy array form. """ if isinstance(freqs, np.ndarray): freqs = freqs elif isinstance(freqs, list): freqs = np.array(freqs) else: try: freqs = np.array([freqs]) except: print( """Could not make freqs into an ndarray, check input.""" ) raise if isinstance(freqs, u.Quantity): if freqs.isscalar: return freqs[None] return freqs return freqs * u.GHz