Source code for soxs.background.point_sources

import numpy as np
from astropy.io import ascii
from astropy.table import Table
from scipy.interpolate import InterpolatedUnivariateSpline
from scipy.special import erf

from soxs.constants import erg_per_keV, keV_per_erg
from soxs.simput import SimputCatalog, SimputPhotonList
from soxs.spatial import construct_wcs
from soxs.spectra import get_tbabs_absorb, get_wabs_absorb
from soxs.utils import mylog, parse_prng, parse_value, soxs_cfg

# Function for computing spectral index of AGN sources
# a fit to the data from Figure 13a of Hickox & Markevitch 2006
# http://adsabs.harvard.edu/abs/2006ApJ...645...95H

# Parameters

aa = -14.0
bb = 0.5
cc = 0.5
dd = 1.8


# Here x = log10(flux)
def get_agn_index(x):
    y = (x - aa) / bb
    return cc * erf(y) + dd


# Index for galaxies
gal_index = 2.0

fb_emin = 0.5  # keV, low energy bound for the logN-logS flux band
fb_emax = 2.0  # keV, high energy bound for the logN-logS flux band

spec_emin = 0.1  # keV, minimum energy of mock spectrum
spec_emax = 10.0  # keV, max energy of mock spectrum


def get_flux_scale(ind, fb_emin, fb_emax, spec_emin, spec_emax):
    f_g = np.log(spec_emax / spec_emin) * np.ones(ind.size)
    f_E = np.log(fb_emax / fb_emin) * np.ones(ind.size)
    n1 = ind != 1.0
    n2 = ind != 2.0
    f_g[n1] = (spec_emax ** (1.0 - ind[n1]) - spec_emin ** (1.0 - ind[n1])) / (
        1.0 - ind[n1]
    )
    f_E[n2] = (fb_emax ** (2.0 - ind[n2]) - fb_emin ** (2.0 - ind[n2])) / (
        2.0 - ind[n2]
    )
    fscale = f_g / f_E
    return fscale


def generate_fluxes(fov, prng):
    from soxs.data import cdf_agn, cdf_fluxes, cdf_gal

    prng = parse_prng(prng)

    fov = parse_value(fov, "arcmin")

    logf = np.log10(cdf_fluxes)

    n_gal = np.rint(cdf_gal[-1])
    n_agn = np.rint(cdf_agn[-1])
    F_gal = cdf_gal / cdf_gal[-1]
    F_agn = cdf_agn / cdf_agn[-1]
    f_gal = InterpolatedUnivariateSpline(F_gal, logf)
    f_agn = InterpolatedUnivariateSpline(F_agn, logf)

    fov_area = fov**2

    n_gal = int(n_gal * fov_area / 3600.0)
    n_agn = int(n_agn * fov_area / 3600.0)
    mylog.debug("%s AGN, %s galaxies in the FOV.", n_agn, n_gal)

    randvec1 = prng.uniform(size=n_agn)
    agn_fluxes = 10 ** f_agn(randvec1)

    randvec2 = prng.uniform(size=n_gal)
    gal_fluxes = 10 ** f_gal(randvec2)

    return agn_fluxes, gal_fluxes


def generate_positions(num, fov, sky_center, prng):
    fov *= 60.0  # convert to arcsec
    x = prng.uniform(low=-0.5 * fov, high=0.5 * fov, size=num)
    y = prng.uniform(low=-0.5 * fov, high=0.5 * fov, size=num)
    w = construct_wcs(*sky_center)
    ra0, dec0 = w.wcs_pix2world(x, y, 1)
    return ra0, dec0


def generate_sources(fov, sky_center, prng=None):
    r"""
    Make a catalog of point sources.

    Parameters
    ----------
    fov : float, (value, unit) tuple, or :class:`~astropy.units.Quantity`
        The field of view in arcminutes.
    sky_center : array-like
        The center RA, Dec of the field of view in degrees.
    prng : :class:`~numpy.random.RandomState` object, integer, or None
        A pseudo-random number generator. Typically will only
        be specified if you have a reason to generate the same
        set of random numbers, such as for a test. Default is None,
        which sets the seed based on the system time.
    """
    prng = parse_prng(prng)

    fov = parse_value(fov, "arcmin")

    agn_fluxes, gal_fluxes = generate_fluxes(fov, prng)

    fluxes = np.concatenate([agn_fluxes, gal_fluxes])

    ind = np.concatenate(
        [get_agn_index(np.log10(agn_fluxes)), gal_index * np.ones(gal_fluxes.size)]
    )

    ra0, dec0 = generate_positions(fluxes.size, fov, sky_center, prng)

    return ra0, dec0, fluxes, ind


def make_ptsrc_background(
    exp_time,
    fov,
    sky_center,
    absorb_model=None,
    nH=None,
    area=40000.0,
    input_sources=None,
    output_sources=None,
    dump_fluxes_band=None,
    diffuse_unresolved=True,
    prng=None,
):
    r"""
    Make a point-source background.

    Parameters
    ----------
    exp_time : float, (value, unit) tuple, or :class:`~astropy.units.Quantity`
        The exposure time of the observation in seconds.
    fov : float, (value, unit) tuple, or :class:`~astropy.units.Quantity`
        The field of view in arcminutes.
    sky_center : array-like
        The center RA, Dec of the field of view in degrees.
    absorb_model : string, optional
        The absorption model to use, "wabs" or "tbabs".
        Defaults to the value in the SOXS configuration file.
    nH : float, (value, unit) tuple, or :class:`~astropy.units.Quantity`, optional
        The hydrogen column in units of 10**22 atoms/cm**2.
        Defaults to the value in the SOXS configuration file.
    area : float, (value, unit) tuple, or :class:`~astropy.units.Quantity`, optional
        The effective area in cm**2. It must be large enough
        so that a sufficiently large sample is drawn for the
        ARF. Default: 40000.
    input_sources : string, optional
        If set to a filename, input the source positions, fluxes,
        and spectral indices from an ASCII table instead of generating
        them. Default: None
    output_sources : string, optional
        If set to a filename, output the properties of the sources
        within the field of view to a file. Default: None
    diffuse_unresolved : boolean, optional
        Add a diffuse component across the entire field of view to represent
        the unresolved flux from sources at very small fluxes. Default: True
    prng : :class:`~numpy.random.RandomState` object, integer, or None
        A pseudo-random number generator. Typically will only
        be specified if you have a reason to generate the same
        set of random numbers, such as for a test. Default is None,
        which sets the seed based on the system time.
    """
    prng = parse_prng(prng)

    exp_time = parse_value(exp_time, "s")
    fov = parse_value(fov, "arcmin")
    if nH is not None:
        nH = parse_value(nH, "1.0e22*cm**-2")
    area = parse_value(area, "cm**2")
    if input_sources is None:
        ra0, dec0, fluxes, ind = generate_sources(fov, sky_center, prng=prng)
    else:
        mylog.info("Reading in point-source properties from %s.", input_sources)
        t = ascii.read(input_sources)
        ra0 = t["RA"].data
        dec0 = t["Dec"].data
        fluxes = t["flux_0.5_2.0_keV"].data
        ind = t["index"].data
    num_sources = fluxes.size

    mylog.debug("Generating spectra from %d sources.", num_sources)

    # If requested, output the source properties to a file
    if output_sources is not None:
        t = Table(
            [ra0, dec0, fluxes, ind], names=("RA", "Dec", "flux_0.5_2.0_keV", "index")
        )
        t["RA"].unit = "deg"
        t["Dec"].unit = "deg"
        t["flux_0.5_2.0_keV"].unit = "erg/(cm**2*s)"
        t["index"].unit = ""
        t.write(output_sources, format="ascii.ecsv", overwrite=True)

    # Add in the diffuse component for completely unresolved sources
    if diffuse_unresolved:
        """
        We add a diffuse background across the field of view, from completely
        unresolved sources, to match the fluxes in Table 3 of Hickox & Markevitch
        2006, ApJ, 645, 95. The first number "F12" is in erg/s/cm**2/deg**2, and
        is the difference between the total flux in sources below 5e-17 erg/s/cm**2
        and the "average" value in Table 3 for the unresolved flux in the 1-2 keV
        band of 1.04e-12 erg/s/cm**2/deg**2. The factor of 2.0 comes from the
        scaling between the fluxes in the 1-2 keV and 0.5-2 keV bands, assuming
        a spectral index of 2. The whole thing finally needs to be scaled by the
        FOV at the end.
        """
        F12 = 0.676e-12  # erg/s/cm**2/deg**2 in 1-2 keV band
        diffuse_flux = 2.0 * F12 * (fov / 60) ** 2
        fluxes = np.append(fluxes, diffuse_flux)
        ind = np.append(ind, 2.0)

    # Pre-calculate for optimization
    eratio = spec_emax / spec_emin
    oma = 1.0 - ind
    invoma = 1.0 / oma
    invoma[oma == 0.0] = 1.0
    fac1 = spec_emin**oma
    fac2 = spec_emax**oma - fac1
    fluxscale = get_flux_scale(ind, fb_emin, fb_emax, spec_emin, spec_emax)

    # Using the energy flux, determine the photon flux by simple scaling
    ref_ph_flux = fluxes * fluxscale * keV_per_erg
    # Now determine the number of photons we will generate
    n_photons = prng.poisson(ref_ph_flux * exp_time * area)

    all_energies = []
    all_ra = []
    all_dec = []
    detected = []

    for i, nph in enumerate(n_photons):
        if nph > 0:
            # Generate the energies in the source frame
            u = prng.uniform(size=nph)
            if ind[i] == 1.0:
                energies = spec_emin * (eratio**u)
            else:
                energies = fac1[i] + u * fac2[i]
                if invoma[i] == -1.0:
                    energies = 1.0 / energies
                else:
                    energies **= invoma[i]

            # Assign positions for this source
            if i < num_sources:
                ra = ra0[i] * np.ones(nph)
                dec = dec0[i] * np.ones(nph)
            else:
                # this is the diffuse source--fill the whole FOV
                ra, dec = generate_positions(nph, fov, sky_center, prng)

            all_energies.append(energies)
            all_ra.append(ra)
            all_dec.append(dec)
            detected.append(i)

    mylog.debug("Finished generating spectra.")

    if dump_fluxes_band is not None:
        dfluxes = np.zeros_like(ref_ph_flux)
        for idx, e in zip(detected, all_energies):
            eidxs = (e > dump_fluxes_band[0]) & (e < dump_fluxes_band[1])
            dfluxes[idx] = np.sum(e[eidxs]) * erg_per_keV / area / exp_time
        filename = f"dump_fluxes_{dump_fluxes_band[0]}_{dump_fluxes_band[1]}.dat"
        np.savetxt(filename, dfluxes)

    all_energies = np.concatenate(all_energies)
    all_ra = np.concatenate(all_ra)
    all_dec = np.concatenate(all_dec)

    # Remove some photons due to Galactic foreground absorption.
    # We will throw a lot of stuff away, but this is more general and still
    # faster.
    if nH is None:
        nH = float(soxs_cfg.get("soxs", "bkgnd_nH"))
    if nH > 0.0:
        if absorb_model is None:
            absorb_model = soxs_cfg.get("soxs", "bkgnd_absorb_model")
        if absorb_model == "wabs":
            absorb = get_wabs_absorb(all_energies, nH)
        elif absorb_model == "tbabs":
            absorb = get_tbabs_absorb(all_energies, nH)
        randvec = prng.uniform(size=all_energies.size)
        all_energies = all_energies[randvec < absorb]
        all_ra = all_ra[randvec < absorb]
        all_dec = all_dec[randvec < absorb]
        all_nph = all_energies.size
        mylog.debug("%d photons remain after foreground galactic absorption.", all_nph)

    all_flux = np.sum(all_energies) * erg_per_keV / (exp_time * area)

    output_events = {
        "ra": all_ra,
        "dec": all_dec,
        "energy": all_energies,
        "flux": all_flux,
    }

    return output_events


