Repository: StingraySoftware/stingray
Pull request #225 (2017-07-26 10:31:48)

Added the possibility of deorbiting event lists (with nanograv/PINT)
+ Adds an optional dependency to nanograv/pint

+ Uses it to deorbit files if a pulsar parameter file is provided


Patch:
diff --git a/.travis.yml b/.travis.yml
index 654ec5a3..0209d94c 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -44,6 +44,11 @@ env:
 matrix:
     fast_finish: true
     include:
+        - python: 3.5
+          env: PIP_DEPENDENCIES='git+https://git@github.com/nanograv/pint.git jplephem emcee statsmodels corner' SETUP_CMD='test --remote-data --coverage'
+
+        - python: 2.7
+          env: PIP_DEPENDENCIES='git+https://git@github.com/nanograv/pint.git jplephem emcee statsmodels corner' SETUP_CMD='test --remote-data'
 
         #Try without importing h5py but importing numba
         - python: 2.7
@@ -148,4 +153,4 @@ after_success:
     # If coveralls.io is set up for this package, uncomment the line
     # below and replace "packagename" with the name of your package.
     # The coveragerc file may be customized as needed for your package.
-    - if [[ $SETUP_CMD == 'test --coverage' ]]; then coveralls --rcfile='stingray/tests/coveragerc'; fi
+    - if [[ $SETUP_CMD == *'--coverage'* ]]; then coveralls --rcfile='stingray/tests/coveragerc'; fi
diff --git a/setup.py b/setup.py
index 5a54e696..5d0d10ba 100755
--- a/setup.py
+++ b/setup.py
@@ -76,6 +76,7 @@ package_info = get_package_info()
 
 # Add the project-global data
 package_info['package_data'].setdefault(PACKAGENAME, [])
+package_info['package_data'][PACKAGENAME].append('*/tests/data/*')
 package_info['package_data'][PACKAGENAME].append('data/*')
 package_info['package_data'][PACKAGENAME].append('tests/data/*')
 package_info['package_data'][PACKAGENAME].append('datasets/*')
diff --git a/stingray/pulse/pulsar.py b/stingray/pulse/pulsar.py
index acea44cb..1d05151a 100644
--- a/stingray/pulse/pulsar.py
+++ b/stingray/pulse/pulsar.py
@@ -6,6 +6,12 @@ import numpy as np
 import collections
 from ..utils import simon, jit
 from scipy.optimize import minimize, basinhopping, curve_fit
