redshifting.py

# -*- coding: utf-8 -*-
# Copyright (C) 2014 Yannick Roehlly, Médéric Boquien, Denis Burgarella
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
# Author: Yannick Roehlly, Médéric Boquien, Denis Burgarella, Guang Yang

"""
Redshifting module
==================

This module implements the redshifting of a SED. The SED must be rest-frame
or the module will raise en exception when processing it.

Note that this module, contrary to the other SED creation modules, actually
changes the individual luminosity contributions as it redshifts everyone.
Also note that doing this, this module does not rely on the SED object
interface but on its inner implementations. That means that if the SED object
is changed, this module may need to be adapted.

"""

from collections import OrderedDict

import numpy as np
from scipy.constants import parsec
from scipy.special import factorial
from ..utils.cosmology import age, luminosity_distance

from . import SedModule

def igm_transmission(wavelength, redshift):
    """Intergalactic transmission (Meiksin, 2006)

    Compute the intergalactic transmission as described in Meiksin, 2006.

    Parameters
    ----------
    wavelength: array like of floats
        The wavelength(s) in nm.
    redshift: float
        The redshift. Must be strictly positive.

    Returns
    -------
    igm_transmission: numpy array of floats
        The intergalactic transmission at each input wavelength.

    """
    n_transitions_low = 10
    n_transitions_max = 31
    gamma = 0.2788  # Gamma(0.5,1) i.e., Gamma(2-beta,1) with beta = 1.5
    n0 = 0.25
    lambda_limit = 91.2  # Lyman limit in nm

    lambda_n = np.empty(n_transitions_max)
    z_n = np.empty((n_transitions_max, len(wavelength)))
    for n in range(2, n_transitions_max):
        lambda_n[n] = lambda_limit / (1. - 1. / float(n*n))
        z_n[n, :] = (wavelength / lambda_n[n]) - 1.

    # From Table 1 in Meiksin (2006), only n >= 3 are relevant.
    # fact has a length equal to n_transitions_low.
    fact = np.array([1., 1., 1., 0.348, 0.179, 0.109, 0.0722, 0.0508, 0.0373,
                     0.0283])

    # First, tau_alpha is the mean Lyman alpha transmitted flux,
    # Here n = 2 => tau_2 = tau_alpha
    tau_n = np.zeros((n_transitions_max, len(wavelength)))
    if redshift <= 4:
        tau_a = 0.00211 * np.power(1. + redshift,  3.7)
        tau_n[2, :] = 0.00211 * np.power(1. + z_n[2, :], 3.7)
    elif redshift > 4:
        tau_a = 0.00058 * np.power(1. + redshift,  4.5)
        tau_n[2, :] = 0.00058 * np.power(1. + z_n[2, :], 4.5)

    # Then, tau_n is the mean optical depth value for transitions
    # n = 3 - 9 -> 1
    for n in range(3, n_transitions_max):
        if n <= 5:
            w = z_n[n, :] < 3
            tau_n[n, w] = (tau_a * fact[n] *
                           np.power(0.25 * (1. + z_n[n, w]), (1. / 3.)))
            w = z_n[n, :] >= 3
            tau_n[n, w] = (tau_a * fact[n] *
                           np.power(0.25 * (1. + z_n[n, w]), (1. / 6.)))
        elif 5 < n <= 9:
            tau_n[n, :] = (tau_a * fact[n] *
                           np.power(0.25 * (1. + z_n[n, :]), (1. / 3.)))
        else:
            tau_n[n, :] = (tau_n[9, :] * 720. /
                           (float(n) * (float(n*n - 1.))))

    for n in range(2, n_transitions_max):
        # If z_n>=redshift or z_n<0, the photon cannot be absorbed by Lyman n->1
        w = (z_n[n, :] >= redshift) | (z_n[n, :] < 0)
        tau_n[n, w] = 0.

    z_l = wavelength / lambda_limit - 1.
    w = z_l < redshift

    tau_l_igm = np.zeros_like(wavelength)
    tau_l_igm[w] = (0.805 * np.power(1. + z_l[w], 3) *
                    (1. / (1. + z_l[w]) - 1. / (1. + redshift)))

    term1 = gamma - np.exp(-1.)

    n = np.arange(n_transitions_low - 1)
    term2 = np.sum(np.power(-1., n) / (factorial(n) * (2*n - 1)))

    term3 = ((1.+redshift) * np.power(wavelength[w]/lambda_limit, 1.5) -
             np.power(wavelength[w]/lambda_limit, 2.5))

    term4 = np.sum(np.array(
        [((2.*np.power(-1., n) / (factorial(n) * ((6*n - 5)*(2*n - 1)))) *
          ((1.+redshift) ** (2.5-(3 * n)) *
           (wavelength[w]/lambda_limit) ** (3*n) -
           (wavelength[w]/lambda_limit) ** 2.5))
         for n in np.arange(1, n_transitions_low)]), axis=0)

    tau_l_lls = np.zeros_like(wavelength)
    tau_l_lls[w] = n0 * ((term1 - term2) * term3 - term4)

    # Reset for short wavelength (z_l<0)
    w = z_l<0
    # Get the normalization factor at z_l=0
    tau_norm_l_igm = np.interp(0, z_l, tau_l_igm)
    tau_norm_l_lls = np.interp(0, z_l, tau_l_lls)
    # Calculate tau_l_igm & tau_l_lls, assuming cross section ~ lambda^2.75 (O'Meara et al. 2013)
    tau_l_igm[w] = tau_norm_l_igm * (z_l[w]+1)**2.75
    tau_l_lls[w] = tau_norm_l_lls * (z_l[w]+1)**2.75

    tau_taun = np.sum(tau_n[2:n_transitions_max, :], axis=0)

    igm_transmission = np.exp(-tau_taun-tau_l_igm-tau_l_lls)

    return igm_transmission


