Repository: StingraySoftware/HENDRICS

Issue:
HENaccelsearch and --delta-z argument
Command HENaccelsearch used with fractional value of --delta-z gives the following error: 

HENaccelsearch --fmin 0.33 --fmax 0.37 --delta-z 0.5 Mur19_xmm_epn_det04_ev.nc --nproc 4 --interbin
/home/gianluca/stingray/lib64/python3.6/site-packages/hendrics/efsearch.py:1681: UserWarning: The accelsearch functionality is experimental. Use with care,  and feel free to report any issues.
  "The accelsearch functionality is experimental. Use with care, "
usage: HENaccelsearch [-h] [--outfile OUTFILE] [--emin EMIN] [--emax EMAX]
                      [--fmin FMIN] [--fmax FMAX] [--nproc NPROC]
                      [--zmax ZMAX] [--delta-z DELTA_Z] [--interbin]
                      [--pad-to-double] [--loglevel LOGLEVEL] [--debug]
                      fname
HENaccelsearch: error: argument --delta-z: invalid int value: '0.5'


Discussion:

Patch:
diff --git a/docs/scripts/cli.rst b/docs/scripts/cli.rst
index b582d67..89b89d7 100644
--- a/docs/scripts/cli.rst
+++ b/docs/scripts/cli.rst
@@ -48,7 +48,7 @@ HENaccelsearch
       --fmax FMAX          Maximum frequency to search, in Hz
       --nproc NPROC        Number of processors to use
       --zmax ZMAX          Maximum acceleration (in spectral bins)
-      --delta-z DELTA_Z    Fdot step for search
+      --delta-z DELTA_Z    Fdot step for search (1 is the default resolution)
       --interbin           Use interbinning
       --pad-to-double      Pad to the double of bins (sort-of interbinning)
       --loglevel LOGLEVEL  use given logging level (one between INFO, WARNING,
@@ -753,8 +753,8 @@ HENreadevents
     usage: HENreadevents [-h] [--noclobber] [-g] [--discard-calibration]
                          [-l LENGTH_SPLIT] [--min-length MIN_LENGTH]
                          [--gti-string GTI_STRING] [--randomize-by RANDOMIZE_BY]
-                         [-o OUTFILE] [--loglevel LOGLEVEL] [--debug]
-                         [--nproc NPROC]
+                         [--additional ADDITIONAL [ADDITIONAL ...]] [-o OUTFILE]
+                         [--loglevel LOGLEVEL] [--debug] [--nproc NPROC]
                          files [files ...]
 
     Read a cleaned event files and saves the relevant information in a standard
@@ -778,6 +778,8 @@ HENreadevents
       --randomize-by RANDOMIZE_BY
                             Randomize event arrival times by this amount (e.g. it
                             might be the 0.073-s frame time in XMM)
+      --additional ADDITIONAL [ADDITIONAL ...]
+                            Additional columns to be read from the FITS file
       -o OUTFILE, --outfile OUTFILE
                             Output file
       --loglevel LOGLEVEL   use given logging level (one between INFO, WARNING,
@@ -928,7 +930,7 @@ HENsumfspec
       -h, --help            show this help message and exit
       -o OUTNAME, --outname OUTNAME
                             Output file name for summed (C)PDS. Default:
-                            tot_{c,}pds.{nc,p}
+                            tot_(c)pds.p
 
 
 HENvarenergy
diff --git a/hendrics/efsearch.py b/hendrics/efsearch.py
index b19290d..79de737 100644
--- a/hendrics/efsearch.py
+++ b/hendrics/efsearch.py
@@ -1486,7 +1486,7 @@ def _common_main(args, func):
             emax=args.emax,
             mjdref=mjdref,
             pepoch=mjdref + ref_time / 86400,
-            oversample=args.oversample
+            oversample=args.oversample,
         )
         efperiodogram.upperlim = pf_from_ssig(np.max(stats), events.time.size)
 
@@ -1722,7 +1722,10 @@ def main_accelsearch(args=None):
         help="Maximum acceleration (in spectral bins)",
     )
     parser.add_argument(
-        "--delta-z", default=1, type=int, help="Fdot step for search"
+        "--delta-z",
+        default=1,
+        type=float,
+        help="Fdot step for search (1 is the default resolution)",
     )
     parser.add_argument(
         "--interbin",
diff --git a/hendrics/io.py b/hendrics/io.py
index f07cfe8..e31c75d 100644
--- a/hendrics/io.py
+++ b/hendrics/io.py
@@ -841,7 +841,7 @@ def _split_high_precision_number(varname, var, probesize):
             var_log10 = np.floor(np.log10(dum))
 
         var = np.asarray(var) / (10.0 ** var_log10)
-        var_I = np.floor(var).astype(np.long)
+        var_I = np.floor(var).astype(int)
         var_F = np.array(var - var_I, dtype=np.double)
     else:
         if np.abs(var) < 1 and np.abs(var) > 0.0:
diff --git a/hendrics/phaseogram.py b/hendrics/phaseogram.py
index 9adbcde..bab778b 100644
--- a/hendrics/phaseogram.py
+++ b/hendrics/phaseogram.py
@@ -148,8 +148,13 @@ class BasePhaseogram(object):
         self.fddot = fddot
         self.nt = nt
         self.nph = nph
-        self.mjdref = mjdref
+        if mjdref is None:
+            warnings.warn("MJDREF not set. All MJD values will be incorrect.")
+        self.mjdref = assign_value_if_none(mjdref, 0.)
+        self.ev_times = np.asarray(ev_times)
         self.gti = gti
+        if gti is None:
+            self.gti = np.array([[ev_times[0], ev_times[-1]]])
         self.label = label
         self.test = test
 
@@ -158,7 +163,6 @@ class BasePhaseogram(object):
         self.time_corr = assign_value_if_none(
             time_corr, np.zeros_like(ev_times)
         )
-        self.ev_times = ev_times
         self.freq = freq
         self.norm = norm
         self.position = position
@@ -172,7 +176,7 @@ class BasePhaseogram(object):
             fill_value="extrapolate",
         )
         self.time_corr_mjd_fun = interp1d(
-            self.ev_times / 86400 + mjdref,
+            self.ev_times / 86400 + self.mjdref,
             self.time_corr / 86400,
             bounds_error=False,
             fill_value="extrapolate",
@@ -202,7 +206,7 @@ class BasePhaseogram(object):
             fdot=fdot,
             fddot=fddot,
             plot=False,
-            pepoch=pepoch,
+            pepoch=self.pepoch,
         )
 
         self.phaseogr, phases, times, additional_info = normalized_phaseogram(
@@ -215,7 +219,7 @@ class BasePhaseogram(object):
             fdot=fdot,
             fddot=fddot,
             plot=False,
-            pepoch=pepoch,
+            pepoch=self.pepoch,
         )
 
         self.phases, self.times = phases, times
