Binary systems

The class ixpeobssim.srcmodel.roi.xBinarySource allows to simulate a point-like periodic source in a binary system. The basic ingredients to describe the system are the same as in the simple periodic source configuration but here the input can be in the more generic form of a parameter file including the system ephemeris. The spin and orbital ephemeris are handled by the ixpeobssim.srcmodel.roi.xOrbitalEphemeris class.

The toy_binary configuration file distributed with ixpeobssim describes a pulsar in a binary system with a power-law photon spectrum with normalization that varies as a function of the pulse phase. To illustrate the flexibility of the code, the polarization degree also varies as a function of both energy and pulsar phase.

import os
import numpy

from ixpeobssim.config import file_path_to_model_name
from ixpeobssim.utils.matplotlib_ import plt, setup_gca
from ixpeobssim.srcmodel.roi import xROIModel, xBinarySource
from ixpeobssim.srcmodel.spectrum import power_law
from ixpeobssim.srcmodel.polarization import constant
from ixpeobssim.utils.fmtaxis import fmtaxis
from ixpeobssim import IXPEOBSSIM_SRCMODEL
from ixpeobssim.srcmodel.ephemeris import xOrbitalEphemeris

__model__ = file_path_to_model_name(__file__)
ra, dec = 18.14, -35.02

def pl_norm(phase):
    """Energy spectrum: power-law normalization as a function of the pulse
    phase.
    """
    return 1.25 + numpy.cos(4 * numpy.pi * phase)

pl_index = 2.
spec = power_law(pl_norm, pl_index)

def pol_deg(E, phase, ra=None, dec=None):
    """Polarization degree as a function of the dynamical variables.

    Since we're dealing with a point source the sky direction (ra, dec) is
    irrelevant and, as they are not used, defaulting the corresponding arguments
    to None allows to call the function passing the energy and phase only.
    """
    norm = numpy.clip(E / 10., 0., 1.)
    return norm * (0.5 + 0.25 * numpy.cos(4 * numpy.pi * (phase - 0.25)))

pol_ang = constant(numpy.radians(65.))
file_path = os.path.join(IXPEOBSSIM_SRCMODEL, 'parfiles', 'SAXJ1808.4-3658.par')
ephemeris = xOrbitalEphemeris.from_file(file_path)

src = xBinarySource('Periodic source', ra, dec, spec, pol_deg, pol_ang,
                           ephemeris)

ROI_MODEL = xROIModel(ra, dec, src)

def display(emin=1., emax=12.):
    """Display the source model.
    """
    energy = numpy.linspace(emin, emax, 100)
    phase = numpy.linspace(0., 1., 100)

    # Pulse profile: power-law normalization.
    plt.figure('%s PL normalization' % __model__)
    plt.plot(phase, pl_norm(phase))
    setup_gca(ymin=0., ymax=2.5, **fmtaxis.pp_pl_norm)

    # Pulse profile: polarization degree.
    plt.figure('%s polarization degree' % __model__)
    for E in [2., 5., 8.]:
        plt.plot(phase, pol_deg(E, phase), label='Energy = %.2f keV' % E)
    setup_gca(ymin=0., ymax=1., legend=True, **fmtaxis.pp_pol_deg)

    # Energy spectrum at different phase values.
    plt.figure('%s spectrum' % __model__)
    for p in [0.00, 0.05, 0.10, 0.15, 0.20, 0.25]:
        plt.plot(energy, spec(energy, p), label='Phase = %.2f' % p)
    setup_gca(xmin=emin, xmax=emax, logx=True, logy=True, legend=True,
              grids=True, **fmtaxis.spec)


if __name__ == '__main__':
    from ixpeobssim.config import bootstrap_display
    bootstrap_display()