+try:
+    import pint.toa as toa
+    import pint
+    HAS_PINT = True
+except ImportError:
+    HAS_PINT = False
 
 
 __all__ = ['pulse_phase', 'phase_exposure', 'fold_events', 'stat',
@@ -585,3 +591,57 @@ def get_TOA(prof, period, tstart, template=None, additional_phase=0,
     toa = tstart + phase_res * period
     toaerr = phase_res_err * period
     return toa, toaerr
+
+
+def get_orbital_correction_from_ephemeris_file(mjdstart, mjdstop, parfile,
+                                               ntimes=1000, ephem="DE405"):
+    """Get a correction for orbital motion from pulsar parameter file.
+
+    Parameters
+    ----------
+    mjdstart, mjdstop : float
+        Start and end of the time interval where we want the orbital solution
+    parfile : str
+        Any parameter file understood by PINT (Tempo or Tempo2 format)
+
+    Other parameters
+    ----------------
+    ntimes : int
+        Number of time intervals to use for interpolation. Default 1000
+
+    Returns
+    -------
+    correction_sec : function
+        Function that accepts in input an array of times in seconds and a
+        floating-point MJDref value, and returns the deorbited times
+    correction_mjd : function
+        Function that accepts times in MJDs and returns the deorbited times.
+    """
+    from scipy.interpolate import interp1d
+    simon("Assuming events are already referred to the solar system "
+          "barycenter (timescale is TDB)")
+    if not HAS_PINT:
+        raise ImportError("You need the optional dependency PINT to use this "
+                          "functionality: github.com/nanograv/pint")
+
+    mjds = np.linspace(mjdstart, mjdstop, ntimes)
+    toalist = [None] * len(mjds)
+    for i, m in enumerate(mjds):
+        toalist[i] = toa.TOA(m, obs='Barycenter', scale='tdb')
+
+    toalist = toa.TOAs(toalist = toalist)
+    if 'tdb' not in toalist.table.colnames:
+        toalist.compute_TDBs()
+    if 'ssb_obs_pos' not in toalist.table.colnames:
+        toalist.compute_posvels(ephem, False)
+
+    m = pint.models.get_model(parfile)
+    tdbtimes = m.get_barycentric_toas(toalist.table)
+
+    correction_mjd = interp1d(mjds, tdbtimes)
+
+    def correction_sec(times, mjdref):
+        deorb_mjds = correction_mjd(times / 86400 + mjdref)
+        return np.array((deorb_mjds - mjdref) * 86400)
+
+    return correction_sec, correction_mjd
\ No newline at end of file


Test patch:
diff --git a/stingray/pulse/tests/data/example_pint.par b/stingray/pulse/tests/data/example_pint.par
new file mode 100644
index 00000000..54d50d0b
--- /dev/null
+++ b/stingray/pulse/tests/data/example_pint.par
@@ -0,0 +1,149 @@
+PSRJ           1855+09
+RAJ             18:57:36.3932884         1  0.00002602730280675029
+DECJ           +09:43:17.29196           1  0.00078789485676919773
+F0             186.49408156698235146     1  0.00000000000698911818
+F1             -6.2049547277487420583e-16 1  1.7380934373573401505e-20
+PEPOCH         49453
+POSEPOCH       49453
+DMEPOCH        49453
+DM             13.29709
+PMRA           -2.5054345161030380639    1  0.03104958261053317181
+PMDEC          -5.4974558631993817232    1  0.06348008663748286318
+PX             1.2288569063263405232     1  0.21243361289239687251
+SINI           0.99741717335200923866    1  0.00182023515130851988
+BINARY         DD
+PB             12.327171194774200418     1  0.00000000079493185824
+T0             49452.940695077335647     1  0.00169031830532837251
+A1             9.2307804312998001928     1  0.00000036890718667634
+OM             276.55142180589701234     1  0.04936551005019605698
+ECC            2.1745265668236919017e-05 1  0.00000004027191312623
+M2             0.26111312480723428917    1  0.02616161008932908066
+START          53358.726464889485214
+FINISH         55108.922917417192366
+TZRMJD         54177.508359343262555
+TZRFRQ         424
+TZRSITE        ao
+TRES           0.395
+EPHVER         5
+DMX_0001       0
+DMX_0002       0.00011110286020705287208 1  0.00004673570459450457
+DMX_0003       -4.4655555822498924595e-05 1  0.00004202379265494554
+DMX_0004       -3.1723662424549210086e-05 1  0.00004222071152827488
+DMX_0005       -2.6615937544541523929e-05 1  0.00003676165288001985
+DMX_0006       7.1454456174372302607e-05 1  0.00004374167394142721
+DMX_0007       7.7436702748396629494e-06 1  0.00005131995439632110
+DMX_0008       6.6283145761418474104e-05 1  0.00003755691800563188
+DMX_0009       9.9608002224418395616e-05 1  0.00003619457239625261
+DMX_0010       0.00021384397943417332134 1  0.00003970891234128656
+DMX_0011       0.00016549149164987531535 1  0.00004615618107309406
+DMX_0012       0.00025014821251322403243 1  0.00003617384078038137
+DMX_0013       0.0003254428344575823765  1  0.00003525448646108801
+DMX_0014       0.00070209493875510397669 1  0.00003567637835272617
+DMX_0015       0.00083029061827721808596 1  0.00003306863035583663
+DMX_0016       0.00094566920797150631083 1  0.00004382957805166126
+DMX_0017       0.0010188311608163321283  1  0.00003607409255560491
+DMX_0018       0.001089116586071231421   1  0.00004813789826343072
+DMX_0019       0.0010386415586231195997  1  0.00004612675176906664
+DMX_0020       0.0013195672245089193941  1  0.00004616167182274142
+DMX_0021       0.0012154222610258824143  1  0.00004294821352263614
+DMX_0022       0.0013377609996199927778  1  0.00004221028001939009
+DMX_0023       0.0016163009324016204224  1  0.00005156225843276498
+DMX_0024       0.0016669838878444672594  1  0.00006395311039024736
+DMX_0025       0.00045683635569106320759 1  0.00004912260448808334
+DMX_0026       0.00052771780683627349866 1  0.00004152704579342970
+DMX_0027       0.00072894864454871476041 1  0.00004234732022848593
+DMX_0028       0.0011420785526064839561  1  0.00005990586419386911
+DMX_0029       0.0011050194788476920155  1  0.00004148068669700297
+DMX_0030       0.001691787104074557445   1  0.00006464664732840410
+DMXR1_0001     53358.7273
+DMXR1_0002     53420.553999999999998
+DMXR1_0003     53448.4787
+DMXR1_0004     53477.401299999999999
+DMXR1_0005     53532.232600000000001
+DMXR1_0006     53603.036100000000001
+DMXR1_0007     53628.9651
+DMXR1_0008     53686.7974
+DMXR1_0009     53715.735700000000001
+DMXR1_0010     53750.622499999999999
+DMXR1_0011     53798.5061
+DMXR1_0012     53851.371600000000001
+DMXR1_0013     53891.252300000000002
+DMXR1_0014     54043.841899999999999
+DMXR1_0015     54092.717700000000001
+DMXR1_0016     54135.582599999999999
+DMXR1_0017     54519.529199999999999
+DMXR1_0018     54569.4146
+DMXR1_0019     54678.100000000000001
+DMXR1_0020     54819.719799999999999
+DMXR1_0021     54862.5884
+DMXR1_0022     54925.432400000000001
+DMXR1_0023     54981.2806
+DMXR1_0024     54998.209199999999999
+DMXR1_0025     53926.1631
+DMXR1_0026     53968.0623
+DMXR1_0027     54008.953499999999998
+DMXR1_0028     54177.478000000000002
+DMXR1_0029     54472.662600000000001
+DMXR1_0030     55108.9038
+DMXR2_0001     53358.773299999999999
+DMXR2_0002     53420.581099999999999
+DMXR2_0003     53448.514500000000002
+DMXR2_0004     53477.429500000000001
+DMXR2_0005     53532.2715
+DMXR2_0006     53603.084599999999998
+DMXR2_0007     53628.983499999999999
+DMXR2_0008     53686.837
+DMXR2_0009     53715.761700000000001
+DMXR2_0010     53750.662199999999999
+DMXR2_0011     53798.534099999999999
+DMXR2_0012     53851.3992
+DMXR2_0013     53891.280100000000001
+DMXR2_0014     54043.8674
+DMXR2_0015     54092.742599999999999
+DMXR2_0016     54135.6066
+DMXR2_0017     54519.5562
+DMXR2_0018     54569.466899999999999
+DMXR2_0019     54678.128299999999999
+DMXR2_0020     54819.746200000000002
+DMXR2_0021     54862.6145
+DMXR2_0022     54925.461300000000001
+DMXR2_0023     54981.325199999999999
+DMXR2_0024     54998.226800000000001
+DMXR2_0025     53926.192399999999999
+DMXR2_0026     53968.0901
+DMXR2_0027     54008.986599999999999
+DMXR2_0028     54177.508600000000001
+DMXR2_0029     54472.6889
+DMXR2_0030     55108.9221
+CLK            TT(TAI)
+MODE 1
+UNITS          TDB
+T2CMETHOD      TEMPO
+NE_SW          0.000
+CORRECT_TROPOSPHERE  N
+EPHEM          DE405
+NITS           1
+NTOA           702
+TIMEEPH        FB90
+CHI2R          2.1896 637
+JUMP -chanid asp_424 7.6456527699426e-07 1
+JUMP -chanid asp_428 1.6049580247793e-06 1
+JUMP -chanid asp_432 2.677273589972e-06 1
+JUMP -chanid asp_436 3.0814628244857e-06 1
+JUMP -chanid asp_440 4.136918843609e-06 1
+JUMP -chanid asp_1382 -5.9069311360427e-05 1
+JUMP -chanid asp_1386 -5.9582397512386e-05 1
+JUMP -chanid asp_1390 -5.9756103518774e-05 1
+JUMP -chanid asp_1394 -5.9347801609886e-05 1
+JUMP -chanid asp_1398 -5.9009509604032e-05 1
+JUMP -chanid asp_1402 -5.886527352183e-05 1
+JUMP -chanid asp_1406 -5.9015609469307e-05 1
+JUMP -chanid asp_1410 -5.8771703168481e-05 1
+JUMP -chanid asp_1414 -5.8756443357489e-05 1
+JUMP -chanid asp_1418 -5.8960922044628e-05 1
+JUMP -chanid asp_1422 -5.8883992500888e-05 1
+JUMP -chanid asp_1426 -5.8918166650158e-05 1
+JUMP -chanid asp_1430 -5.8867635126778e-05 1
+JUMP -chanid asp_1434 -5.9028297376503e-05 1
+JUMP -chanid asp_1438 -5.9000827072833e-05 1
+JUMP -chanid asp_1442 -5.8537946239722e-05 1
diff --git a/stingray/pulse/tests/data/example_pint.tim b/stingray/pulse/tests/data/example_pint.tim
new file mode 100644
index 00000000..c4af4d3d
--- /dev/null
+++ b/stingray/pulse/tests/data/example_pint.tim
@@ -0,0 +1,4 @@
+FORMAT 1
+53358.000056.3.000.000.tsum 420.000 53358.7731394424088 1.196 ao -fe 430G -be ASP -B 430 -bw 4.0 -tobs 1806.363 -pta NANOGrav -proc dfg+12 -chanid asp_420
+53420.000025.3.000.000.tsum 420.000 53420.5541580125302 1.835 ao -fe 430G -be ASP -B 430 -bw 4.0 -tobs 1806.363 -pta NANOGrav -proc dfg+12 -chanid asp_420
+53448.000028.3.000.000.tsum 420.000 53448.4789265955255 1.675 ao -fe 430G -be ASP -B 430 -bw 4.0 -tobs 1806.363 -pta NANOGrav -proc dfg+12 -chanid asp_420
diff --git a/stingray/pulse/tests/test_pulse.py b/stingray/pulse/tests/test_pulse.py
index 87e646cf..48144d21 100644
--- a/stingray/pulse/tests/test_pulse.py
+++ b/stingray/pulse/tests/test_pulse.py
@@ -1,17 +1,75 @@
 from __future__ import division, print_function, absolute_import
+
 import numpy as np
 from stingray.pulse.pulsar import fold_events, get_TOA
 from stingray.pulse.pulsar import stat, z_n, pulse_phase, phase_exposure
 from stingray.pulse.pulsar import fold_detection_level, z2_n_detection_level
 from stingray.pulse.pulsar import fold_profile_probability, z2_n_probability
+from stingray.pulse.pulsar import get_orbital_correction_from_ephemeris_file
+from ..pulsar import HAS_PINT
+from astropy.tests.helper import remote_data
+import pytest
+import os
 
 def _template_fun(phase, ph0, amplitude, baseline=0):
     return baseline + amplitude * np.cos((phase - ph0) * 2 * np.pi)
 
 
 class TestAll(object):
-
     """Unit tests for the stingray.pulsar module."""
+    @classmethod
+    def setup_class(cls):
+        cls.curdir = os.path.abspath(os.path.dirname(__file__))
+        cls.datadir = os.path.join(cls.curdir, 'data')
+
+    @remote_data
+    @pytest.mark.skipif('not HAS_PINT')
+    def test_pint_installed_correctly(self):
+        import pint.toa as toa
+        from pint.residuals import resids
+        import pint.models.model_builder as mb
+        import astropy.units as u
+        parfile = os.path.join(self.datadir, 'example_pint.par')
+        timfile = os.path.join(self.datadir, 'example_pint.tim')
+
+        toas = toa.get_TOAs(timfile, ephem="DE405",
+                            planets=False, include_bipm=False)
+        model = mb.get_model(parfile)
+
+        pint_resids_us = resids(toas, model, False).time_resids.to(u.s)
+
+        # Due to the gps2utc clock correction. We are at 3e-8 seconds level.
+        assert np.all(np.abs(pint_resids_us.value) < 3e-6)
+
+    @remote_data
+    @pytest.mark.skipif('not HAS_PINT')
+    def test_orbit_from_parfile(self):
+        import pint.toa as toa
+        parfile = os.path.join(self.datadir, 'example_pint.par')
+        timfile = os.path.join(self.datadir, 'example_pint.tim')
+
+        toas = toa.get_TOAs(timfile, ephem="DE405",
+                            planets=False, include_bipm=False)
+        mjds = np.array([m.value for m in toas.get_mjds(high_precision=True)])
+
+        mjdstart, mjdstop = mjds[0] - 1, mjds[-1] + 1
+
+        correction_sec, correction_mjd = \
+            get_orbital_correction_from_ephemeris_file(mjdstart, mjdstop,
+                                                       parfile, ntimes=1000)
+
+        mjdref = 50000
+        toa_sec = (mjds - mjdref) * 86400
+        corr = correction_mjd(mjds)
+        corr_s = correction_sec (toa_sec, mjdref)
+        assert np.allclose(corr, corr_s / 86400 + mjdref)
+
+    @pytest.mark.skipif('HAS_PINT')
+    def test_orbit_from_parfile_raises(self):
+        print("Doesn't have pint")
+        with pytest.raises(ImportError):
+            get_orbital_correction_from_ephemeris_file(0, 0,
+                                                       parfile='ciaociao')
 
     def test_stat(self):
         """Test pulse phase calculation, frequency only."""