class Redshifting(SedModule):
    """Redshift a SED

    This module redshift a rest-frame SED. If the SED is already redshifted, an
    exception is raised.

    """

    parameter_list = OrderedDict([
        ("redshift", (
            "cigale_list(minvalue=0.)",
            "Redshift of the objects. Leave empty to use the redshifts from the"
            " input file.",
            None
        ))
    ])

    def _init_code(self):
        """Compute the age of the Universe at a given redshift
        """
        self.redshift = float(self.parameters["redshift"])

        # Raise an error when applying a negative redshift. This module is
        # not for blue-shifting.
        if self.redshift < 0.:
            raise Exception(f"The redshift provided is negative "
                            f"({self.redshift}).")

        self.universe_age = age(self.redshift)
        self.luminosity_distance = luminosity_distance(self.redshift)
        # We do not define the values of the IGM attenuation component yet.
        # This is because we need the wavelength grid for that first. This
        # will be assigned on the first call.
        self.igm_attenuation = {}

    def process(self, sed):
        """Redshift the SED

        Parameters
        ----------
        sed: pcigale.sed.SED object

        """
        redshift = self.redshift

        # If the SED is already redshifted, raise an error.
        if ('universe.redshift' in sed.info and
            sed.info['universe.redshift'] > 0.):
            raise Exception(f"The SED is already redshifted (z="
                            f"{sed.info['universe.redshift']}).")

        if redshift > 0.:
            # We redshift directly the SED wavelength grid
            sed.wavelength_grid *= 1. + redshift

            # We modify each luminosity contribution to keep energy constant
            sed.luminosities *= 1. / (1. + redshift)
            sed.luminosity *= 1. / (1. + redshift)

        sed.add_info("universe.redshift", redshift)
        sed.add_info("universe.luminosity_distance", self.luminosity_distance)
        sed.add_info("universe.age", self.universe_age)

        # We identify the right grid from the length of the wavelength array.
        # It is not completely foolproof but it is good enough. We need to do
        # that in case two models do not have the same wavelength sampling.
        # This is the case for instance if some but not all models have an AGN
        # fraction of 0.
        key = sed.wavelength_grid.size
        if key not in self.igm_attenuation:
            self.igm_attenuation[key] = igm_transmission(sed.wavelength_grid,
                                                         redshift) - 1.
        sed.add_contribution('igm', sed.wavelength_grid,
                             self.igm_attenuation[key] * sed.luminosity)
        sed.add_module(self.name, self.parameters)


# SedModule to be returned by get_module
Module = Redshifting