Skip to content

Commit a99cdf1

Browse files
bmerrybmerry
committed
Improve UVW coordinate calculations
See the comments added to the code for detail about why this approach improves things. As a side benefit, the special case handling for the poles can be removed. This will change the results at the poles, in that orientation will be determined by the given RA. The golden RA/DEC in the unit tests changed. These are just a paste of what is coming out, rather than independently verified results.
1 parent 40d69c7 commit a99cdf1

2 files changed

Lines changed: 61 additions & 19 deletions

File tree

‎katpoint/target.py‎

Lines changed: 31 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -618,24 +618,42 @@ def uvw_basis(self, timestamp=None, antenna=None):
618618
dimension to the timestamp.
619619
"""
620620
timestamp, antenna = self._set_timestamp_antenna_defaults(timestamp, antenna)
621-
# NCP vector is J2000 NCP
622-
ncp = construct_radec_target(0.0, np.pi / 2.0)
623-
# Get J2000 NCP az-el vector at current epoch pointed to by reference antenna
624-
ncp_az, ncp_el = ncp.azel(timestamp, antenna)
621+
if is_iterable(timestamp) and self.body_type != 'radec':
622+
# Some calculations depend on ra/dec in a way that won't easily
623+
# vectorise.
624+
bases = [self.uvw_basis(t, antenna) for t in timestamp]
625+
return np.stack(bases, axis=-1)
626+
627+
# Offset the target slightly in declination to approximate the
628+
# derivative of ENU in the direction of increasing declination. This
629+
# used to just use the NCP, but the astrometric-to-topocentric
630+
# conversion doesn't simply rotate the celestial sphere, but also
631+
# distorts it, and so that introduced errors.
632+
#
633+
# To avoid issues close to the poles, we always offset towards the
634+
# equator. We also can't offset by too little, as ephem uses only
635+
# single precision and this method suffers from loss of precision.
636+
# 0.03 was found by experimentation (albeit on a single data set) to
637+
# to be large enough to avoid the numeric instability.
638+
if is_iterable(timestamp):
639+
# Due to the test above, this is a radec target and so timestamp
640+
# doesn't matter. But we want a scalar.
641+
ra, dec = self.radec(None, antenna)
642+
else:
643+
ra, dec = self.radec(timestamp, antenna)
644+
offset_sign = -1 if dec > 0 else -1
645+
offset = construct_radec_target(ra, dec + 0.03 * offset_sign)
646+
# Get offset az-el vector at current epoch pointed to by reference antenna
647+
offset_az, offset_el = offset.azel(timestamp, antenna)
625648
# Obtain direction vector(s) from reference antenna to target
626649
az, el = self.azel(timestamp, antenna)
627650
# w axis points toward target
628651
w = np.array(azel_to_enu(az, el))
629-
# enu vector pointing from reference antenna to north celestial pole
630-
z = np.array(azel_to_enu(ncp_az, ncp_el))
631-
# u axis is orthogonal to z and w, and row_stack makes it 2-D array of column vectors (for finding poles)
632-
u = np.row_stack(np.cross(z, w, axis=0))
652+
# enu vector pointing from reference antenna to offset point
653+
z = np.array(azel_to_enu(offset_az, offset_el))
654+
# u axis is orthogonal to z and w, and row_stack makes it 2-D array of column vectors
655+
u = np.row_stack(np.cross(z, w, axis=0)) * offset_sign
633656
u_norm = np.sqrt(np.sum(u ** 2, axis=0))
634-
# If the target is a celestial pole (so that w equals z or -z), u and v become degenerate
635-
poles = u_norm < 1e-12
636-
# Arbitrarily pick east vector of ENU system as u in this case
637-
u[:, poles] = [[1.0], [0.0], [0.0]]
638-
u_norm[poles] = 1.0
639657
# Ensure that u and w (and therefore v) have the same shape to handle scalar vs array output correctly
640658
u = u.reshape(w.shape) / u_norm
641659
v = np.cross(w, u, axis=0)

‎katpoint/test/test_target.py‎

Lines changed: 30 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -164,6 +164,7 @@ def setUp(self):
164164
self.ant1 = katpoint.Antenna('A1, -31.0, 18.0, 0.0, 12.0, 0.0 0.0 0.0')
165165
self.ant2 = katpoint.Antenna('A2, -31.0, 18.0, 0.0, 12.0, 10.0 -10.0 0.0')
166166
self.ts = katpoint.Timestamp('2013-08-14 08:25')
167+
self.uvw = [10.822902530598125, -9.103009437168168, -2.220446049250313e-16]
167168

168169
def test_coords(self):
169170
"""Test coordinate conversions for coverage."""
@@ -185,16 +186,39 @@ def test_delay(self):
185186
def test_uvw(self):
186187
"""Test uvw calculation."""
187188
u, v, w = self.target.uvw(self.ant2, self.ts, self.ant1)
188-
np.testing.assert_almost_equal(u, 10.821750916197391, decimal=5)
189-
np.testing.assert_almost_equal(v, -9.1043784587765906, decimal=5)
190-
np.testing.assert_almost_equal(w, 4.7781625336985198e-10, decimal=5)
189+
np.testing.assert_almost_equal(u, self.uvw[0], decimal=5)
190+
np.testing.assert_almost_equal(v, self.uvw[1], decimal=5)
191+
np.testing.assert_almost_equal(w, self.uvw[2], decimal=5)
191192

192193
def test_uvw_array(self):
193194
"""Test uvw calculation on an array."""
194195
u, v, w = self.target.uvw(self.ant2, np.array([self.ts, self.ts]), self.ant1)
195-
np.testing.assert_array_almost_equal(u, np.array([10.821750916197391] * 2), decimal=5)
196-
np.testing.assert_array_almost_equal(v, np.array([-9.1043784587765906] * 2), decimal=5)
197-
np.testing.assert_array_almost_equal(w, np.array([4.7781625336985198e-10] * 2), decimal=5)
196+
np.testing.assert_array_almost_equal(u, np.array([self.uvw[0]] * 2), decimal=5)
197+
np.testing.assert_array_almost_equal(v, np.array([self.uvw[1]] * 2), decimal=5)
198+
np.testing.assert_array_almost_equal(w, np.array([self.uvw[2]] * 2), decimal=5)
199+
200+
def test_uvw_array_radec(self):
201+
"""Test uvw calculation on an array when the target is a radec target."""
202+
ra, dec = self.target.radec(self.ts, self.ant1)
203+
target = katpoint.construct_radec_target(ra, dec)
204+
u, v, w = target.uvw(self.ant2, np.array([self.ts, self.ts]), self.ant1)
205+
np.testing.assert_array_almost_equal(u, np.array([self.uvw[0]] * 2), decimal=5)
206+
np.testing.assert_array_almost_equal(v, np.array([self.uvw[1]] * 2), decimal=5)
207+
np.testing.assert_array_almost_equal(w, np.array([self.uvw[2]] * 2), decimal=5)
208+
209+
def test_uvw_hemispheres(self):
210+
"""Test uvw calculation near the equator.
211+
212+
The implementation behaves differently depending on the sign of
213+
declination. This test is to catch sign flip errors.
214+
"""
215+
target1 = katpoint.construct_radec_target(0.0, -1e-9)
216+
target2 = katpoint.construct_radec_target(0.0, +1e-9)
217+
u1, v1, w1 = target1.uvw(self.ant2, self.ts, self.ant1)
218+
u2, v2, w2 = target2.uvw(self.ant2, self.ts, self.ant1)
219+
np.testing.assert_almost_equal(u1, u2, decimal=5)
220+
np.testing.assert_almost_equal(v1, v2, decimal=5)
221+
np.testing.assert_almost_equal(w1, w2, decimal=5)
198222

199223
def test_lmn(self):
200224
"""Test lmn calculation."""

0 commit comments

Comments
 (0)