Merge branch 'main' into adding-lunar-coords

.github/workflows/deploy.yaml

@@ -31,4 +31,4 @@ jobs:
           python -m build
 
       - name: Publish to PyPI
-        uses: pypa/[email protected].2
+        uses: pypa/[email protected].3

.github/workflows/run-docs-code.yaml

@@ -41,6 +41,17 @@ jobs:
           pip install papermill ipykernel
           conda list
 
+      - name: Install Poetry
+        run: |
+          pipx install poetry
+          poetry --version
+      - name: Install ska-ost-array
+        run: |
+          git clone https://gitlab.com/ska-telescope/ost/ska-ost-array-config.git
+          cd ska-ost-array-config
+          poetry install
+          cd ..
+
       - name: Install ipykernel
         run: python -m ipykernel install --user --name sense --display-name "sense"
 

.github/workflows/testsuite.yaml

@@ -40,6 +40,16 @@ jobs:
         run: |
           pip install .[test]
           conda list
+      - name: Install Poetry
+        run: |
+          pipx install poetry
+          poetry --version
+      - name: Install ska-ost-array
+        run: |
+          git clone https://gitlab.com/ska-telescope/ost/ska-ost-array-config.git
+          cd ska-ost-array-config
+          poetry install
+          cd ..
 
       - name: Run Tests
         run: |

.pre-commit-config.yaml

@@ -24,7 +24,7 @@ repos:
 
 - repo: https://github.com/astral-sh/ruff-pre-commit
   # Ruff version.
-  rev: v0.8.2
+  rev: v0.9.2
   hooks:
     # Run the linter.
     - id: ruff

.readthedocs.yml

@@ -17,6 +17,9 @@ build:
 conda:
   environment: docs/environment.yml
 
+sphinx:
+  configuration: docs/conf.py
+
 # Optionally set the version of Python and requirements required to build your docs
 python:
   install:

docs/tutorials/using_builtin_observatories.ipynb

@@ -146,6 +146,47 @@
     "print(obs.frequency)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can also create an observatory from a specific SKA configuration provided by `ska-ost-array-config` package (see https://gitlab.com/ska-telescope/ost/ska-ost-array-config for installation instructions) using the `from_ska` method:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "obs = Observatory.from_ska(\n",
+    "    subarray_type=\"AA*\", array_type=\"low\", Trcv=100.0 * un.K, frequency=75.0 * un.MHz\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "You can create a custom array by passing additional parameters to the array class (see `ska-ost-array-config` documentation). For example, to select all core stations, add 6 stations in the E1 cluster, and remove two stations from the core array:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "custom_obs = Observatory.from_ska(\n",
+    "    subarray_type=\"custom\",\n",
+    "    array_type=\"low\",\n",
+    "    Trcv=100.0 * un.K,\n",
+    "    frequency=75.0 * un.MHz,\n",
+    "    custom_stations=\"C*, E1-*\",\n",
+    "    exclude_stations=\"C1,C2\",\n",
+    ")"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},

src/py21cmsense/_utils.py