@@ -233,13 +237,13 @@ class BasePhaseogram(object):
         ax.set_xlabel("Phase")
 
         def s2d(x):
-            return (x - pepoch) / 86400
+            return (x - self.pepoch) / 86400
 
         def d2s(x):
-            return (x - mjdref) * 86400
+            return (x - self.mjdref) * 86400
 
         secax = ax.secondary_yaxis("left", functions=(s2d, d2s))
-        secax.set_ylabel(f"d from MJD {pepoch / 86400 + mjdref}")
+        secax.set_ylabel(f"d from MJD {self.pepoch / 86400 + self.mjdref}")
 
         plt.setp(ax.get_yticklabels(), visible=False)
         ax.set_yticks([])
@@ -515,7 +519,10 @@ class BasePhaseogram(object):
             tm_string += "T0(MET)        {}\n".format(self.t0)
             tm_string += "PB(s)          {}\n".format(self.orbital_period)
 
-        tm_string += "PEPOCH(MET)    {}\n".format(self.pepoch)
+        tm_string += "# PEPOCH(MET)  {}\n".format(self.pepoch)
+        start, stop = self.gti.min(), self.gti.max()
+        tm_string += "START         {}\n".format(start / 86400 + self.mjdref)
+        tm_string += "FINISH        {}\n".format(stop / 86400 + self.mjdref)
         return tm_string
 
 
@@ -725,15 +732,13 @@ class BinaryPhaseogram(BasePhaseogram):
         self.asini = 0
         self.t0 = None
 
-        if "orbital_period" in kwargs:
-            self.orbital_period = kwargs["orbital_period"]
-            kwargs.pop("orbital_period")
-        if "asini" in kwargs:
-            self.asini = kwargs["asini"]
-            kwargs.pop("asini")
-        if "t0" in kwargs:
-            self.t0 = kwargs["t0"]
-            kwargs.pop("t0")
+        try:
+            self.orbital_period = kwargs.pop("orbital_period")
+            self.asini = kwargs.pop("asini")
+            self.t0 = kwargs.pop("t0")
+        except KeyError:
+            raise RuntimeError("Please specify all binary parameters")
+
         BasePhaseogram.__init__(self, *args, **kwargs)
 
     def _construct_widgets(self):
