Counts of sources from MOJAVE beamed LF (See Cara&Lister 2008)

mojave_vlba_counts.py

import numpy as np
from scipy.integrate import quad
from astropy.cosmology import LambdaCDM

# cm in pc
pc = 3.086 * 1E+18


def beta(Gamma):
    """
    Velocity in units of speed of light [c].

    :param Gamma:
        Lorentz factor.
    """
    return np.sqrt(Gamma**2.-1.)/Gamma


def beamed_LF(alpha=-2.53, m=1.4, z0=1.29, sigma=0.76, n0=2.87*10**(-10),
              l1=10**22.2, l2=10**29.1, p=2.0, gamma1=1.25, gamma2=32, k=-1.5,
              L_star=10.0**27, L1=10.0**25, L2=1.1*10**29):
    """
    Returns function of MOJAVE beamed LF with arguments (L, z), where L -
    luminosity in W/Hz and z - redshift.
    """

    delta_min = 1./gamma2
    delta_max = gamma2+np.sqrt(gamma2**2-1.0)

    # Domain where this LF is valid: obs. lum. L1 < L < L2, z1 < z < z2
    if L1 is None:
        L1 = delta_min**p*l1
    if L2 is None:
        L2 = delta_max**p*l2

    # Determine normalizing constant C
    C = (k+1.)/(gamma2**(k+1)-gamma1**(k+1))
    from scipy import special
    from scipy.integrate import quad

    def beta_paper(z, a, b):
        return special.betainc(a, b, z)*special.beta(a, b)

    def d1(l):
        return min(delta_max, max(delta_min, (l/l2)**(1./p)))

    def d2(l):
        return max(delta_min, min(delta_max, (l/l1)**(1./p)))

    def f(delta):
        if 1./gamma2 < delta < 1./gamma1:
            return 1./delta
        elif 1./gamma1 < delta < gamma1*(1.+beta(gamma1)):
            return gamma1
        elif gamma1*(1.+beta(gamma1)) < delta < delta_max:
            return (1.+delta**2)/(2.*delta)
        else:
            return gamma2

    def P_delta(delta):
        return C/(2.*delta**2)*(beta_paper(1.-1./gamma2**2, 0.5, -0.5*k)-
                                beta_paper(1.-1./f(delta)**2, 0.5, -0.5*k))

    def f_D(z):
        return z**m*np.exp(-0.5*((z-z0)/sigma)**2)

    def LF(L, z):
        if not L1 < L < L2:
            return 0
        return n0/L_star*(L/L_star)**alpha*f_D(z)*\
               quad(lambda x: P_delta(x)*x**(-1.0*p*(alpha+1.0)), d1(L), d2(L))[
                   0]

    return LF


def get_number_of_sources(zmin=0.04, zmax=4, n_z=20, F_min=1.5):
    """
    Return number of sources brighter then ``F_min`` per steradian.

    :param zmin: (optional)
        Minimal redhsift. (default: ``0.04``)
    :param zmax: (optional)
        Maximum redhsift. (default: ``4``)
    :param n_z: (optional)
        NUmber of bins across z. TUse larger values for more precise estimate.
        (default: ``20``)
    :param F_min: (optional)
        Minimal observed Flux [Jy]. (default: ``1.5``)

    :return:
        NUmber of sources brighter the ``F_min`` per steradian.
    """
    cosmo = LambdaCDM(H0=70, Om0=0.3, Ode0=0.7)

    lf = beamed_LF()

    def lf_(l, z):
        return cosmo.differential_comoving_volume(z).value*lf(l, z)

    n_zs = list()
    zs = np.linspace(zmin, zmax, n_z)
    dz = zs[1]-zs[0]
    for z in zs:
        dl_cm = cosmo.luminosity_distance(z).value*10**6*pc
        # -23 to convert Jy to erg/Hz/s and -7 to make it W/Hz
        lum_min = 4.*np.pi*dl_cm**2*F_min*(1+z)**(-1)*10**(-23)*10**(-7)
        l1 = max(10**22.2, lum_min)
        n_z = quad(lf_, l1, 1.0*10**29.1, args=(z,))[0]
        print("at z={} n={}".format(z, n_z*dz))
        n_zs.append(n_z)

    return sum(n_zs)*dz


if __name__ == "__main__":
    n = get_number_of_sources(F_min=0.001)
    print(n)