diff --git a/eventdisplay_anasum/light_curve_analysis.C b/eventdisplay_anasum/light_curve_analysis.C
index fef766f..112080c 100644
--- a/eventdisplay_anasum/light_curve_analysis.C
+++ b/eventdisplay_anasum/light_curve_analysis.C
@@ -102,16 +102,12 @@ void light_curve_analysis(
      string iAnaSumFile = "anasum/anasum.combined.root",
      string iMJDIntervalFile = "HESSJ0632p057.MJD.v3.txt",
      double iEnergyThreshold = 0.35,
+     double iSpectralIndex = -2.5,
      double iSig_UL = -5., double iSig_minEvents = -10.)
 {
-     // energy threshold in TeV
-     double i_fixed_Emin = iEnergyThreshold;
-     // assumed spectral index:
-     double i_fixed_Index = -2.5;
-
      cout << "# Light curve generation ";
-     cout << "E_min > " << i_fixed_Emin << " TeV; ";
-     cout << "Spectral index " << i_fixed_Index << endl;
+     cout << "E_min > " << iEnergyThreshold << " TeV; ";
+     cout << "Spectral index " << iSpectralIndex << endl;
 
      vector<pair<double, double>> min_max;
      size_t n_intervals = 0;
@@ -138,8 +134,8 @@ void light_curve_analysis(
                min_max[i].second
           );
           a.setSignificanceParameters( iSig_UL, iSig_minEvents, 0.95 );
-          a.setSpectralParameters( i_fixed_Emin, 1., i_fixed_Index );
-          a.calculateIntegralFlux( i_fixed_Emin );
+          a.setSpectralParameters( iEnergyThreshold, 1., iSpectralIndex);
+          a.calculateIntegralFlux( iEnergyThreshold );
           a.printECSVLine();
      }
 }
diff --git a/eventdisplay_anasum/light_curve_analysis.sh b/eventdisplay_anasum/light_curve_analysis.sh
index 71b1615..0760d82 100755
--- a/eventdisplay_anasum/light_curve_analysis.sh
+++ b/eventdisplay_anasum/light_curve_analysis.sh
@@ -2,8 +2,8 @@
 # Light-curve analysis using Eventdisplay container and podman
 #
 
-if [ "$#" -ne 2 ]; then
-    echo "Usage: $0 <anasum_file> <time_interval_file>"
+if [ "$#" -ne 4 ]; then
+    echo "Usage: $0 <anasum_file> <time_interval_file> <energy threshold> <spectral index>"
     echo "  (use 'RUNWISE' for time interval file for run-wise analysis)"
     echo
     exit 1
@@ -22,9 +22,8 @@ light_curves()
 }
 
 echo "Running light-curve analysis..."
-ETRESH="0.3"
 MINSIG="-5."
 MINEVT="-10."
-MINSIG="2."
-MINEVT="2."
-light_curves "$1" "$2" "$ETRESH" "$MINSIG" "$MINEVT"
+# MINSIG="2."
+# MINEVT="2."
+light_curves "$1" "$2" "$3" "$4" "$MINSIG" "$MINEVT"
diff --git a/v2dl5/binaries.py b/v2dl5/binaries.py
index 79575e1..6bbb838 100644
--- a/v2dl5/binaries.py
+++ b/v2dl5/binaries.py
@@ -21,6 +21,7 @@ def binary_properties():
             "name": "HESS J0632+057",
             "orbital_period": 317.3,
             "mjd_0": 54857.0,  # Bongiorno et al 2011
