web2flux.py

#! /usr/bin/env python
from __future__ import print_function, division

import glob

from astropy.cosmology import Planck15
import astropy.units as units
import matplotlib.pyplot as plt
import numpy as np

from astrobits.transforms import radec_to_cartesian, cartesian_to_radec, rotate_about_y


def flux(L, z, alpha):
    """
    L = 4 pi S D^2 / (1 + z)^1+alpha
    """
    return (L * (1 + z)**(1 + alpha)) / (4 * np.pi * Planck15.luminosity_distance(z).to_value(units.meter)**2)


try:
    data = np.load('catalogue-web.npy')
except Exception:
    alpha = -1

    ras, decs, S154s = [], [], []
    for f in glob.glob('*.dat'):
        print(f)
        redshift_slice = float(f[10:][:-11])
        if redshift_slice > 0.3:
            continue

        cat = np.loadtxt(f).T
        z = cat[1]
        L200 = cat[0] / 1E7  # Luminosity @ 200 MHz (ergs / s / Hz -> Watts / Hz)

        # TEMPORARY FIX
        # To fix bug in Franco's code. Has been fixed for redshfit slices <= 0.3
        if redshift_slice > 0.3:
            P *= (1400 / 200)**2
            P /= 1E6

        S200 = flux(L200, z, alpha)          # Flux @ 200 MHz
        S200 = S200 * 1E26                # Watts / m^2 / Hz -> Jy
        S154 = S200 * (154 / 200)**alpha  # Flux @ 154 MHz

        ras.extend(cat[2])
        decs.extend(cat[3])
        S154s.extend(S154)

    ras, decs = np.radians(ras), np.radians(decs)
    x, y, z = radec_to_cartesian(ras, decs)
    x, y, z = rotate_about_y(x, y, z, np.radians(27))  # Rotate into EoR0 field
    ras, decs = cartesian_to_radec(x, y, z)

    data = np.array([
        ras,                 # Ra
        decs,                # Dec
        [154E6] * len(ras),  # Reference freq
        S154s,               # Flux @ reference freq
        [-1] * len(ras)      # Spectral index
    ]).T
    np.save('catalogue-web.npy', data)

print("Total flux across the field:", data[:, 3].sum())

plt.hist(np.log10(data[:, 3]), bins=100, log=True)
plt.show()

plt.scatter(data[:, 0], data[:, 1], s=3, alpha=0.1, c=np.log10(data[:, 3]))
plt.show()