luminosity_density.py

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

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

redshift_slices = [0.01, 0.02, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3]
volumes = []

for i, label in enumerate(redshift_slices[:-1]):
    if i == 0:
        z_min = 0
    else:
        z_min = redshift_slices[i]

    z_max = redshift_slices[i + 1]

    print(label, z_min, z_max)

    ratio = 1 * np.pi / 129600
    volume = (Planck15.comoving_volume(z_max) - Planck15.comoving_volume(z_min))
    volume *= ratio
    print(volume.to(units.Mpc**3))
    volumes.append(volume.to(units.Mpc**3).value)

volumes = np.array(volumes)

plt.figure()
for label in [0, 1, 2, 4, 5]: # Skip 3 as it is broken according to Franco
    L_totals = []
    for z in redshift_slices[:-1]:
        print("Processing:", label, z)
        data = np.loadtxt('%d/cone_1x1_z%gtxt_web%d.dat' % (label, z, label))
        L_totals.append(data[:, 0].sum())

    plt.scatter(redshift_slices[:-1], np.array(L_totals) / volumes, label='Model %d' % label)

plt.yscale('log')
plt.savefig('luminosity_per_zslice.png', dpi=600, transparent=True)
plt.show()