Repository: StingraySoftware/stingray
Pull request #197 (2017-06-19 11:43:38)

Add MJDREF information
In real-life applications, when loading observations from different instruments, the reference MJD will always be needed. It is also useful in order to save the light curves with a clear time reference.

<!-- Reviewable:start -->
---
This change is [<img src="https://reviewable.io/review_button.svg" height="34" align="absmiddle" alt="Reviewable"/>](https://reviewable.io/reviews/stingraysoftware/stingray/197)
<!-- Reviewable:end -->



Patch:
diff --git a/stingray/crossspectrum.py b/stingray/crossspectrum.py
index 04e3c238..5c890dd0 100644
--- a/stingray/crossspectrum.py
+++ b/stingray/crossspectrum.py
@@ -133,6 +133,7 @@ class Crossspectrum(object):
         self.gti = gti
         self.lc1 = lc1
         self.lc2 = lc2
+
         self._make_crossspectrum(lc1, lc2)
 
     def _make_crossspectrum(self, lc1, lc2):
@@ -144,6 +145,9 @@ class Crossspectrum(object):
         if not isinstance(lc2, Lightcurve):
             raise TypeError("lc2 must be a lightcurve.Lightcurve object")
 
+        if self.lc2.mjdref != self.lc1.mjdref:
+            raise ValueError("MJDref is different in the two light curves")
+
         # Then check that GTIs make sense
         if self.gti is None:
             self.gti = cross_two_gtis(lc1.gti, lc2.gti)
diff --git a/stingray/events.py b/stingray/events.py
index 8faeb79d..15df218e 100644
--- a/stingray/events.py
+++ b/stingray/events.py
@@ -122,7 +122,8 @@ class EventList(object):
             tseg = self.gti[-1][1] - tstart
 
         return Lightcurve.make_lightcurve(self.time, dt, tstart=tstart,
-                                          gti=self.gti, tseg=tseg)
+                                          gti=self.gti, tseg=tseg,
+                                          mjdref=self.mjdref)
 
     @staticmethod
     def from_lc(lc):
diff --git a/stingray/lightcurve.py b/stingray/lightcurve.py
index b2354751..788a6a46 100644
--- a/stingray/lightcurve.py
+++ b/stingray/lightcurve.py
@@ -21,7 +21,7 @@ valid_statistics = ["poisson", "gauss", None]
 
 class Lightcurve(object):
     def __init__(self, time, counts, err=None, input_counts=True,
-                 gti=None, err_dist='poisson'):
+                 gti=None, err_dist='poisson', mjdref=0):
         """
         Make a light curve object from an array of time stamps and an
         array of counts.
@@ -62,6 +62,9 @@ class Lightcurve(object):
             uncertainties and other statistical values apropriately.
             Default makes no assumptions and keep errors equal to zero.
 
+        mjdref: float
+            MJD reference (useful in most high-energy mission data)
+
 
         Attributes
         ----------
@@ -98,6 +101,10 @@ class Lightcurve(object):
         dt: float
             The time resolution of the light curve.
 
+        mjdref: float
+            MJD reference date (tstart / 86400 gives the date in MJD at the
+            start of the observation)
+
         tseg: float
             The total duration of the light curve.
 
@@ -155,6 +162,7 @@ class Lightcurve(object):
                       "Sorry for the inconvenience.")
                 err = np.zeros_like(counts)
 
+        self.mjdref = mjdref
         self.time = np.asarray(time)
         self.dt = time[1] - time[0]
 
@@ -194,6 +202,22 @@ class Lightcurve(object):
                                               self.tstart + self.tseg]]))
         check_gtis(self.gti)
 
+    def change_mjdref(self, new_mjdref):
+        """Change the MJDREF of the light curve.
+
+        Times will be now referred to this new MJDREF
+
+        Parameters
+        ----------
+        new_mjdref : float
+            New MJDREF
+        """
+        time_shift = (new_mjdref - self.mjdref) * 86400
+
+        new_lc = self.shift(time_shift)
+        new_lc.mjdref = new_mjdref
+        return new_lc
+
     def shift(self, time_shift):
         """
         Shift the light curve and the GTIs in time.
