diff --git a/.github/workflows/deploy.yaml b/.github/workflows/deploy.yaml index ada9892..21c46d7 100644 --- a/.github/workflows/deploy.yaml +++ b/.github/workflows/deploy.yaml @@ -31,4 +31,4 @@ jobs: python -m build - name: Publish to PyPI - uses: pypa/gh-action-pypi-publish@v1.12.2 + uses: pypa/gh-action-pypi-publish@v1.12.3 diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 6cff94d..37d7b37 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -24,7 +24,7 @@ repos: - repo: https://github.com/astral-sh/ruff-pre-commit # Ruff version. - rev: v0.8.3 + rev: v0.8.4 hooks: # Run the linter. - id: ruff diff --git a/src/py21cmsense/_utils.py b/src/py21cmsense/_utils.py index 1b4ad4b..b8badc4 100644 --- a/src/py21cmsense/_utils.py +++ b/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,17 +79,27 @@ 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__"): @@ -86,8 +107,8 @@ def phase_past_zenith( 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)) diff --git a/src/py21cmsense/observation.py b/src/py21cmsense/observation.py index f012873..39d3117 100644 --- a/src/py21cmsense/observation.py +++ b/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 diff --git a/src/py21cmsense/observatory.py b/src/py21cmsense/observatory.py index 6ab6521..fe7d084 100644 --- a/src/py21cmsense/observatory.py +++ b/src/py21cmsense/observatory.py @@ -284,7 +284,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). @@ -300,6 +303,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 ------- @@ -313,7 +319,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: @@ -472,6 +483,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. @@ -538,9 +550,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)) diff --git a/src/py21cmsense/units.py b/src/py21cmsense/units.py index 279136b..1276dae 100644 --- a/src/py21cmsense/units.py +++ b/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 = [ ( diff --git a/tests/test_observation.py b/tests/test_observation.py index 3e163ca..abcac47 100644 --- a/tests/test_observation.py +++ b/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)