[docs] def make_point_sources_file( filename, name, exp_time, fov, sky_center, absorb_model=None, nH=None, area=40000.0, append=False, overwrite=False, src_filename=None, input_sources=None, output_sources=None, diffuse_unresolved=True, prng=None, ): """ Make a SIMPUT catalog made up of contributions from point sources. Parameters ---------- filename : string The filename for the SIMPUT catalog. name : string The name of the SIMPUT photon list. exp_time : float, (value, unit) tuple, or :class:`~astropy.units.Quantity` The exposure time of the observation in seconds. fov : float, (value, unit) tuple, or :class:`~astropy.units.Quantity` The field of view in arcminutes. sky_center : array-like The center RA, Dec of the field of view in degrees. absorb_model : string, optional The absorption model to use, "wabs" or "tbabs". Defaults to the value in the SOXS configuration file. nH : float, (value, unit) tuple, or :class:`~astropy.units.Quantity`, optional The hydrogen column in units of 10**22 atoms/cm**2. Defaults to the value in the SOXS configuration file. area : float, (value, unit) tuple, or :class:`~astropy.units.Quantity`, optional The effective area in cm**2. It must be large enough so that a sufficiently large sample is drawn for the ARF. Default: 40000. append : boolean, optional If True, the photon list source will be appended to an existing SIMPUT catalog. Default: False overwrite : boolean, optional Set to True to overwrite previous files. Default: False src_filename : string, optional If set, this will be the filename to write the source to. By default, the source will be written to the same file as the SIMPUT catalog. input_sources : string, optional If set to a filename, input the source positions, fluxes, and spectral indices from an ASCII table instead of generating them. Default: None output_sources : string, optional If set to a filename, output the properties of the sources within the field of view to a file. Default: None diffuse_unresolved : boolean, optional Add a diffuse component across the entire field of view to represent the unresolved flux from sources at very small fluxes. Default: True prng : :class:`~numpy.random.RandomState` object, integer, or None A pseudo-random number generator. Typically will only be specified if you have a reason to generate the same set of random numbers, such as for a test. Default is None, which sets the seed based on the system time. """ events = make_ptsrc_background( exp_time, fov, sky_center, absorb_model=absorb_model, nH=nH, area=area, input_sources=input_sources, output_sources=output_sources, prng=prng, diffuse_unresolved=diffuse_unresolved, ) phlist = SimputPhotonList( events["ra"], events["dec"], events["energy"], events["flux"], name=name ) if append: cat = SimputCatalog.from_file(filename) cat.append(phlist, src_filename=src_filename, overwrite=overwrite) else: cat = SimputCatalog.from_source( filename, phlist, src_filename=src_filename, overwrite=overwrite ) return cat
[docs] def make_point_source_list(output_file, fov, sky_center, prng=None): r""" Make a list of point source properties and write it to an ASCII table file. Parameters ---------- output_file : string The ASCII table file to write the source properties to. fov : float, (value, unit) tuple, or :class:`~astropy.units.Quantity` The field of view in arcminutes. sky_center : array-like The center RA, Dec of the field of view in degrees. prng : :class:`~numpy.random.RandomState` object, integer, or None A pseudo-random number generator. Typically will only be specified if you have a reason to generate the same set of random numbers, such as for a test. Default is None, which sets the seed based on the system time. """ ra0, dec0, fluxes, ind = generate_sources(fov, sky_center, prng=prng) t = Table( [ra0, dec0, fluxes, ind], names=("RA", "Dec", "flux_0.5_2.0_keV", "index") ) t["RA"].unit = "deg" t["Dec"].unit = "deg" t["flux_0.5_2.0_keV"].unit = "erg/(cm**2*s)" t["index"].unit = "" t.write(output_file, format="ascii.ecsv", overwrite=True)