web2flux.py

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

import argparse

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


parser = argparse.ArgumentParser()
parser.add_argument('datafiles', nargs='+')
args = parser.parse_args()


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, S155s, alphas, redshifts = [], [], [], [], []
    for f in args.datafiles:
        print(f)

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

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

        ras.extend(cat[2])
        decs.extend(cat[3])
        S155s.extend(S155)
        redshifts.extend(z)

    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
        S155s,               # Flux @ reference freq
        [-1] * len(ras),     # Spectral index
        redshifts,           # Redshifts
    ]).T
    np.save('catalogue-web.npy', data)

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