+            "mjd_orbit_count": 52953.0,  # Adams et al 2021
         },
         "LS I +61 303": {
             "name": "LS I +61 303",
diff --git a/v2dl5/light_curves/binary_plotting.py b/v2dl5/light_curves/binary_plotting.py
index 5cb58d1..02fe8d2 100644
--- a/v2dl5/light_curves/binary_plotting.py
+++ b/v2dl5/light_curves/binary_plotting.py
@@ -158,6 +158,9 @@ def _get_number_columns_and_rows(self, number):
             n_columns = 4
             fig_size = (10, 8)
         if number > 16:
+            n_columns = 5
+            fig_size = (16, 10)
+        if number > 21:
             n_columns = 8
             fig_size = (16, 10)
         return n_columns, math.ceil(number / n_columns), fig_size
@@ -236,7 +239,7 @@ def plot_flux_vs_orbit_number(
 
         for i in range(phase_bins):
             axes = plt.subplot(n_rows, n_columns, i + 1)
-            axes.set_xlim([min(orbits), max(orbits)])
+            axes.set_xlim([min(orbits)-1., max(orbits)+1.])
             axes.set_ylim([y_min, y_max])
 
             phase_min = i * (1.0 / phase_bins)
@@ -304,11 +307,11 @@ def _get_light_curve_in_mjd_limits(
         data,
         y_key,
         time_axis,
-        mjd_min,
-        mjd_max,
-        orbit_number,
-        phase_min,
-        phase_max,
+        mjd_min=None,
+        mjd_max=None,
+        orbit_number=None,
+        phase_min=None,
+        phase_max=None,
         min_significance=None,
     ):
         """
@@ -341,6 +344,9 @@ def _get_light_curve_in_mjd_limits(
             Light-curve data as lists.
         """
         x, y, e, x_ul, y_ul = [], [], [], [], []
+        if y_key not in data:
+            self._logger.warning(f"Y-axis {y_key} not found in data")
+            return x, y, e, x_ul, y_ul
         mjd = data["MJD"]
         for i, t in enumerate(mjd):
             if (
@@ -401,6 +407,158 @@ def get_marker_and_color(self, idx):
         )
         return color, marker
 
+    def plot_live_time_vs_phase_bin(
+        self,
+        instrument,
+        phase_bins=10,
+        file_type=".pdf",
+        figure_dir="./figures/",
+    ):
+        """Plot live time histogram vs phase bin."""
+        self._logger.info("Plotting live time histogram vs phase bin")
+
+        ax = plotting_utilities.paper_figures(None, None)
+        ax.set_xlim([0, 1])
+
+        for idx, (instrument, data) in enumerate(self.data.items()):
+            if self.config[idx].get("plot_live_time_histogram", False) is False:
+                continue
+            color, _ = self.get_marker_and_color(idx)
+
+            x, y, _, _, _ = self._get_light_curve_in_mjd_limits(
+                data, "live_time", "orbital phase",
+            )
+            if len(x) == 0:
+                continue
+            x = np.array(x)
+            y = np.array(y)
+            bin_edges = np.linspace(0, 1, phase_bins + 1)
+            bin_width = bin_edges[1] - bin_edges[0]
+            bin_heights, _ = np.histogram(x, bins=bin_edges, weights=y)
+            bin_centers = (bin_edges[:-1] + bin_edges[1:]) / 2
+            plt.bar(
+                bin_centers,
+                bin_heights,
+                width=bin_width*0.9,
+                label=(
+                    instrument
+                    if self.config[idx].get("plot_label") is None
+                    else self.config[idx]["plot_label"]
+                ),
+                color=color,
+                alpha=0.7,
+            )
+
+        plt.xlabel("Phase")
+        plt.ylabel("Live Time (h)")
+
+        plt.legend()
+        plotting_utilities.print_figure(
+            f"Light-Curve-{self.binary['name']}-Live-Time-vs-Phase-Bin",
+            file_type=file_type,
+            figure_dir=figure_dir,
+        )
+
+    def plot_index_vs_flux(
+        self,
+        instrument,
+        file_type=".pdf",
+        figure_dir="./figures/",
+    ):
+        """Plot spectral index vs flux."""
+        self._logger.info("Plotting spectral index vs flux")
+
+        _ = plotting_utilities.paper_figures(None, None)
+
+        for idx, (instrument, data) in enumerate(self.data.items()):
+            if self.config[idx].get("plot_flux_vs_index", False) is False:
+                continue
+            color, _ = self.get_marker_and_color(idx)
+
+            _, x, x_e, _, _ = self._get_light_curve_in_mjd_limits(
+                data, "flux", "orbital phase",
+            )
+            _, y, y_e, _, _ = self._get_light_curve_in_mjd_limits(
+                data, "index", "orbital phase",
+            )
+            if len(x) == 0 or len(y) == 0:
+                continue
+            x = np.array(x)
+            y = np.array(y)
+            plt.errorbar(
+                x,
+                y,
+                xerr=x_e,
+                yerr=y_e,
+                label=(
+                    instrument
+                    if self.config[idx].get("plot_label") is None
+                    else self.config[idx]["plot_label"]
+                ),
+                color=color,
+                linestyle="none",
+                marker='o',
+                alpha=0.7,
+            )
+
+        plt.xlabel(self.config[0].get("flux_axis_label", ""))
+        plt.ylabel(self.config[0].get("index_axis_label", ""))
+
+        plt.legend()
+        plotting_utilities.print_figure(
+            f"Light-Curve-{self.binary['name']}-Flux-vs-Index",
+            file_type=file_type,
+            figure_dir=figure_dir,
+        )
+
+    def plot_distribution(
+        self,
+        instrument,
+        y_axis="flux",
+        file_type=".pdf",
+        figure_dir="./figures/",
+    ):
+        """Plot distribution of y-axis data."""
+        self._logger.info(f"Plotting {y_axis} distribution")
+
+        _ = plotting_utilities.paper_figures(None, None)
+
+        for idx, (instrument, data) in enumerate(self.data.items()):
+            if self.config[idx].get("plot_1d_distribution", False) is False:
+                continue
+            color, _ = self.get_marker_and_color(idx)
+
+            _, y, _, _, _ = self._get_light_curve_in_mjd_limits(
+                data, y_axis, "orbital phase",
+            )
+            if len(y) == 0:
+                continue
+            y = np.array(y)
+            plt.hist(
+                y,
+                bins=25,
+                label=(
+                    instrument
+                    if self.config[idx].get("plot_label") is None
+                    else self.config[idx]["plot_label"]
+                ),
+                color=color,
+                alpha=0.7,
+            )
+
+        plt.xlabel(self.config[0].get(y_axis + "_axis_label", ""))
+        plt.ylabel("Counts")
+
+        plt.legend()
+        plotting_utilities.print_figure(
+            f"Light-Curve-{self.binary['name']}-{y_axis}-Distribution",
+            file_type=file_type,
+            figure_dir=figure_dir,
+        )
+
+
+
+
 
 # Temporary stuff - probably not needed
 #
diff --git a/v2dl5/light_curves/data_reader.py b/v2dl5/light_curves/data_reader.py
index a9cf178..e21b9aa 100644
--- a/v2dl5/light_curves/data_reader.py
+++ b/v2dl5/light_curves/data_reader.py
@@ -70,7 +70,11 @@ def _read_fluxes_from_file(self, data_config):
 
         if data_config["file_name"].endswith((".csv", ".ecsv")):
             self._logger.info("Reading data from %s", data_config["file_name"])
-            return self._read_fluxes_from_ecsv_file(data_config)
+            return self._read_fluxes_from_ecsv_file(
+                file_name=data_config["file_name"],
+                mjd_min=data_config.get("mjd_min", -1.0),
+                mjd_max=data_config.get("mjd_max", -1.0),
+            )
 
         return None
 
@@ -114,7 +118,11 @@ def _add_orbital_parameters(self, data):
         data["phase_err"] = [data["phase_err_low"], data["phase_err_hig"]]
 
         data["orbit_number"] = [
-            orbit.get_orbit_number(mjd=mjd, orbital_period=orbital_period, mjd_0=mjd_0)
+            orbit.get_orbit_number(
+                mjd=mjd,
+                orbital_period=orbital_period,
+                mjd_orbit_count=float(self.binary.get("mjd_orbit_count", self.binary.get("mjd_0")))
+            )
             for mjd in data["MJD"]
         ]
 
@@ -126,15 +134,15 @@ def convert_photon_to_energy_flux(self, c_e, e_0, gamma):
         return [v * f for v in self], [e * f for e in c_e]
 
     def _read_fluxes_from_ecsv_file(
-        self, data_config, time_min_max=True, mjd_min=-1.0, mjd_max=-1.0
+        self, file_name, time_min_max=True, mjd_min=-1.0, mjd_max=-1.0
     ):
         """
         Read gamma-ray fluxes from ecsv file (open gamma-ray format).
 
         Parameters
         ----------
-        data_config: dict
-            configuration dictionary
+        file_name: str, Path
+            Name of the file to read.
         time_min_max: bool
             flag to read time_min and time_max
         mjd_min: float
@@ -143,43 +151,41 @@ def _read_fluxes_from_ecsv_file(
             MJD max value for MJD cut
 
         """
-        table = Table.read(data_config["file_name"])
+        table = Table.read(file_name)
         f = {}
-        try:
-            if not time_min_max:
-                table["time_min"] = table["time"].data
-                table["time_max"] = table["time"].data + 0.1
-
-            # MJD filter
-            condition = np.ones(len(table), dtype=bool)
-            if mjd_min > -1:
-                condition &= table["time_min"] > mjd_min
-            if mjd_max > -1:
-                condition &= table["time_max"] < mjd_max
-            table = table[condition]
-
-            f = {}
-            f["time_min"] = table["time_min"].data.tolist()
-            f["time_max"] = table["time_max"].data.tolist()
-            f["flux"] = table["flux"].data.flatten().tolist()
-            if "flux_err" in table.colnames:
-                f["flux_err"] = table["flux_err"].data.flatten().tolist()
-            else:
-                up = table["flux_up"].data.flatten().tolist()
-                down = table["flux_down"].data.flatten().tolist()
-                f["flux_err"] = [0.5 * abs(u - d) for u, d in zip(up, down)]
-            f["MJD"] = [0.5 * (a + b) for a, b in zip(f["time_min"], f["time_max"])]
-            f["MJD_err"] = [0.5 * (b - a) for a, b in zip(f["time_min"], f["time_max"])]
-            if "flux_ul" in table.colnames:
-                flux_ul = table["flux_ul"].data.flatten().tolist()
-                is_ul = table["is_ul"].data.flatten().tolist()
-                f["flux_ul"] = [flux if is_ul else -1.0 for flux, is_ul in zip(flux_ul, is_ul)]
-            else:
-                f["flux_ul"] = [-1.0 if fe > 0 else fl for fe, fl in zip(f["flux_err"], f["flux"])]
-            for column_name in ["significance", "index", "index_err", "live_time"]:
-                if column_name in table.colnames:
-                    f[column_name] = table[column_name].data.flatten().tolist()
-        except KeyError:
-            self._logger.error(f"Incomplete data file; key not found in {data_config['file_name']}")
-            raise KeyError
+
+        if not time_min_max:
+            table["time_min"] = table["time"].data
+            table["time_max"] = table["time"].data + 0.1
+
+        # MJD filter
+        condition = np.ones(len(table), dtype=bool)
+        if mjd_min > -1:
+            condition &= table["time_min"] > mjd_min
+        if mjd_max > -1:
+            condition &= table["time_max"] < mjd_max
+        table = table[condition]
+
+        f = {}
+        f["time_min"] = table["time_min"].data.tolist()
+        f["time_max"] = table["time_max"].data.tolist()
+        f["flux"] = table["flux"].data.flatten().tolist()
+        if "flux_err" in table.colnames:
+            f["flux_err"] = table["flux_err"].data.flatten().tolist()
+        else:
+            up = table["flux_up"].data.flatten().tolist()
+            down = table["flux_down"].data.flatten().tolist()
+            f["flux_err"] = [0.5 * abs(u - d) for u, d in zip(up, down)]
+        f["MJD"] = [0.5 * (a + b) for a, b in zip(f["time_min"], f["time_max"])]
+        f["MJD_err"] = [0.5 * (b - a) for a, b in zip(f["time_min"], f["time_max"])]
+        if "flux_ul" in table.colnames:
+            flux_ul = table["flux_ul"].data.flatten().tolist()
+            is_ul = table["is_ul"].data.flatten().tolist()
+            f["flux_ul"] = [flux if is_ul else -1.0 for flux, is_ul in zip(flux_ul, is_ul)]
+        else:
+            f["flux_ul"] = [-1.0 if fe > 0 else fl for fe, fl in zip(f["flux_err"], f["flux"])]
+        for column_name in ["significance", "index", "index_err", "live_time"]:
+            if column_name in table.colnames:
+                f[column_name] = table[column_name].data.flatten().tolist()
+
         return f
diff --git a/v2dl5/orbital_phase.py b/v2dl5/orbital_phase.py
index 05b3ded..938cdc2 100644
--- a/v2dl5/orbital_phase.py
+++ b/v2dl5/orbital_phase.py
@@ -122,7 +122,7 @@ def get_orbital_phase_range(mjd_min, mjd_max, phase_mean, orbital_period, mjd_0,
     return ph_err
 
 
-def get_orbit_number(mjd, orbital_period, mjd_0):
+def get_orbit_number(mjd, orbital_period, mjd_orbit_count):
     """
     Return orbit number for a given MJD since MJD0.
 
@@ -132,15 +132,15 @@ def get_orbit_number(mjd, orbital_period, mjd_0):
         MJD
     orbital_period: float
         Orbital period (in units of days)
-    mjd_0: float
-        Reference MJD
+    mjd_orbit_count: float
+        Reference MJD for counting orbits
 
     Returns
     -------
     int
         Orbit number
     """
-    return math.ceil((mjd - mjd_0) / orbital_period)
+    return math.ceil((mjd - mjd_orbit_count) / orbital_period)
 
 
 # def getMJDOrbitZeroPhase(mjd, object="HESS J0632+057", orbital_period=317.3):
diff --git a/v2dl5/scripts/plot_binary_light_curves.py b/v2dl5/scripts/plot_binary_light_curves.py
index 8cecf9e..3872e19 100644
--- a/v2dl5/scripts/plot_binary_light_curves.py
+++ b/v2dl5/scripts/plot_binary_light_curves.py
@@ -123,6 +123,25 @@ def main():
             file_type=args.plot_type,
             figure_dir=args.figure_dir,
         )
+        plotter.plot_distribution(
+            instrument=args.instrument,
+            y_axis=y_key,
+            file_type=args.plot_type,
+            figure_dir=args.figure_dir,
+        )
+
+    plotter.plot_index_vs_flux(
+        instrument=args.instrument,
+        file_type=args.plot_type,
+        figure_dir=args.figure_dir,
+    )
+
+    plotter.plot_live_time_vs_phase_bin(
+        instrument=args.instrument,
+        phase_bins=args.orbital_bins,
+        file_type=args.plot_type,
+        figure_dir=args.figure_dir,
+    )
 
 
 if __name__ == "__main__":