@@ -204,7 +228,7 @@ class Lightcurve(object):
             The amount of time that the light curve will be shifted
         """
         new_lc = Lightcurve(self.time + time_shift, self.counts,
-                            gti=self.gti + time_shift)
+                            gti=self.gti + time_shift, mjdref=self.mjdref)
         new_lc.countrate = self.countrate
         new_lc.counts = self.counts
         new_lc.counts_err = self.counts_err
@@ -243,6 +267,9 @@ class Lightcurve(object):
             raise ValueError("Time arrays of both light curves must be "
                              "of same dimension and equal.")
 
+        if self.mjdref != other.mjdref:
+            raise ValueError("MJDref is different in the two light curves")
+
         new_counts = np.add(self.counts, other.counts)
 
         if self.err_dist.lower() != other.err_dist.lower():
@@ -261,7 +288,8 @@ class Lightcurve(object):
         common_gti = cross_two_gtis(self.gti, other.gti)
 
         lc_new = Lightcurve(self.time, new_counts,
-                            err=new_counts_err, gti=common_gti)
+                            err=new_counts_err, gti=common_gti,
+                            mjdref=self.mjdref)
 
         return lc_new
 
@@ -297,6 +325,10 @@ class Lightcurve(object):
             raise ValueError("Time arrays of both light curves must be "
                              "of same dimension and equal.")
 
+        if self.mjdref != other.mjdref:
+            raise ValueError("MJDref is different in the two light curves")
+
+
         new_counts = np.subtract(self.counts, other.counts)
 
         if self.err_dist.lower() != other.err_dist.lower():
@@ -315,7 +347,8 @@ class Lightcurve(object):
         common_gti = cross_two_gtis(self.gti, other.gti)
 
         lc_new = Lightcurve(self.time, new_counts,
-                            err=new_counts_err, gti=common_gti)
+                            err=new_counts_err, gti=common_gti,
+                            mjdref=self.mjdref)
 
         return lc_new
 
@@ -338,7 +371,8 @@ class Lightcurve(object):
         array([100, 100, 100])
         """
         lc_new = Lightcurve(self.time, -1*self.counts,
-                            err=self.counts_err, gti=self.gti)
+                            err=self.counts_err, gti=self.gti,
+                            mjdref=self.mjdref)
 
         return lc_new
 
@@ -388,7 +422,7 @@ class Lightcurve(object):
         elif isinstance(index, slice):
             new_counts = self.counts[index.start:index.stop:index.step]
             new_time = self.time[index.start:index.stop:index.step]
-            return Lightcurve(new_time, new_counts)
+            return Lightcurve(new_time, new_counts, mjdref=self.mjdref)
         else:
             raise IndexError("The index must be either an integer or a slice "
                              "object !")
@@ -413,7 +447,7 @@ class Lightcurve(object):
         return baseline
 
     @staticmethod
-    def make_lightcurve(toa, dt, tseg=None, tstart=None, gti=None):
+    def make_lightcurve(toa, dt, tseg=None, tstart=None, gti=None, mjdref=0):
 
         """
         Make a light curve out of photon arrival times.
@@ -482,7 +516,7 @@ class Lightcurve(object):
 
         counts = np.asarray(counts)
 
-        return Lightcurve(time, counts, gti=gti)
+        return Lightcurve(time, counts, gti=gti, mjdref=mjdref)
 
     def rebin(self, dt_new, method='sum'):
         """
@@ -516,7 +550,8 @@ class Lightcurve(object):
             utils.rebin_data(self.time, self.counts, dt_new,
                              yerr=self.counts_err, method=method)
 
-        lc_new = Lightcurve(bin_time, bin_counts, err=bin_err)
+        lc_new = Lightcurve(bin_time, bin_counts, err=bin_err,
+                            mjdref=self.mjdref)
         return lc_new
 
     def join(self, other):
