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 numpy as np


def flux(P, z, alpha):
    return (P * (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)
        cat = np.loadtxt(f).T
        z = cat[1]
        P = cat[0] * 1E-7  # ergs -> Watt

        # TEMPORARY FIX
        P *= (1400 / 200)**2

        S200 = flux(P, z, alpha)
        S200 = S200 / 1E-26  # Watts / m2 / Hz-> Jy
        S154 = S200 * (154 / 200)**alpha

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

    data = np.array([
        np.radians(ras),
        np.radians(np.array(decs) - 27),  # EoR0 field center
        [154E6] * len(ras),               # Reference frequency
        S154s,                            # Flux @ ref freq
        [-1] * len(ras)                   # Spectral index
    ]).T  # Has structure [[ra, dec, refref, flux, alpha]]
    np.save('catalogue-web.npy', data)