From ddcba2517b463c0b40d5389062b757c57012a116 Mon Sep 17 00:00:00 2001
From: Alan Loh <alan.loh.up7@gmail.com>
Date: Wed, 24 Jan 2024 17:40:40 +0100
Subject: [PATCH] Source contamination - set the array factor to NaN if the
 pointing is below the horizon

---
 nenupy/__init__.py               |  2 +-
 nenupy/schedule/contamination.py | 20 +++++++++++++++-----
 2 files changed, 16 insertions(+), 6 deletions(-)

diff --git a/nenupy/__init__.py b/nenupy/__init__.py
index 63983d5..744f7d7 100644
--- a/nenupy/__init__.py
+++ b/nenupy/__init__.py
@@ -5,7 +5,7 @@
 __copyright__ = "Copyright 2023, nenupy"
 __credits__ = ["Alan Loh"]
 __license__ = "MIT"
-__version__ = "2.6.2"
+__version__ = "2.6.3"
 __maintainer__ = "Alan Loh"
 __email__ = "alan.loh@obspm.fr"
 
diff --git a/nenupy/schedule/contamination.py b/nenupy/schedule/contamination.py
index 24d1ed7..bfb02e3 100644
--- a/nenupy/schedule/contamination.py
+++ b/nenupy/schedule/contamination.py
@@ -84,7 +84,7 @@ def lst_time(self):
 
     # --------------------------------------------------------- #
     # ------------------------ Methods ------------------------ #
-    def plot(self, time_unit : str = 'utc', **kwargs):
+    def plot(self, time_unit: str = "utc", **kwargs):
         """
         :param time_unit:
             To select in which units the time axe is to be displayed : 'utc' or 'lst'. Default is 'utc'. 
@@ -95,13 +95,15 @@ def plot(self, time_unit : str = 'utc', **kwargs):
         """
         fig = plt.figure(figsize=kwargs.get("figsize", (10, 6)))
 
-        if time_unit == 'utc' :
+        if time_unit == "utc":
             time_to_plot = self.time.datetime
             time_label = f"Time (UTC since {self.time[0].isot.split('.')[0]})"
-        elif time_unit == 'lst' :
+        elif time_unit == "lst":
             time_to_plot = self.lst_time.rad
             time_label = "Local Sidereal Time (rad)"
-        
+
+
+        # TODO : roll if there is a gap in lst and or utc
         plt.pcolormesh(
             time_to_plot,
             self.frequency.to(u.MHz).value,
@@ -118,7 +120,7 @@ def plot(self, time_unit : str = 'utc', **kwargs):
 
         # Format of the time axis tick labels
         ax = plt.gca()
-        if time_unit == 'utc' :
+        if time_unit == "utc":
             ax.xaxis.set_major_formatter(
                 mdates.DateFormatter("%H:%M:%S")
             )
@@ -506,6 +508,8 @@ def _compute_array_factor(self, use_antenna_gain: bool = True) -> np.ndarray:
         # Mini-Array selection
         ma_mask = np.isin(MA_ROTATIONS, self.miniarray_rotations)
         af = 0
+        below_horizon_mask = self.pointing.horizontal_coordinates.alt.deg <= 0
+
         # Compute the array factor for each Mini-Array
         for m in MA_INDICES[ma_mask]:
             ma = MiniArray(index=m)
@@ -523,6 +527,12 @@ def _compute_array_factor(self, use_antenna_gain: bool = True) -> np.ndarray:
         log.info("Computing the array factor...")
         with ProgressBar() if log.getEffectiveLevel() <= logging.INFO else DummyCtMgr():
             array_factor = af.compute()
+
+        # Set the array factor to NaN if the pointing is below the horizon
+        if np.any(below_horizon_mask):
+            array_factor[below_horizon_mask, ...] = np.nan
+            log.warning("Some time samples match a sub-horizon pointing, the corresponding array-factor is set to NaN.")
+
         return array_factor / np.max(array_factor, axis=3)[:, :, :, None]