@@ -556,6 +591,9 @@ class Lightcurve(object):
         >>> lc.counts
         array([ 300,  100,  400,  600, 1200,  800])
         """
+        if self.mjdref != other.mjdref:
+            raise ValueError("MJDref is different in the two light curves")
+
         if self.dt != other.dt:
             utils.simon("The two light curves have different bin widths.")
 
@@ -618,7 +656,8 @@ class Lightcurve(object):
 
         gti = join_gtis(self.gti, other.gti)
 
-        lc_new = Lightcurve(new_time, new_counts, err=new_counts_err, gti=gti)
+        lc_new = Lightcurve(new_time, new_counts, err=new_counts_err, gti=gti,
+                            mjdref=self.mjdref)
 
         return lc_new
 
@@ -894,7 +933,8 @@ class Lightcurve(object):
             stop = stop_bins[i]
             # Note: GTIs are consistent with default in this case!
             new_lc = Lightcurve(self.time[start:stop], self.counts[start:stop],
-                                err=self.counts_err[start:stop])
+                                err=self.counts_err[start:stop],
+                                mjdref=self.mjdref)
             list_of_lcs.append(new_lc)
 
         return list_of_lcs
diff --git a/stingray/varenergyspectrum.py b/stingray/varenergyspectrum.py
index 79cc417f..ecf1fa6d 100644
--- a/stingray/varenergyspectrum.py
+++ b/stingray/varenergyspectrum.py
@@ -116,7 +116,8 @@ class VarEnergySpectrum(object):
                                              self.bin_time,
                                              tstart=tstart,
                                              tseg=tstop - tstart,
-                                             gti=gti)
+                                             gti=gti,
+                                             mjdref=self.events1.mjdref)
 
         if exclude:
             ref_intervals = self._decide_ref_intervals(channel_band,
@@ -125,7 +126,7 @@ class VarEnergySpectrum(object):
             ref_intervals = [self.ref_band]
 
         ref_lc = Lightcurve(base_lc.time, np.zeros_like(base_lc.counts),
-                            gti=base_lc.gti)
+                            gti=base_lc.gti, mjdref=base_lc.mjdref)
 
         for i in ref_intervals:
             good = (energies2 >= i[0]) & (energies2 < i[1])
@@ -133,7 +134,8 @@ class VarEnergySpectrum(object):
                                                 self.bin_time,
                                                 tstart=tstart,
                                                 tseg=tstop - tstart,
-                                                gti=gti)
+                                                gti=gti,
+                                                mjdref=self.events2.mjdref)
             ref_lc = ref_lc + new_lc
 
         return base_lc, ref_lc


Test patch:
diff --git a/stingray/tests/test_lightcurve.py b/stingray/tests/test_lightcurve.py
index e7611480..bbff3fb2 100644
--- a/stingray/tests/test_lightcurve.py
+++ b/stingray/tests/test_lightcurve.py
@@ -238,6 +238,12 @@ class TestLightcurve(object):
         lc = lc1 + lc2
         np.testing.assert_almost_equal(lc.gti, self.gti)
 
+    def test_add_with_different_mjdref(self):
+        lc1 = Lightcurve(self.times, self.counts, gti=self.gti, mjdref=57000)
+        lc2 = Lightcurve(self.times, self.counts, gti=self.gti, mjdref=57001)
+        with pytest.raises(ValueError):
+            lc1 + lc2
+
     def test_add_with_different_gtis(self):
         gti = [[0., 3.5]]
         lc1 = Lightcurve(self.times, self.counts, gti=self.gti)
@@ -264,6 +270,7 @@ class TestLightcurve(object):
 
         assert np.all(lc.counts == lc1.counts + lc2.counts)
         assert np.all(lc.countrate == lc1.countrate + lc2.countrate)
+        assert lc1.mjdref == lc.mjdref
 
     def test_sub_with_diff_time_arrays(self):
         _times = [1, 2, 3, 4, 5]
@@ -283,6 +290,12 @@ class TestLightcurve(object):
             lc = lc1 - lc2
             assert "ightcurves have different statistics" in str(w[0].message)
 
+    def test_sub_with_different_mjdref(self):
+        lc1 = Lightcurve(self.times, self.counts, gti=self.gti, mjdref=57000)
+        lc2 = Lightcurve(self.times, self.counts, gti=self.gti, mjdref=57001)
+        with pytest.raises(ValueError):
+            lc1 - lc2
+
     def test_subtraction(self):
         _counts = [3, 4, 5, 6]
 
@@ -293,6 +306,7 @@ class TestLightcurve(object):
 
         expected_counts = np.array([1, 2, 3, 4])
         assert np.all(lc.counts == expected_counts)
+        assert lc1.mjdref == lc.mjdref
 
     def test_negation(self):
         lc = Lightcurve(self.times, self.counts)
@@ -300,6 +314,7 @@ class TestLightcurve(object):
         _lc = lc + (-lc)
 
         assert not np.all(_lc.counts)
+        assert _lc.mjdref == lc.mjdref
 
     def test_len_function(self):
         lc = Lightcurve(self.times, self.counts)
@@ -324,6 +339,8 @@ class TestLightcurve(object):
         assert np.all(lc[:2].counts == np.array([2, 2]))
         assert np.all(lc[2:].counts == np.array([2, 2]))
         assert np.all(lc[:].counts == np.array([2, 2, 2, 2]))
+        assert lc[:].mjdref == lc.mjdref
+
 
     def test_slicing_index_error(self):
         lc = Lightcurve(self.times, self.counts)
@@ -342,6 +359,12 @@ class TestLightcurve(object):
             lc1.join(lc2)
             assert "different bin widths" in str(w[0].message)
 
+    def test_join_with_different_mjdref(self):
+        lc1 = Lightcurve(self.times, self.counts, gti=self.gti, mjdref=57000)
+        lc2 = Lightcurve(self.times, self.counts, gti=self.gti, mjdref=57001)
+        with pytest.raises(ValueError):
+            lc1.join(lc2)
+
     def test_join_disjoint_time_arrays(self):
         _times = [5, 6, 7, 8]
         _counts = [2, 2, 2, 2]
@@ -353,6 +376,7 @@ class TestLightcurve(object):
 
         assert len(lc.counts) == len(lc.time) == 8
         assert np.all(lc.counts == 2)
+        assert lc.mjdref == lc1.mjdref
 
     def test_join_overlapping_time_arrays(self):
         _times = [3, 4, 5, 6]
@@ -375,11 +399,13 @@ class TestLightcurve(object):
         assert np.all(lc1.time == np.array([2, 3, 4]))
         assert np.all(lc1.counts == np.array([2, 2, 2]))
         np.testing.assert_almost_equal(lc1.gti[0][0], 1.5)
+        assert lc1.mjdref == lc.mjdref
 
         lc2 = lc.truncate(stop=2)
         assert np.all(lc2.time == np.array([1, 2]))
         assert np.all(lc2.counts == np.array([2, 2]))
         np.testing.assert_almost_equal(lc2.gti[-1][-1], 2.5)
+        assert lc2.mjdref == lc.mjdref
 
     def test_truncate_by_time_stop_less_than_start(self):
         lc = Lightcurve(self.times, self.counts)
@@ -399,26 +425,31 @@ class TestLightcurve(object):
         assert np.all(lc1.time == np.array([1, 2, 3, 4]))
         assert np.all(lc1.counts == np.array([2, 2, 2, 2]))
         np.testing.assert_almost_equal(lc1.gti[0][0], 0.5)
+        assert lc1.mjdref == lc.mjdref
 
         lc2 = lc.truncate(stop=3, method='time')
         assert np.all(lc2.time == np.array([1, 2]))
         assert np.all(lc2.counts == np.array([2, 2]))
         np.testing.assert_almost_equal(lc2.gti[-1][-1], 2.5)
+        assert lc2.mjdref == lc.mjdref
 
     def test_sort(self):
         _times = [1, 2, 3, 4]
         _counts = [40, 10, 20, 5]
-        lc = Lightcurve(_times, _counts)
+        lc = Lightcurve(_times, _counts, mjdref=57000)
+        mjdref = lc.mjdref
 
         lc.sort()
 
         assert np.all(lc.counts == np.array([5, 10, 20, 40]))
         assert np.all(lc.time == np.array([4, 2, 3, 1]))
+        assert lc.mjdref == mjdref
 
         lc.sort(reverse=True)
 
         assert np.all(lc.counts == np.array([40, 20, 10,  5]))
         assert np.all(lc.time == np.array([1, 3, 2, 4]))
+        assert lc.mjdref == mjdref
 
     def test_plot_matplotlib_not_installed(self):
         try: