import math
import numpy
from os import path
import datetime
from astropy.coordinates import Angle
from astropy import units as u
from astropy.time import Time
from casacore.measures import measures
from casacore.quanta import quantity
from casacore import tables
[docs]def time_to_quantity(in_time):
"""
I have found that for some reason I get either a string or a byte array as returned isot format from astropy.time
There seems to be no reason to it - but it does seem to be a function of the value. Some dates are always byte arrays
and some are always strings - this helper function just tests what it going on tries to return a valid casacore.quanta
:param astropy.Time in_time: astropy Time object
:returns: a casacore.quantity of the same value
"""
in_time.format = 'isot'
if len(in_time.value) > 1:
time_datetime = datetime.datetime.fromisoformat(in_time.value)
else:
time_datetime = datetime.datetime.fromisoformat(in_time.value[0])
out_time = quantity(time_datetime.isoformat())
return out_time
def _subtable(t, name, readonly=True):
return tables.table(t + "/" + name, readonly=readonly)
[docs]def get_r(h, delta) -> numpy.array:
"""
Return the matrix that will project XYZ into UVW. This is the matrix that is typically used to get the
transformation from a geocentric position to a UVW at the epoch of date.
See A.R. Thompson, J.M. Moran, and G.W. Swenson Jr.,
Interferometry and Synthesis in Radio Astronomy (equation 4.3 in the Third Edition)
:param double h: -- the Greenwich hour angle in radian
:param double delta: -- the declination in radian
:returns: A numpy array of the matrix R
:rtype: numpy.array
"""
sin_h = math.sin(h)
cos_h = math.cos(h)
sin_d = math.sin(delta)
cos_d = math.cos(delta)
matrix_r = numpy.array([[sin_h, cos_h, 0],
[-1 * sin_d * cos_h, sin_d * sin_h, cos_d],
[cos_d * cos_h, -1 * cos_d * sin_h, sin_d]])
return matrix_r
[docs]def get_uvw(xyposA: list, xyposB: list, epoch: Time, ra: Angle, decl: Angle, swap_baselines=False) -> numpy.array:
"""
Return the baseline vector for the 2 given positions in epoch of date
:param list xyposA: position of antenna 1 (vector geocentric XYZ in m)
:param list xyposB: position of antenna 2 (vector geocentric XYZ in m)
:param astropy.Time epoch: Time for which to calculate the UVW
:param astropy.Angle ra: Right Ascension (J2000)
:param astropy.Angle decl: Declination (J2000)
:param Bool swap_baselines: (xyposB - xyposA) is assumed - if the reverse is required set this to True
:returns: The uvw baseline at Epoch of Date
:rtype: numpy.array
"""
if swap_baselines:
swap = -1.0
else:
swap = 1.0
hour = epoch.sidereal_time('apparent', 'greenwich') - ra
R = get_r(hour.radian, decl.radian)
dx = swap * (xyposB[0] - xyposA[0])
dy = swap * (xyposB[1] - xyposA[1])
dz = swap * (xyposB[2] - xyposA[2])
baseline = [[dx], [dy], [dz]]
uvw = numpy.matmul(R, baseline)
return numpy.array([uvw[0], uvw[1], uvw[2]])
[docs]def get_uvw_J2000(xyposA: numpy.array, xyposB: numpy.array, refpos: numpy.array, epoch: Time, ra: Angle, decl: Angle,
position_frame='itrf', swap_baselines=False) -> numpy.array:
"""
Return the baseline vector for the 2 given positions at the J2000 Epoch. This ensures that any
image formed by the Fourier transform of a UVW cube in this frame will itself be in the ICRS frame.
Note: I believe that the ICRS and J2000 are aligned <at> J2000 - but caveat emptor until I've sorted that out.
Note: This assumes you have a reference position at which all the sidereal angles are calculated but it is more likely
that you want a "per antenna" calculation. So be sure you want this before you use it.
:param list xyposA: position of antenna 1 (vector geocentric XYZ in m)
:param list xyposB: position of antenna 2 (vector geocentric XYZ in m)
:param astropy.Time epoch: Time for which to calculate the UVW
:param astropy.Angle ra: Right Ascension (J2000)
:param astropy.Angle decl: Declination (J2000)
:param str position_frame: The frame of reference for the positions ITRF is default - WGS84 is also supported
:param Bool swap_baselines: (xyposB - xyposA) is assumed - if the reverse is required set this to True
:returns: The uvw baseline at J2000
:rtype: numpy.array [u,v,w]
"""
if swap_baselines:
swap = -1.0
else:
swap = 1.0
dm = measures()
refant = dm.position(position_frame, quantity(refpos[0], 'm'), quantity(refpos[1], 'm'), quantity(refpos[2], 'm'))
refant_itrf = dm.measure(refant, 'itrf')
"""Reference Position"""
dm.do_frame(
refant_itrf) # where are we probably should just be the observatory - although antenna based UVW are probably best
"""Reference Direction"""
ra_measure = quantity(ra.to_string(unit=u.rad, decimal=True) + " rad")
dec_measure = quantity(decl.to_string(unit=u.rad, decimal=True) + " rad")
source = dm.direction('j2000', ra_measure, dec_measure)
dm.do_frame(source) # where are we looking
"""Time"""
epoch_val = time_to_quantity(epoch)
epoch_measure = dm.epoch('ut1', epoch_val)
dm.do_frame(epoch_measure) # what time is it
"""Baseline Epoch of Date"""
dx = quantity(swap * (xyposB[0] - xyposA[0]), 'm')
dy = quantity(swap * (xyposB[1] - xyposA[1]), 'm')
dz = quantity(swap * (xyposB[2] - xyposA[2]), 'm')
b = dm.baseline('itrf', dx, dy, dz)
"""Get UVW in the J2000 reference frame"""
d = dm.to_uvw(b)
uvw = numpy.array(d.get("xyz").get_value('m'))
return uvw
[docs]def get_antenna_uvw(xyposA: numpy.array, epoch: Time, ra: Angle, decl: Angle, position_frame='itrf',
epoch_frame='J2000', swap_baselines=False) -> numpy.array:
"""
Return the geocentric uvw vector for the the given position at the given epoch.
Note: This actually just calls get_uvw_J2000 with the geocentre as the other antenna and this antenna as the reference
position.
Note: This is the per antenna calculation that you probably want
:param list xyposA: position of the antenna (vector geocentric XYZ in m)
:param astropy.Time epoch: Time for which to calculate the UVW
:param str epoch_frame: Whether you want a catalog (J2000) frome or epoch of date (default: 'J2000')
:param str position_frame: The frame of the input positions
(generally WGS84 or ITRF) we are using casacore for
this and these are the two frames supported
:param astropy.Angle ra: Right Ascension (J2000)
:param astropy.Angle decl: Declination (J2000)
:param Bool swap_baselines: (xyposB - xyposA) is assumed - if the reverse is required set this to True
:returns: The geocentric UVW baseline
:rtype: numpy.array [u,v,w]
"""
xyposB = numpy.array([0, 0, 0])
if epoch_frame == 'J2000' or epoch_frame == 'j2000':
uvw = get_uvw_J2000(xyposA, xyposB, xyposA, epoch, ra, decl, position_frame, swap_baselines)
else:
uvw = get_uvw(xyposA, xyposB, epoch, ra, decl, swap_baselines)
return uvw