Improve performance by avoiding copying the skycoord object and inste…

…ad using indexing to filter it

README.md

@@ -26,23 +26,26 @@ See the [documentation](https://cosmoglobe.github.io/zodipy/) for more informati
 ## A simple example
 ```python
 import astropy.units as u
+from astropy.coordinates import SkyCoord
 from astropy.time import Time
 
-from zodipy import Zodipy
+from zodipy import Model
 
+model = Model(25*u.micron)
 
-model = Zodipy(model="dirbe")
-
-emission = model.get_emission_ang(
-    25 * u.micron,
-    theta=[10, 10.1, 10.2] * u.deg,
-    phi=[90, 89, 88] * u.deg,
-    obs_time=Time("2022-01-01 12:00:00"),
-    obs="earth",
+lon = [10, 10.1, 10.2] * u.deg
+lat = [90, 89, 88] * u.deg
+skycoord = SkyCoord(
+    lon,
+    lat,
+    obstime=Time("2022-01-01 12:00:00"),
+    frame="galactic",
 )
 
+emission = model.evaluate(skycoord, obspos="earth")
+
 print(emission)
-#> [15.35392831 15.35495051 15.35616009] MJy / sr
+#> [27.52410841 27.66581351 27.81270207] MJy / sr
 ```
 
 ## Related scientific papers

docs/index.md

@@ -21,27 +21,29 @@ ZodiPy is an [Astropy affiliated](https://www.astropy.org/affiliated/) package f
 ## A simple example
 ```python
 import astropy.units as u
+from astropy.coordinates import SkyCoord
 from astropy.time import Time
 
-from zodipy import Zodipy
+from zodipy import Model
 
+model = Model(25*u.micron)
 
-model = Zodipy("dirbe")
-
-emission = model.get_emission_ang(
-    25 * u.micron,
-    theta=[10, 10.1, 10.2] * u.deg,
-    phi=[90, 89, 88] * u.deg,
-    obs_time=Time("2022-01-01 12:00:00"),
-    obs="earth",
-    lonlat=True,
+lon = [10, 10.1, 10.2] * u.deg
+lat = [90, 89, 88] * u.deg
+skycoord = SkyCoord(
+    lon,
+    lat,
+    obstime=Time("2022-01-01 12:00:00"),
+    frame="galactic",
 )
 
+emission = model.evaluate(skycoord, obspos="earth")
+
 print(emission)
-#> [15.53095493 15.52883577 15.53121942] MJy / sr
+#> [27.52410841 27.66581351 27.81270207] MJy / sr
 ```
 
-What's going on here:
+What's going on here: THIS NEEDS TO BE UPDATED
 
 - We start by initializing the [`Zodipy`][zodipy.zodipy.Zodipy] class, which is our interface, where we specify that we want to use the DIRBE interplanetary dust model.
 - We use the [`get_emission_ang`][zodipy.zodipy.Zodipy.get_emission_ang] method which is a method to simulate emission from angular sky coordinates (see the [reference](reference.md) for other available simulation methods).

zodipy/zodipy.py

@@ -123,23 +123,19 @@ def evaluate(
             emission: Simulated zodiacal light in units of 'MJy/sr'.
 
         """
-        (unique_lon, unique_lat), indicies = np.unique(
-            np.vstack([skycoord.spherical.lon.value, skycoord.spherical.lat.value]),
-            return_inverse=True,
-            axis=1,
-        )
-
         if skycoord.obstime is None:
             msg = "The `obstime` attribute of the `SkyCoord` object must be set."
             raise ValueError(msg)
 
-        skycoord = coords.SkyCoord(
-            unique_lon,
-            unique_lat,
-            unit=units.deg,
-            frame=skycoord.frame,
-            obstime=skycoord.obstime,
+        # Pick out unique coordinates, and only calculate the emission for these. and return the
+        # inverse indices to map the output back to the original coordinates.
+        _, index, inverse = np.unique(
+            [skycoord.spherical.lon, skycoord.spherical.lat],
+            return_index=True,
+            return_inverse=True,
+            axis=1,
         )
+        skycoord = skycoord[index]
 
         earth_skycoord = get_earth_skycoord(skycoord.obstime, ephemeris=self._ephemeris)
         obs_skycoord = get_obs_skycoord(
@@ -216,7 +212,7 @@ def evaluate(
 
         else:
             integrated_comp_emission = np.zeros(
-                (self._interplanetary_dust_model.ncomps, skycoord.size)
+                (self._interplanetary_dust_model.ncomps, inverse.size)
             )
 
             for idx, comp_label in enumerate(self._interplanetary_dust_model.comps.keys()):
@@ -235,8 +231,8 @@ def evaluate(
                     * (stop[comp_label] - start[comp_label])
                 )
 
-        emission = np.zeros((self._interplanetary_dust_model.ncomps, indicies.size))
-        emission = integrated_comp_emission[:, indicies] << (units.MJy / units.sr)
+        emission = np.zeros((self._interplanetary_dust_model.ncomps, inverse.size))
+        emission = integrated_comp_emission[:, inverse] << (units.MJy / units.sr)
 
         return emission if return_comps else emission.sum(axis=0)