Update visibility determination

Fixes #28

gravitational_wave_toy/observe.py

@@ -4,7 +4,6 @@
 
 import astropy.units as u
 
-# %matplotlib inline
 import matplotlib.pyplot as plt
 import numpy as np
 from astropy.coordinates import Distance
@@ -503,6 +502,7 @@ def observe(
         max_energy: u.Quantity = None,
         max_time: u.Quantity = 12 * u.hour,
         target_precision: u.Quantity = 1 * u.s,
+        n_time_steps: int = 100,
         max_iter=100,
         sensitivity_mode: Literal["sensitivity", "photon_flux"] = "sensitivity",
     ):
@@ -545,26 +545,73 @@ def observe(
             # start the procedure
             self.start_time = start_time
             delay = start_time
-
-            # check maximum time
-            visible = self.check_if_visible(delay, max_time + delay, sensitivity, sensitivity_mode=sensitivity_mode)
-
-            # not visible even after maximum observation time
-            if not visible:
+
+            ''' four cases:
+            1. average fluence is always above the sensitivity curve -> immediately detectable
+            2. average fluence is monotonically increasing -> intersection point is moment of detectability
+            3. average fluence is always below the sensitivity curve -> never detectable
+            4. average fluence crosses sensitivity curve multiple times -> first point where it is above sensitivity is the moment of detectability
+            '''
+
+            # Check immediate detectability
+            visible = self.check_if_visible(delay, target_precision + delay, sensitivity, sensitivity_mode=sensitivity_mode)
+
+            if visible:
+                self.end_time = target_precision + delay
+                self.obs_time = target_precision
+                self.seen = True
+                # print("Immediate detectability")
                 return self.output()
-
+
+            # Perform an incremental search
+            time_steps = np.logspace(
+                np.log10(target_precision.value),
+                np.log10(max_time.value),
+                n_time_steps,
+            ) * u.s
+
+            detected_in_scan = False
+
+            for i, time_step in enumerate(time_steps):
+
+                previous_time = start_time + time_step
+
+                if i == (n_time_steps - 1):
+                    break
+                else:
+                    current_time = start_time + time_steps[i + 1]
+
+                visible = self.check_if_visible(start_time, previous_time, sensitivity, sensitivity_mode=sensitivity_mode)
+                # print(f"Start time: {start_time}, Current time: {current_time}, Previous time: {previous_time}, {visible}")
+
+                if visible:
+                    # print(f"Detected between {previous_time} and {current_time}")
+                    detected_in_scan = True
+                    break
+
+            if not detected_in_scan:
+                self.end_time = -1 * u.s
+                self.obs_time = -1 * u.s
+                self.seen = False
+                # print("Never detectable")
+                return self.output()
+
+            # Fine-tune search with bisection
+            # print(f"Fine-tuning search with bisection between {start_time} and {current_time}")
             end_time, seen = self._bisect_find_zeros(
-                sensitivity, start_time, max_time + start_time, target_precision, max_iter=max_iter
+                sensitivity, start_time, current_time, target_precision, max_iter=max_iter, lowest_possible=previous_time, highest_possible=current_time
             )
 
             if seen:
                 self.end_time = end_time
                 self.obs_time = end_time - start_time
                 self.seen = True
+                # print(f"Detected at {end_time}")
             else:
                 self.end_time = -1 * u.s
                 self.obs_time = -1 * u.s
                 self.seen = False
+                # print("Never detectable, but weirdly")
 
             # just return dict now
             return self.output()