Source code for common.sky.cat_builder

"""
Created on Feb. 25, 2019.

Functions to read an astrophysical catalogue (e.g. Swift/BAT), compute count rate through the ARF and save infos in a fits file that will be read by the simulator.

@author: Coleiro Alexis, Univ. Paris Diderot, APC

"""

import numpy as np
from astropy.io import fits
import common.sky.catalog as cata
import astropy.units as u
from common.instru.model_effect import ECLAIRsDetectorEffectDefault

    

def catalog_with_channel(cat_filename, dpix, verbose=False):
    """return a catalog with mean flux in channel 
    
    :param cat_filename: path to fits catalog
    :type cat_filename: string
    :param dpix: object calibration ECLAIRs
    :type dpix: 
    :param verbose:
    :type verbose:
    """
    assert isinstance(dpix, ECLAIRsDetectorEffectDefault)    
    f_fits = fits.open(cat_filename)
    sources_cat = f_fits[1].data
    f_fits.close()
    # in keV. Energy band of the Swift/BAT catalogue. Assumed to be always the same.
    band_bat = np.array([14., 195.])
    spectral_index = sources_cat['GAMMA']
    #(given in 1e-12 erg/s/cm^2 in the Swift/BAT catalogue)
    flux_integrated_kev  = ((sources_cat['FLUX']*1e-12)*u.erg).to('keV').value
    catalog_of_sources = cata.CatalogAstroWithEnergySpecSampling(dpix.chan_center)    
    for idx_src in range(len(sources_cat[1:])):
        spec_index = spectral_index[idx_src]
        src_name = sources_cat['COUNTERPART_NAME'][idx_src]
        fact_indice = (2 - spec_index)
        if np.isclose(spec_index, 2, atol=1e-3):
            phi0 = flux_integrated_kev[idx_src]/np.log(band_bat[1]/band_bat[0])            
        else:
            band_power = np.power(band_bat, 2 - spec_index)                
            phi0 = fact_indice*flux_integrated_kev[idx_src]/(band_power[1] - band_power[0])
        #spec_index
        fact_indice = (1 - spec_index)
        if np.isclose(spec_index, 1, atol=1e-3):
            boundary_log = np.log(dpix.chan_boundary)
            flux_in_channel = phi0*(boundary_log[1:] - boundary_log[:-1])
            print(f"{src_name} closed to 1 {spec_index}")          
        else:
            boundary_pow = np.power(dpix.chan_boundary, 1 - spec_index)      
            flux_in_channel = phi0*(boundary_pow[1:] - boundary_pow[:-1])/fact_indice            
        #print(src_name)
        try:
            assert (flux_in_channel > 0).all()
        except AssertionError:
            print("src name ", src_name, spec_index)
            raise AssertionError
            
        my_src = (sources_cat['RA'][idx_src], sources_cat['DEC'][idx_src], 
                  flux_in_channel, src_name)
        catalog_of_sources.add_src(my_src)
    return catalog_of_sources