@@ -41,7 +41,12 @@ def find_nearest(array, value):
 
 @un.quantity_input
 def phase_past_zenith(
-        time_past_zenith: un.day, bls_enu: np.ndarray, latitude, world: str = "earth", use_apparent: bool = True
+    time_past_zenith: un.hour,
+    bls_enu: np.ndarray,
+    latitude: float,
+    world: str = "earth",
+    phase_center_dec: un.rad = None,
+    use_apparent: bool = True,
 ):
     """Compute UVWs phased to a point rotated from zenith by a certain amount of time.
 
@@ -52,14 +57,21 @@ def phase_past_zenith(
     Parameters
     ----------
     time_past_zenith
-        The time passed since the point was at zenith. If float, assumed to be in units
-        of days.
-    uvws0 : array
-        The UVWs when phased to zenith.
+        The time passed since the point was at its closest to zenith. Must be a
+        quantity with time units.
+    bls_enu
+        An (Nbls, 3)-shaped array containing the ENU coordinates of the baseline
+        vectors (equivalent to the UVWs if phased to zenith).
     latitude
         The latitude of the center of the array, in radians.
     world
         Whether the telescope is on the Earth or Moon.
+    phase_center_dec
+        If given, the declination of the phase center. If not given, it is set to the
+        latitude of the array (i.e. the phase center passes through zenith).
+    use_apparent
+        Whether to use the apparent coordinates of the phase center (i.e. after
+        accounting for nutation and precession etc.)
 
     Returns
     -------
@@ -77,35 +89,49 @@ def phase_past_zenith(
     jd = 2454600
 
     if world == "earth":
-        zenith_coord = SkyCoord(
+        tm = Time(jd, format="jd")
+
+        phase_center_coord = SkyCoord(
             alt=90 * un.deg,
             az=0 * un.deg,
-            obstime=Time(jd, format="jd"),
+            obstime=tm,
             frame="altaz",
             location=telescope_location,
         )
     else:
-        zenith_coord = LunarSkyCoord(
+        tm = LTime(jd, format='jd')
+        phase_center_coord = LunarSkyCoord(
                 alt=90 * un.deg,
                 az=0 * un.deg,
-                obstime=LTime(jd, format="jd"),
+                obstime=tm,
                 frame="lunartopo",
                 location=telescope_location,
             )
 
-    zenith_coord = zenith_coord.transform_to("icrs")
+    phase_center_coord = phase_center_coord.transform_to("icrs")
+
+    if phase_center_dec is not None:
+        phase_center_coord = phase_center_coord.__class__(
+            ra=phase_center_coord.ra,
+            dec=phase_center_dec,
+            obstime=tm,
+            frame=phase_center_coord.frame,
+            location=telescope_location,
+        )
+
+    phase_center_coord.obstime.location = telescope_location
+
 
-    zenith_coord.obstime.location = telescope_location
-    obstimes = zenith_coord.obstime + time_past_zenith
+    obstimes = phase_center_coord.obstime + time_past_zenith
     lsts = obstimes.sidereal_time("apparent", longitude=0.0).rad
 
     if not hasattr(lsts, "__len__"):
         lsts = np.array([lsts])
 
     if use_apparent:
         app_ra, app_dec = uvutils.phasing.calc_app_coords(
-            lon_coord=zenith_coord.ra.to_value("rad"),
-            lat_coord=zenith_coord.dec.to_value("rad"),
+            lon_coord=phase_center_coord.ra.to_value("rad"),
+            lat_coord=phase_center_coord.dec.to_value("rad"),
             time_array=obstimes.utc.jd,
             telescope_loc=telescope_location,
         )
@@ -114,8 +140,8 @@ def phase_past_zenith(
         app_dec = np.tile(app_dec, len(bls_enu))
 
     else:
-        app_ra = zenith_coord.ra.to_value("rad") * np.ones(len(bls_enu) * len(lsts))
-        app_dec = zenith_coord.dec.to_value("rad") * np.ones(len(bls_enu) * len(lsts))
+        app_ra = phase_center_coord.ra.to_value("rad") * np.ones(len(bls_enu) * len(lsts))
+        app_dec = phase_center_coord.dec.to_value("rad") * np.ones(len(bls_enu) * len(lsts))
 
     # Now make everything nbls * ntimes big.
     _lsts = np.tile(lsts, len(bls_enu))

src/py21cmsense/baseline_filters.py

@@ -58,7 +58,7 @@ class BaselineRange(BaselineFilter):
     def _bl_max_vld(self, att, val):
         if val <= self.bl_min:
             raise ValueError(
-                "bl_max must be greater than bl_min, got " f"bl_min={self.bl_min} and bl_max={val}"
+                f"bl_max must be greater than bl_min, got bl_min={self.bl_min} and bl_max={val}"
             )
 
     def __call__(self, bl: tp.Length) -> bool:

src/py21cmsense/observation.py

@@ -88,6 +88,11 @@ class Observation:
         that use astropy are preferred.
     cosmo : LambdaCDM
         An astropy cosmology object to use.
+    phase_center_dec
+        The declination of the phase center (typically used for a tracked scan
+        observation). By default, this is set to the latitude of the observatory
+        (i.e. the phase-center passes through zenith, like a typical drift-scan
+        observation).
     """
 
     observatory: obs.Observatory = attr.ib(validator=vld.instance_of(obs.Observatory))
@@ -128,6 +133,7 @@ class Observation:
     tsky_ref_freq: tp.Frequency = attr.ib(default=150 * un.MHz, validator=ut.positive)
     use_approximate_cosmo: bool = attr.ib(default=False, converter=bool)
     cosmo: LambdaCDM = attr.ib(default=Planck15, converter=Planck15.from_format)
+    phase_center_dec = attr.ib(validator=(tp.vld_physical_type("angle")))
 
     @classmethod
     def from_yaml(cls, yaml_file):
@@ -210,6 +216,10 @@ def _lst_bin_size_default(self):
         else:
             return self.observatory.observation_duration
 
+    @phase_center_dec.default
+    def _phase_center_dec_default(self):
+        return self.observatory.latitude
+
     @cached_property
     def beam_crossing_time(self):
         """The time it takes for a source to cross the beam.
@@ -292,6 +302,7 @@ def uv_coverage(self) -> np.ndarray:
             integration_time=self.integration_time,
             observation_duration=self.lst_bin_size,
             ndecimals=self.redundancy_tol,
+            phase_center_dec=self.phase_center_dec,
         )
 
     @cached_property