diff --git a/hendrics/read_events.py b/hendrics/read_events.py
index 73764b1..b9e0f24 100644
--- a/hendrics/read_events.py
+++ b/hendrics/read_events.py
@@ -25,7 +25,7 @@ def treat_event_file(
     length_split=None,
     randomize_by=None,
     discard_calibration=False,
-    additional_columns=None
+    additional_columns=None,
 ):
     """Read data from an event file, with no external GTI information.
 
@@ -52,8 +52,12 @@ def treat_event_file(
     # gtistring = assign_value_if_none(gtistring, "GTI,GTI0,STDGTI")
     log.info("Opening %s" % filename)
     try:
-        events = EventList.read(filename, format_="hea", gtistring=gtistring,
-                                additional_columns=additional_columns)
+        events = EventList.read(
+            filename,
+            format_="hea",
+            gtistring=gtistring,
+            additional_columns=additional_columns,
+        )
     except TypeError:  # pragma: no cover
         evtdata = load_events_and_gtis(filename, gtistring=gtistring)
         events = evtdata.ev_list


Test patch:
diff --git a/hendrics/tests/test_efsearch.py b/hendrics/tests/test_efsearch.py
index 98f1bc8..60940cf 100644
--- a/hendrics/tests/test_efsearch.py
+++ b/hendrics/tests/test_efsearch.py
@@ -433,7 +433,7 @@ class TestEFsearch:
                 "--fast",
                 "--find-candidates",
                 "--oversample",
-                "4"
+                "4",
             ]
         )
         outfile = "events_Z22_9.85-9.95Hz_fast" + HEN_FILE_EXTENSION
@@ -478,16 +478,7 @@ class TestEFsearch:
         evfile = self.dum
         with pytest.warns(UserWarning) as record:
             _ = main_zsearch(
-                [
-                    evfile,
-                    "-f",
-                    "9.85",
-                    "-F",
-                    "9.95",
-                    "-n",
-                    "2",
-                    "--fast",
-                ]
+                [evfile, "-f", "9.85", "-F", "9.95", "-n", "2", "--fast",]
             )
         assert np.any(
             [
@@ -714,7 +705,17 @@ class TestEFsearch:
         evfile = self.dum
         with pytest.warns(UserWarning) as record:
             outfile = main_accelsearch(
-                [evfile, "--fmin", "1", "--fmax", "10", "--zmax", "1"]
+                [
+                    evfile,
+                    "--fmin",
+                    "1",
+                    "--fmax",
+                    "10",
+                    "--zmax",
+                    "1",
+                    "--delta-z",
+                    "0.5",
+                ]
             )
         assert np.any(
             [
diff --git a/hendrics/tests/test_fspec.py b/hendrics/tests/test_fspec.py
index 4c5ec5a..0b4ffa1 100644
--- a/hendrics/tests/test_fspec.py
+++ b/hendrics/tests/test_fspec.py
@@ -262,14 +262,12 @@ class TestFullRun(object):
 
     def test_cpds_wrong_norm(self):
         """Test CPDS production."""
-        command = (
-            "{0} {1} -f 128 --save-dyn -k CPDS --norm blablabla -o {2}".format(
-                os.path.join(self.datadir, "monol_testA_E3-50_lc")
-                + HEN_FILE_EXTENSION,
-                os.path.join(self.datadir, "monol_testB_E3-50_lc")
-                + HEN_FILE_EXTENSION,
-                os.path.join(self.datadir, "monol_test_E3-50"),
-            )
+        command = "{0} {1} -f 128 --save-dyn -k CPDS --norm blablabla -o {2}".format(
+            os.path.join(self.datadir, "monol_testA_E3-50_lc")
+            + HEN_FILE_EXTENSION,
+            os.path.join(self.datadir, "monol_testB_E3-50_lc")
+            + HEN_FILE_EXTENSION,
+            os.path.join(self.datadir, "monol_test_E3-50"),
         )
         with pytest.warns(UserWarning) as record:
             hen.fspec.main(command.split())
diff --git a/hendrics/tests/test_phaseogram.py b/hendrics/tests/test_phaseogram.py
index 427431f..929127e 100644
--- a/hendrics/tests/test_phaseogram.py
+++ b/hendrics/tests/test_phaseogram.py
@@ -5,6 +5,7 @@ import numpy as np
 from hendrics.io import save_events, HEN_FILE_EXTENSION, load_folding
 from hendrics.efsearch import main_zsearch
 from hendrics.phaseogram import main_phaseogram, run_interactive_phaseogram
+from hendrics.phaseogram import InteractivePhaseogram, BinaryPhaseogram
 from hendrics.base import hen_root
 from hendrics.fold import HAS_PINT
 from hendrics.plot import plot_folding
@@ -113,6 +114,19 @@ class TestPhaseogram:
             ]
         )
 
+    def test_phaseogram_all_defaults(self):
+        times = [0, 3., 4.]
+        with pytest.warns(UserWarning) as record:
+            _ = InteractivePhaseogram(times, 1., test=True)
+        assert np.any(["MJDREF not set." in r.message.args[0] for r in record])
+
+    def test_binary_all_no_orb(self):
+        times = [0, 3., 4.]
+        with pytest.raises(RuntimeError) as excinfo:
+            # Missing t0
+            _ = BinaryPhaseogram(times, 1., orbital_period=3., asini=3.)
+        assert "Please specify all binary parameters" in str(excinfo.value)
+
     def test_phaseogram_input_norm_invalid(self):
         evfile = self.dum
         with pytest.warns(UserWarning) as record: