igm_utils.py 4.8 KB
Newer Older
1
# -*- coding: utf-8 -*-
Yannick Roehlly's avatar
Yannick Roehlly committed
2
# Copyright (C) 2014 Centre de données Astrophysiques de Marseille
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 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
# Licensed under the CeCILL-v2 licence - see Licence_CeCILL_V2-en.txt
# Authors: Denis Burgarella, Yannick Roehlly

import math
import numpy as np


def ig_transmission_meiksin(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.

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

    """

    redshift = float(redshift)
    wavelength = np.array(wavelength, dtype=float)

    # For redshift below 2, we assume that there is no IGM attenuation.
    if redshift <= 2.:
        ig_tranmission = np.ones(len(wavelength))
        return ig_tranmission

    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.zeros(n_transitions_low)
    z_n = np.zeros((n_transitions_low, len(wavelength)))
    for n in range(2, n_transitions_low):
        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_transistions_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_low):
        if n <= 5:
            for l, index in enumerate(wavelength):
                if z_n[n, l] < 3:
                    tau_n[n, l] = (tau_a * fact[n] *
                                   np.power(0.25 * (1. + z_n[n, l]),
                                            (1. / 3.)))
                else:
                    tau_n[n, l] = (tau_a * fact[n] *
                                   np.power(0.25 * (1. + z_n[n, l]),
                                            (1. / 6.)))

        elif 5 < n <= 9:
            for l, index in enumerate(wavelength):
                tau_n[n, l] = (tau_a * fact[n] *
                               np.power(0.25 * (1. + z_n[n, l]),
                                        (1. / 3.)))

        else:
            for l, index in enumerate(wavelength):
                tau_n[n, l] = (tau_n[9, l] * 720. /
                               (float(n) * (float(n*n - 1.))))

    for n in range(2, n_transitions_low):
        for l, index in enumerate(wavelength):
            if z_n[n, l] >= redshift:
                tau_n[n, l] = 0.

    z_l = wavelength / lambda_limit - 1.
    tau_l_igm = np.zeros(len(wavelength))
    for l, index in enumerate(wavelength):
        if z_l[l] >= redshift:
            tau_l_igm[l] = 0.
        else:
            tau_l_igm[l] = (0.805 * np.power(1. + z_l[l], 3) *
                            (1. / (1. + z_l[l]) - 1. / (1. + redshift)))

    term1 = gamma - np.exp(-1.)
    term2 = 0.
    for n in range(n_transitions_max):
        term2 = term2 + np.power(-1., n) / (math.factorial(n) * float(2*n - 1))
    term3 = ((1.+redshift) * np.power(wavelength/lambda_limit, 1.5) -
             np.power(wavelength/lambda_limit, 2.5))
    term4 = np.zeros(len(wavelength))
    for n in range(1, n_transitions_max):
        term4 = (term4 + (2.*np.power(-1., n) / (math.factorial(n) *
                                                 float((6*n - 5)*(2*n - 1)))) *
                 (np.power(1.+redshift, 2.5-float(3 * n)) *
                  np.power(wavelength/lambda_limit, 3*n) -
                  np.power(wavelength/lambda_limit, 2.5)))

    tau_l_lls = np.zeros(len(wavelength))
    for l, index in enumerate(wavelength):
        if z_l[l] >= redshift:
            tau_l_lls[l] = 0.
        else:
            tau_l_lls[l] = n0 * ((term1 - term2) * term3[l] - term4[l])

    tau_taun = 0.
    for n in range(n_transitions_max):
        tau_taun = tau_taun + tau_n[n, :]

    weight = np.zeros(len(wavelength))
    for l, index in enumerate(wavelength):
        weight[l] = 0.5*(1.+math.erf(0.01*(wavelength[l]-(1+redshift)*60.)))

    tau = tau_taun + tau_l_igm + tau_l_lls
    ig_tranmission = np.exp(- tau) * weight

    return ig_tranmission