Fix phase-wrapped binning and add negative phase support

eclipsebin/binning.py

@@ -75,6 +75,22 @@ def __init__(
             "atol_secondary": None,
         }
 
+        # Detect and store original phase range for denormalization later
+        self._original_phase_min = np.min(phases)
+        self._original_phase_max = np.max(phases)
+        self._original_phase_range = self._original_phase_max - self._original_phase_min
+
+        # Normalize phases to [0, 1] if not already
+        if self._original_phase_min < 0 or self._original_phase_max > 1:
+            self._needs_denormalization = True
+            normalized_phases = self._normalize_phases(phases)
+            sort_idx = np.argsort(normalized_phases)
+            self.data["phases"] = normalized_phases[sort_idx]
+            self.data["fluxes"] = fluxes[sort_idx]
+            self.data["flux_errors"] = flux_errors[sort_idx]
+        else:
+            self._needs_denormalization = False
+
         self.set_atol(primary=atol_primary, secondary=atol_secondary)
 
         # Identify primary and secondary eclipse minima (in original phase space)
@@ -95,6 +111,34 @@ def __init__(
             self.primary_eclipse = self.get_eclipse_boundaries(primary=True)
             self.secondary_eclipse = self.get_eclipse_boundaries(primary=False)
 
+    def _normalize_phases(self, phases):
+        """
+        Normalize phases from original range to [0, 1].
+
+        Args:
+            phases (np.ndarray): Phases in original range
+
+        Returns:
+            np.ndarray: Phases normalized to [0, 1]
+        """
+        # Shift so minimum is at 0, then scale to [0, 1]
+        return (phases - self._original_phase_min) / self._original_phase_range
+
+    def _denormalize_phases(self, phases):
+        """
+        Convert phases from [0, 1] back to original range.
+
+        Args:
+            phases (np.ndarray): Phases in [0, 1] range
+
+        Returns:
+            np.ndarray: Phases in original range
+        """
+        if not self._needs_denormalization:
+            return phases
+        # Scale from [0, 1] back to original range
+        return phases * self._original_phase_range + self._original_phase_min
+
     def find_minimum_flux_phase(self):
         """
         Finds the phase of the minimum flux, corresponding to the primary eclipse.
@@ -426,17 +470,20 @@ def find_bin_edges(self):
 
     def _rewrap_to_original_phase(self, phases_array):
         """
-        Rewrap phases back to original phase space before unwrapping.
+        Rewrap phases back to original phase space before unwrapping,
+        then denormalize if original input had non-standard range.
 
         Args:
-            phases_array (np.ndarray): Array of phases in unwrapped space
+            phases_array (np.ndarray): Array of phases in unwrapped [0, 1] space
 
         Returns:
             np.ndarray: Phases shifted back to original space
         """
-        if self._phase_shift == 0.0:
-            return phases_array
-        return (phases_array - self._phase_shift) % 1.0
+        result = phases_array
+        if self._phase_shift != 0.0:
+            result = (result - self._phase_shift) % 1.0
+        # Denormalize back to original range (e.g., [-0.5, 0.5])
+        return self._denormalize_phases(result)
 
     def calculate_bins(self, return_in_original_phase=True):
         """
@@ -620,7 +667,10 @@ def plot_binned_light_curve(self, bin_centers, bin_means, bin_stds):
         )
         plt.xlabel("Phases", fontsize=14)
         plt.ylabel("Normalized Flux", fontsize=14)
-        plt.xlim(0, 1)
+        if self._needs_denormalization:
+            plt.xlim(self._original_phase_min, self._original_phase_max)
+        else:
+            plt.xlim(0, 1)
         ylims = plt.ylim()
 
         # Get eclipse boundaries in original phase space
@@ -691,7 +741,10 @@ def plot_unbinned_light_curve(self):
             label="Secondary Eclipse",
         )
         plt.ylim(ylims)
-        plt.xlim(0, 1)
+        if self._needs_denormalization:
+            plt.xlim(self._original_phase_min, self._original_phase_max)
+        else:
+            plt.xlim(0, 1)
         plt.ylabel("Normalized Flux", fontsize=14)
         plt.xlabel("Phases", fontsize=14)
         plt.legend()

setup.py

@@ -6,7 +6,7 @@
 
 setup(
     name="eclipsebin",
-    version="0.3.3",
+    version="0.4.0",
     author="Jackie Blaum",
     author_email="jackie.blaum@gmail.com",
     description="A specialized binning scheme for eclipsing binary star light curves",

tests/test_eclipsing_binary_binner.py

@@ -527,3 +527,107 @@ def test_both_eclipses_near_boundary(both_near_boundary_light_curve):
     bin_centers, bin_means, bin_errors = binner.bin_light_curve(plot=False)
     assert abs(len(bin_centers) - 100) <= 2  # Allow ±2 bins due to duplicates='drop'
     assert np.all(bin_errors > 0)
+
+
+@pytest.fixture
+def negative_phase_light_curve():
+    """
+    Fixture for light curve with phases in [-0.5, 0.5] range (PHOEBE-style).
+    Primary eclipse at phase 0, secondary eclipse wrapped around ±0.5 boundary.
+    """
+    np.random.seed(99)
+    # Phases from -0.5 to 0.5 (PHOEBE convention)
+    phases = np.linspace(-0.5, 0.4999, 10000)
+    fluxes = np.ones_like(phases)
+
+    # Primary eclipse centered at phase 0
+    # Indices 4500-5500 correspond to phases around 0
+    fluxes[4500:5000] = np.linspace(0.95, 0.8, 500)
+    fluxes[5000:5500] = np.linspace(0.8, 0.95, 500)
+
+    # Secondary eclipse wrapping around ±0.5 boundary
+    # Near phase 0.5 (end of array) and -0.5 (start of array)
+    fluxes[9700:10000] = np.linspace(0.95, 0.9, 300)  # phase ~0.47 to 0.5
+    fluxes[0:300] = np.linspace(0.9, 0.95, 300)  # phase -0.5 to ~-0.47
+
+    flux_errors = np.random.normal(0.01, 0.001, 10000)
+
+    # Random subset
+    random_indices = np.random.choice(range(len(phases)), size=5000, replace=False)
+    return phases[random_indices], fluxes[random_indices], flux_errors[random_indices]
+
+
+def test_negative_phase_input(negative_phase_light_curve):
+    """Test that negative phase inputs (PHOEBE-style [-0.5, 0.5]) are handled correctly."""
+    phases, fluxes, flux_errors = negative_phase_light_curve
+
+    # Verify input has negative phases
+    assert np.min(phases) < 0, "Test fixture should have negative phases"
+
+    # Should not raise an error
+    binner = EclipsingBinaryBinner(
+        phases, fluxes, flux_errors, nbins=100, fraction_in_eclipse=0.2
+    )
+
+    # Binning should succeed
+    bin_centers, bin_means, bin_errors = binner.bin_light_curve(plot=False)
+
+    # Results should be in original phase space [-0.5, 0.5]
+    assert np.min(bin_centers) >= -0.5, "Bin centers should be >= -0.5"
+    assert np.max(bin_centers) <= 0.5, "Bin centers should be <= 0.5"
+
+    # Should have expected number of bins (allow small tolerance)
+    assert abs(len(bin_centers) - 100) <= 2
+
+    # All bin errors should be positive
+    assert np.all(bin_errors > 0)
+
+
+def test_negative_phase_primary_at_zero(negative_phase_light_curve):
+    """Test that primary eclipse detection works when primary is at phase 0."""
+    phases, fluxes, flux_errors = negative_phase_light_curve
+
+    binner = EclipsingBinaryBinner(
+        phases, fluxes, flux_errors, nbins=100, fraction_in_eclipse=0.2
+    )
+
+    # Primary minimum should be near phase 0 in original space
+    # In normalized space it's at 0.5, but we want to verify detection worked
+    primary_min_original = binner._denormalize_phases(
+        np.array([binner.primary_eclipse_min_phase])
+    )[0]
+    assert (
+        -0.1 < primary_min_original < 0.1
+    ), f"Primary eclipse should be near phase 0, got {primary_min_original}"
+
+
+def test_negative_phase_secondary_wrapped():
+    """Test secondary eclipse that wraps around ±0.5 boundary in PHOEBE-style phases."""
+    np.random.seed(101)
+    phases = np.linspace(-0.5, 0.4999, 10000)
+    fluxes = np.ones_like(phases)
+
+    # Primary at phase 0
+    fluxes[4800:5200] = np.linspace(0.95, 0.8, 400)
+
+    # Secondary wrapping around ±0.5 (the boundary in PHOEBE space)
+    fluxes[9800:10000] = np.linspace(0.97, 0.92, 200)  # near +0.5
+    fluxes[0:200] = np.linspace(0.92, 0.97, 200)  # near -0.5
+
+    flux_errors = np.abs(np.random.normal(0.01, 0.001, 10000))
+
+    random_indices = np.random.choice(range(len(phases)), size=5000, replace=False)
+    phases = phases[random_indices]
+    fluxes = fluxes[random_indices]
+    flux_errors = flux_errors[random_indices]
+
+    binner = EclipsingBinaryBinner(
+        phases, fluxes, flux_errors, nbins=100, fraction_in_eclipse=0.2
+    )
+
+    bin_centers, bin_means, bin_errors = binner.bin_light_curve(plot=False)
+
+    # Verify output is in original range
+    assert np.min(bin_centers) >= -0.5
+    assert np.max(bin_centers) <= 0.5
+    assert np.all(bin_errors > 0)