From b7f1486f70446d98c8e72cd3ae986a9cb07137ef Mon Sep 17 00:00:00 2001
From: Metin San <metinisan@gmail.com>
Date: Mon, 16 Aug 2021 14:02:04 +0200
Subject: [PATCH] Set non observed pixels to NAN

---
 zodipy/simulation.py | 34 +++++++++++++++++++++++-----------
 zodipy/zodi.py       |  5 ++---
 2 files changed, 25 insertions(+), 14 deletions(-)

diff --git a/zodipy/simulation.py b/zodipy/simulation.py
index bcc00c9..159b061 100644
--- a/zodipy/simulation.py
+++ b/zodipy/simulation.py
@@ -3,6 +3,7 @@
 from itertools import repeat
 from math import radians
 from typing import Iterable, List
+import warnings
 
 import healpy as hp
 import numpy as np
@@ -52,9 +53,9 @@ def simulate(self, nside: int, freq: float, mask: float) -> np.ndarray:
         freq : float
             Frequency [GHz] at which to evaluate the IPD model.
         mask : float, optional
-            Angle [deg] between observer and the sun for which all pixels 
+            Angle [deg] between observer and the Sun for which all pixels 
             are masked at each observation. A mask of 90 degrees can be 
-            selected to simulate an observer that never looks inwards the sun.
+            selected to simulate an observer that never looks inwards the Sun.
             
         Returns
         -------
@@ -63,10 +64,10 @@ def simulate(self, nside: int, freq: float, mask: float) -> np.ndarray:
         """
 
     @staticmethod
-    def get_unmasked_pixels(
+    def get_observer_pixels(
         X_observer: np.ndarray, X_unit: np.ndarray, ang: float
     ) -> List[np.ndarray]:
-        """Returns the unmasked pixels.
+        """Returns a list of the observed pixels per observation.
         
         All pixels that have an angular distance of larger than some angle
         between the observer and the sun are masked.
@@ -78,7 +79,7 @@ def get_unmasked_pixels(
         X_unit: `np.ndarray`
             Array containing heliocentric unit vectors.
         ang: float
-            Angle for which all pixels are masked.
+            Angle for which all pixels are masked [deg].
         
         Returns
         -------
@@ -117,13 +118,17 @@ def simulate(self, nside: int, freq: float, mask: float) -> np.ndarray:
         X_earth  = self.earth_locations
         X_unit = np.asarray(hp.pix2vec(nside, pixels))
 
+        n_observations = len(X_observer)
+
         if mask is None:
-            pixels = list(repeat(pixels, n_observations := len(X_observer)))
+            pixels = list(repeat(pixels, n_observations))
         else:
             pixels = self.get_unmasked_pixels(X_observer, X_unit, ang=mask)
 
         components = self.model.components
-        emission = np.zeros((n_observations, len(components), npix))
+
+        # The emission is initialized as NANs representing unobserved pixels
+        emission = np.zeros((n_observations, len(components), npix)) + np.NAN
 
         for observation_idx, (observer_pos, earth_pos) in enumerate(zip(X_observer, X_earth)):
             observed_pixels = pixels[observation_idx]
@@ -131,13 +136,13 @@ def simulate(self, nside: int, freq: float, mask: float) -> np.ndarray:
 
             for comp_idx, (comp_name, comp) in enumerate(components.items()):
                 integration_config = self.integration_config[comp_name]
-                R, dR = integration_config.R, integration_config.dR
+                R = integration_config.R
 
                 comp_emission = comp.get_emission(
                     freq, observer_pos, earth_pos, unit_vectors, R
                 )
                 integrated_comp_emission = integration_config.integrator(
-                    comp_emission, R, dx=dR, axis=0
+                    comp_emission, R, dx=integration_config.dR, axis=0
                 )
 
                 comp_emissivity = self.model.emissivities.get_emissivity(
@@ -146,5 +151,12 @@ def simulate(self, nside: int, freq: float, mask: float) -> np.ndarray:
                 integrated_comp_emission *= comp_emissivity
 
                 emission[observation_idx, comp_idx, observed_pixels] = integrated_comp_emission
-        
-        return emission.mean(axis=0) * 1e20
\ No newline at end of file
+
+        with warnings.catch_warnings():
+            # np.nanmean throws a RuntimeWarning if all pixels along an 
+            # axis is NANs. This may occur when parts of the sky is left
+            # unobserver over all observations. Here we manually disable 
+            # the warning in the aforementioned scenario.
+            warnings.filterwarnings("ignore", category=RuntimeWarning)
+
+            return np.nanmean(emission, axis=0) * 1e20
\ No newline at end of file
diff --git a/zodipy/zodi.py b/zodipy/zodi.py
index 8eb46b1..ebee4c7 100644
--- a/zodipy/zodi.py
+++ b/zodipy/zodi.py
@@ -1,6 +1,5 @@
-from collections.abc import Iterable as Iterable_
-from typing import Optional, Union, Iterable
-from datetime import datetime, timedelta
+from typing import Optional, Union
+from datetime import datetime
 
 import astropy.units as u
 import numpy as np