feat: add ability to set declination of phase center

src/py21cmsense/_utils.py

@@ -39,7 +39,11 @@ def find_nearest(array, value):
 
 @un.quantity_input
 def phase_past_zenith(
-    time_past_zenith: un.day, bls_enu: np.ndarray, latitude, use_apparent: bool = True
+    time_past_zenith: un.hour,
+    bls_enu: np.ndarray,
+    latitude: float,
+    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.
 
@@ -50,12 +54,19 @@ 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.
+    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
     -------
@@ -68,26 +79,36 @@ def phase_past_zenith(
 
     # JD is arbitrary
     jd = 2454600
+    tm = Time(jd, format="jd")
 
-    zenith_coord = SkyCoord(
+    phase_center_coord = SkyCoord(
         alt=90 * un.deg,
         az=0 * un.deg,
-        obstime=Time(jd, format="jd"),
+        obstime=tm,
         frame="altaz",
         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 = SkyCoord(
+            ra=phase_center_coord.ra,
+            dec=phase_center_dec,
+            obstime=tm,
+            frame="icrs",
+            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,
         )
@@ -96,8 +117,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/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))
@@ -129,6 +134,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):
@@ -182,6 +188,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.
@@ -264,6 +274,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

src/py21cmsense/observatory.py

@@ -233,7 +233,10 @@ def baselines_metres(self) -> tp.Meters:
         return (self.antpos[np.newaxis, :, :] - self.antpos[:, np.newaxis, :]).to(un.m)
 
     def projected_baselines(
-        self, baselines: tp.Length | None = None, time_offset: tp.Time = 0 * un.hour
+        self,
+        baselines: tp.Length | None = None,
+        time_offset: tp.Time = 0 * un.hour,
+        phase_center_dec: tp.Angle | None = None,
     ) -> np.ndarray:
         """Compute the *projected* baseline lengths (in wavelengths).
 
@@ -249,6 +252,9 @@ def projected_baselines(
         time_offset
             The amount of time elapsed since the phase center was at zenith.
             Assumed to be in days unless otherwise defined. May be negative.
+        phase_center_dec
+            The declination of the phase center of the observation. By default, the
+            same as the latitude of the array.
 
         Returns
         -------
@@ -262,7 +268,12 @@ def projected_baselines(
 
         bl_wavelengths = baselines.reshape((-1, 3)) * self.metres_to_wavelengths
 
-        out = ut.phase_past_zenith(time_offset, bl_wavelengths, self.latitude)
+        out = ut.phase_past_zenith(
+            time_past_zenith=time_offset,
+            bls_enu=bl_wavelengths,
+            latitude=self.latitude,
+            phase_center_dec=phase_center_dec,
+        )
 
         out = out.reshape(*orig_shape[:-1], np.size(time_offset), orig_shape[-1])
         if np.size(time_offset) == 1:
@@ -421,6 +432,7 @@ def grid_baselines(
         baseline_filters: Callable | tuple[Callable] = (),
         observation_duration: tp.Time | None = None,
         ndecimals: int = 1,
+        phase_center_dec: tp.Angle | None = None,
     ) -> np.ndarray:
         """
         Grid baselines onto a pre-determined uvgrid, accounting for earth rotation.
@@ -487,9 +499,9 @@ def grid_baselines(
 
         time_offsets = self.time_offsets_from_obs_int_time(integration_time, observation_duration)
 
-        uvws = self.projected_baselines(baselines, time_offsets).reshape(
-            baselines.shape[0], time_offsets.size, 3
-        )
+        uvws = self.projected_baselines(
+            baselines, time_offsets, phase_center_dec=phase_center_dec
+        ).reshape(baselines.shape[0], time_offsets.size, 3)
 
         # grid each baseline type into uv plane
         dim = len(self.ugrid(bl_max))

src/py21cmsense/units.py

@@ -24,6 +24,7 @@ class UnitError(ValueError):
 TempSquared = un.Quantity[un.get_physical_type("temperature") ** 2]
 Wavenumber = un.Quantity[littleh / un.Mpc]
 Delta = un.Quantity[un.mK**2]
+Angle = un.Quantity["angle"]
 
 time_as_distance = [
     (

tests/test_observation.py

@@ -19,7 +19,9 @@ def bm():
 @pytest.fixture(scope="module")
 def observatory(bm):
     return Observatory(
-        antpos=np.array([[0, 0, 0], [14, 0, 0], [28, 0, 0], [70, 0, 0]]) * units.m,
+        antpos=np.array([[0, 0, 0], [14, 0, 0], [28, 0, 0], [70, 0, 0], [0, 14, 0], [23, -45, 0]])
+        * units.m,
+        latitude=-32 * units.deg,
         beam=bm,
     )
 
@@ -111,3 +113,41 @@ def test_trcv_func(observatory: Observatory):
     )
     assert obs.Trcv.unit == units.K
     assert obs.Trcv == (observatory.beam.frequency / units.MHz) * 0.01 * units.K
+
+
+def test_non_zenith_pointing(bm):
+    """Test that not pointing at zenith gives different results."""
+    observatory = Observatory(
+        antpos=np.array([[0, 0, 0], [14, 0, 0], [28, 0, 0], [70, 0, 0], [0, 14, 0], [23, -45, 0]])
+        * units.m,
+        latitude=-32 * units.deg,
+        beam=bm,
+    )
+    at_zenith = Observation(observatory=observatory)
+
+    almost_zenith = Observation(
+        observatory=observatory, phase_center_dec=observatory.latitude * 0.99
+    )
+    np.testing.assert_allclose(at_zenith.uv_coverage, almost_zenith.uv_coverage)
+
+    not_zenith = Observation(
+        observatory=observatory,
+        phase_center_dec=observatory.latitude + 45 * units.deg,
+    )
+    assert not np.allclose(not_zenith.uv_coverage, at_zenith.uv_coverage)
+
+
+def test_non_zenith_pointing_only_ew(bm):
+    """Test that not pointing at zenith gives the SAME results if all baselines are EW."""
+    observatory = Observatory(
+        antpos=np.array([[0, 0, 0], [14, 0, 0], [28, 0, 0], [70, 0, 0]]) * units.m,
+        latitude=-32 * units.deg,
+        beam=bm,
+    )
+    at_zenith = Observation(observatory=observatory)
+
+    not_zenith = Observation(
+        observatory=observatory,
+        phase_center_dec=observatory.latitude + 45 * units.deg,
+    )
+    np.testing.assert_allclose(not_zenith.uv_coverage, at_zenith.uv_coverage)