# Licensed under a 3-clause BSD style license - see LICENSE.rst
from __future__ import absolute_import, division, print_function, unicode_literals
import numpy as np
from scipy.interpolate import InterpolatedUnivariateSpline
from astropy.utils import lazyproperty
from astropy import units as u
from astropy.table import Table
from ..utils.scripts import make_path
from ..utils.time import time_ref_from_dict
from ..utils.fitting import Parameter, Parameters, Model
__all__ = ["PhaseCurveTableModel", "LightCurveTableModel"]
[docs]class PhaseCurveTableModel(Model):
r"""Temporal phase curve model.
Phase for a given time is computed as:
.. math::
\phi(t) = \phi_0 + f_0(t-t_0) + (1/2)f_1(t-t_0)^2 + (1/6)f_2(t-t_0)^3
Strictly periodic sources such as gamma-ray binaries have ``f1=0`` and ``f2=0``.
Sources like some pulsars where the period spins up or down have ``f1!=0``
and / or ``f2 !=0``. For a binary, ``f0`` should be calculated as 1/T,
where T is the period of the binary in unit of ``seconds``.
The "phase curve", i.e. multiplicative flux factor for a given phase is given
by a `~astropy.table.Table` of nodes ``(phase, norm)``, using linear interpolation
and circular behaviour, where ``norm(phase=0) == norm(phase=1)``.
Parameters
----------
table : `~astropy.table.Table`
A table of 'PHASE' vs 'NORM' should be given
time_0 : float
The MJD value where phase is considered as 0.
phase_0 : float
Phase at the reference MJD
f0, f1, f2 : float
Derivatives of the function phi with time of order 1, 2, 3
in units of ``s^-1, s^-2 & s^-3``, respectively.
Examples
--------
Create an example phase curve object::
from astropy.table import Table
from gammapy.utils.scripts import make_path
from gammapy.time.models import PhaseCurveTableModel
filename = make_path('$GAMMAPY_DATA/tests/phasecurve_LSI_DC.fits')
table = Table.read(str(filename))
phase_curve = PhaseCurveTableModel(table, time_0=43366.275, phase_0=0.0, f0=4.367575e-7, f1=0.0, f2=0.0)
Use it to compute a phase and evaluate the phase curve model for a given time:
>>> phase_curve.phase(time=46300.0)
0.7066006737999402
>>> phase_curve.evaluate_norm_at_time(46300)
0.49059393580053845
"""
def __init__(self, table, time_0, phase_0, f0, f1, f2):
self.table = table
self.parameters = Parameters(
[
Parameter("time_0", time_0),
Parameter("phase_0", phase_0),
Parameter("f0", f0),
Parameter("f1", f1),
Parameter("f2", f2),
]
)
[docs] def phase(self, time):
"""Evaluate phase for a given time.
Parameters
----------
time : array_like
Returns
-------
phase : array_like
"""
pars = self.parameters
time_0 = pars["time_0"].value
phase_0 = pars["phase_0"].value
f0 = pars["f0"].value
f1 = pars["f1"].value
f2 = pars["f2"].value
t = (time - time_0) * u.day.to(u.second)
phase = self._evaluate_phase(t, phase_0, f0, f1, f2)
return np.remainder(phase, 1)
@staticmethod
def _evaluate_phase(t, phase_0, f0, f1, f2):
return phase_0 + t * (f0 + t * (f1 / 2 + f2 / 6 * t))
[docs] def evaluate_norm_at_time(self, time):
"""Evaluate for a given time.
Parameters
----------
time : array_like
Time since the ``reference`` time.
Returns
-------
norm : array_like
"""
phase = self.phase(time)
return self.evaluate_norm_at_phase(phase)
[docs] def evaluate_norm_at_phase(self, phase):
xp = self.table["PHASE"]
fp = self.table["NORM"]
return np.interp(x=phase, xp=xp, fp=fp, period=1)
[docs]class LightCurveTableModel(Model):
"""Temporal light curve model.
The lightcurve is given as a table with columns ``time`` and ``norm``.
The ``norm`` is supposed to be a unite-less multiplicative factor in the model,
to be multiplied with a spectral model.
The model does linear interpolation for times between the given ``(time, norm)`` values.
The implementation currently uses `scipy.interpolate.InterpolatedUnivariateSpline`,
using degree ``k=1`` to get linear interpolation.
This class also contains an ``integral`` method, making the computation of
mean fluxes for a given time interval a one-liner.
Parameters
----------
table : `~astropy.table.Table`
A table with 'TIME' vs 'NORM'
Examples
--------
Read an example light curve object:
>>> from gammapy.time.models import LightCurveTableModel
>>> path = '$GAMMAPY_DATA/tests/models/light_curve/lightcrv_PKSB1222+216.fits'
>>> light_curve = LightCurveTableModel.read(path)
Show basic information about the lightcurve:
>>> print(light_curve)
LightCurve model summary:
Start time: 59000.5 MJD
End time: 61862.5 MJD
Norm min: 0.01551196351647377
Norm max: 1.0
Compute ``norm`` at a given time:
>>> light_curve.evaluate_norm_at_time(46300)
0.49059393580053845
Compute mean ``norm`` in a given time interval:
>>> light_curve.mean_norm_in_time_interval(46300, 46301)
"""
def __init__(self, table):
self.table = table
def __str__(self):
ss = "LightCurveTableModel model summary:\n"
ss += "Start time: {} MJD\n".format(self._time[0].mjd)
ss += "End time: {} MJD\n".format(self._time[-1].mjd)
ss += "Norm min: {}\n".format(self.table["NORM"].min())
ss += "Norm max: {}\n".format(self.table["NORM"].max())
return ss
[docs] @classmethod
def read(cls, path):
"""Read lightcurve model table from FITS file.
TODO: This doesn't read the XML part of the model yet.
"""
path = make_path(path)
table = Table.read(str(path))
return cls(table)
@lazyproperty
def _interpolator(self):
x = self.table["TIME"].data
y = self.table["NORM"].data
return InterpolatedUnivariateSpline(x, y, k=1)
@lazyproperty
def _time_ref(self):
return time_ref_from_dict(self.table.meta)
@lazyproperty
def _time(self):
return self._time_ref + self.table["TIME"].data * u.s
[docs] def evaluate_norm_at_time(self, time):
"""Evaluate for a given time.
Parameters
----------
time : array_like
Time since the ``reference`` time.
Returns
-------
norm : array_like
"""
return self._interpolator(time)
[docs] def mean_norm_in_time_interval(self, time_min, time_max):
"""Compute mean ``norm`` in a given time interval.
TODO: vectorise, i.e. allow arrays of time intervals in a single call.
Parameters
----------
time_min, time_max : float
Time interval
Returns
-------
norm : float
Mean norm
"""
dt = time_max - time_min
integral = self._interpolator.integral(time_min, time_max)
return integral / dt