
""" convert the Labbe+ individual sources into a GSMF """



import numpy as np
from astropy.table import Table, Column, vstack
import astropy.units as u
from astropy.cosmology import Planck18 as cosmo
from flare.stats import poisson_confidence_interval

study_name = 'labbe22'

data = Table.read(f'../../../../Individual/obs/{study_name}.ecsv')

print(data)

output = {'z': [], 'dz': [], 'log10Mstar': [], 'dlog10Mstar': [], 'N': [], 'phi': [], 'phi_err_low': [], 'phi_err_upp': []}

for (z1, z2) in [(8., 9.), (9., 11.)]:

    z = 0.5*(z1+z2)
    dz = z2-z1

    vol = cosmo.comoving_volume(z2) - cosmo.comoving_volume(z1)

    area = 40.*u.arcmin*u.arcmin

    vol_survey = vol*area.to('sr')/(4*np.pi*u.sr)


    dlog10Mstar = 1.0

    for log10Mstar in [10.5]:

        output['z'].append(z)
        output['dz'].append(dz)
        output['log10Mstar'].append(log10Mstar)
        output['dlog10Mstar'].append(dlog10Mstar)

        s = (data['z']>z1)&(data['z']<z2)&(np.fabs(data['log10Mstar']-log10Mstar)<dlog10Mstar/2.)

        N = np.sum(s)
        output['N'].append(N)

        phi = N/vol_survey.to('Mpc^3').value

        output['phi'].append(phi)

        CI = poisson_confidence_interval(N)

        output['phi_err_low'].append(phi - CI[0]/vol_survey.to('Mpc^3').value)
        output['phi_err_upp'].append(CI[1]/vol_survey.to('Mpc^3').value - phi)




table = Table()
table.add_column(Column(data = output['z'], name = 'z'))
table.add_column(Column(data = output['dz'], name = 'deltaz'))
table.add_column(Column(data = output['log10Mstar'], name = 'log10Mstar', unit = 'dex(solMass)'))
table.add_column(Column(data = output['dlog10Mstar'], name = 'deltalog10Mstar'))
table.add_column(Column(data = output['N'], name = 'N'))
table.add_column(Column(data = output['phi'], name = 'phi', unit = '1 / (dex Mpc3)'))
table.add_column(Column(data = output['phi_err_low'], name = 'phi_err_low', unit = '1 / (dex Mpc3)'))
table.add_column(Column(data = output['phi_err_upp'], name = 'phi_err_upp', unit = '1 / (dex Mpc3)'))

table.meta['x'] = 'log10Mstar'
table.meta['y'] = 'phi'
table.meta['name'] = 'Labbe et al. (2022)'
table.meta['type'] = 'binned'
table.meta['redshifts'] = list(set(table['z'].data))
table.meta['references'] = ['2022arXiv220712446L']
table.meta['notes'] = 'Generated by Stephen Wilkins based on 100% completeness'
table.write(f'../binned/{study_name}.ecsv', overwrite = True)
