From 9c21c3f0ed4cbd4893b039ec338d684e9ec7e0e1 Mon Sep 17 00:00:00 2001 From: Steven Murray Date: Tue, 29 Oct 2024 10:18:43 +0100 Subject: [PATCH 1/6] feat: add ability to set declination of phase center --- src/py21cmsense/_utils.py | 37 +++++++++++++++++++++++++----- src/py21cmsense/observation.py | 11 +++++++++ src/py21cmsense/observatory.py | 22 ++++++++++++++---- src/py21cmsense/units.py | 1 + tests/test_observation.py | 42 +++++++++++++++++++++++++++++++++- 5 files changed, 101 insertions(+), 12 deletions(-) diff --git a/src/py21cmsense/_utils.py b/src/py21cmsense/_utils.py index 1b4ad4b..59cc9d5 100644 --- a/src/py21cmsense/_utils.py +++ b/src/py21cmsense/_utils.py @@ -1,11 +1,15 @@ """Utility functions for 21cmSense.""" +from __future__ import annotations + import numpy as np from astropy import units as un from astropy.coordinates import EarthLocation, SkyCoord from astropy.time import Time from pyuvdata import utils as uvutils +from . import units as tp + def between(xmin, xmax): """Return an attrs validation function that checks a number is within bounds.""" @@ -39,7 +43,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: tp.Time, + bls_enu: np.ndarray, + latitude: float, + phase_center_dec: tp.Angle | None = None, + use_apparent: bool = True, ): """Compute UVWs phased to a point rotated from zenith by a certain amount of time. @@ -50,12 +58,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,16 +83,26 @@ def phase_past_zenith( # JD is arbitrary jd = 2454600 + tm = Time(jd, format="jd") zenith_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") + if phase_center_dec is not None: + zenith_coord = SkyCoord( + ra=zenith_coord.ra, + dec=phase_center_dec, + obstime=tm, + frame="icrs", + location=telescope_location, + ) + obstimes = zenith_coord.obstime + time_past_zenith lsts = obstimes.sidereal_time("apparent", longitude=0.0).rad 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 f1bba59..33b6a99 100644 --- a/src/py21cmsense/observatory.py +++ b/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)) 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) From 5ebc804d31973fdf464bc5925afafae090367342 Mon Sep 17 00:00:00 2001 From: Steven Murray Date: Tue, 29 Oct 2024 10:48:39 +0100 Subject: [PATCH 2/6] maint: correct typing for py312 --- src/py21cmsense/_utils.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/src/py21cmsense/_utils.py b/src/py21cmsense/_utils.py index 59cc9d5..637f1de 100644 --- a/src/py21cmsense/_utils.py +++ b/src/py21cmsense/_utils.py @@ -1,6 +1,6 @@ """Utility functions for 21cmSense.""" -from __future__ import annotations +from typing import Optional import numpy as np from astropy import units as un @@ -8,8 +8,6 @@ from astropy.time import Time from pyuvdata import utils as uvutils -from . import units as tp - def between(xmin, xmax): """Return an attrs validation function that checks a number is within bounds.""" @@ -43,10 +41,10 @@ def find_nearest(array, value): @un.quantity_input def phase_past_zenith( - time_past_zenith: tp.Time, + time_past_zenith: un.hour, bls_enu: np.ndarray, latitude: float, - phase_center_dec: tp.Angle | None = None, + phase_center_dec: Optional[un.rad] = None, # noqa use_apparent: bool = True, ): """Compute UVWs phased to a point rotated from zenith by a certain amount of time. From 2ee81ac9bd4fe4c0a35c47a57b81afa1766882ae Mon Sep 17 00:00:00 2001 From: Steven Murray Date: Tue, 29 Oct 2024 11:31:32 +0100 Subject: [PATCH 3/6] fix: typing is weird for astropy --- src/py21cmsense/_utils.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/src/py21cmsense/_utils.py b/src/py21cmsense/_utils.py index 637f1de..9ed0be9 100644 --- a/src/py21cmsense/_utils.py +++ b/src/py21cmsense/_utils.py @@ -1,7 +1,5 @@ """Utility functions for 21cmSense.""" -from typing import Optional - import numpy as np from astropy import units as un from astropy.coordinates import EarthLocation, SkyCoord @@ -44,7 +42,7 @@ def phase_past_zenith( time_past_zenith: un.hour, bls_enu: np.ndarray, latitude: float, - phase_center_dec: Optional[un.rad] = None, # noqa + 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. From c3fa0ef40fcfefb982ed091a171f045cfbb794fb Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 23 Dec 2024 20:16:28 +0000 Subject: [PATCH 4/6] [pre-commit.ci] pre-commit autoupdate MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit updates: - [github.com/astral-sh/ruff-pre-commit: v0.8.3 → v0.8.4](https://github.com/astral-sh/ruff-pre-commit/compare/v0.8.3...v0.8.4) --- .pre-commit-config.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From fe21e58882974dc97c9f801cf5d0c06c5ed86e94 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Wed, 1 Jan 2025 09:35:31 +0000 Subject: [PATCH 5/6] chore(deps): bump pypa/gh-action-pypi-publish from 1.12.2 to 1.12.3 Bumps [pypa/gh-action-pypi-publish](https://github.com/pypa/gh-action-pypi-publish) from 1.12.2 to 1.12.3. - [Release notes](https://github.com/pypa/gh-action-pypi-publish/releases) - [Commits](https://github.com/pypa/gh-action-pypi-publish/compare/v1.12.2...v1.12.3) --- updated-dependencies: - dependency-name: pypa/gh-action-pypi-publish dependency-type: direct:production update-type: version-update:semver-patch ... Signed-off-by: dependabot[bot] --- .github/workflows/deploy.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From 3a84ca8bbf3585352c8acea0d57cf4af0a67dad4 Mon Sep 17 00:00:00 2001 From: DanielaBreitman Date: Wed, 8 Jan 2025 21:23:48 +0100 Subject: [PATCH 6/6] style: rename zenith_coord to phase_center_coord --- src/py21cmsense/_utils.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/src/py21cmsense/_utils.py b/src/py21cmsense/_utils.py index 9ed0be9..b8badc4 100644 --- a/src/py21cmsense/_utils.py +++ b/src/py21cmsense/_utils.py @@ -81,25 +81,25 @@ def phase_past_zenith( jd = 2454600 tm = Time(jd, format="jd") - zenith_coord = SkyCoord( + phase_center_coord = SkyCoord( alt=90 * un.deg, az=0 * un.deg, 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: - zenith_coord = SkyCoord( - ra=zenith_coord.ra, + 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__"): @@ -107,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, ) @@ -117,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))