In [ ]:
%matplotlib inline
import os
import sys
sys.path.append(os.path.join('..', '..'))
from data_models.parameters import arl_path
results_dir = arl_path('test_results')
from matplotlib import pylab
import numpy
from astropy.coordinates import SkyCoord
from astropy import units as u
from astropy.wcs.utils import pixel_to_skycoord
from matplotlib import pyplot as plt
from wrappers.serial.visibility.base import create_blockvisibility
from wrappers.serial.calibration.operations import apply_gaintable
from wrappers.serial.visibility.operations import copy_visibility
from wrappers.serial.calibration.calibration import solve_gaintable
from wrappers.serial.visibility.coalesce import convert_blockvisibility_to_visibility, \
convert_visibility_to_blockvisibility
from wrappers.serial.calibration.operations import create_gaintable_from_blockvisibility
from wrappers.serial.image.operations import show_image
from wrappers.serial.simulation.testing_support import create_test_image, simulate_gaintable
from wrappers.serial.simulation.configurations import create_named_configuration
from wrappers.serial.imaging.base import create_image_from_visibility
from workflows.serial.imaging.imaging_serial import predict_list_serial_workflow
from data_models.polarisation import PolarisationFrame
pylab.rcParams['figure.figsize'] = (8.0, 8.0)
pylab.rcParams['image.cmap'] = 'rainbow'
import logging
log = logging.getLogger()
log.setLevel(logging.DEBUG)
log.addHandler(logging.StreamHandler(sys.stdout))
Construct LOW core configuration
In [ ]:
lowcore = create_named_configuration('LOWBD2-CORE')
We create the visibility. This just makes the uvw, time, antenna1, antenna2, weight columns in a table
In [ ]:
times = numpy.zeros([1])
vnchan = 128
frequency = numpy.linspace(0.8e8, 1.2e8, vnchan)
channel_bandwidth = numpy.array(vnchan*[frequency[1]-frequency[0]])
phasecentre = SkyCoord(ra=+15.0 * u.deg, dec=-45.0 * u.deg, frame='icrs', equinox='J2000')
bvt = create_blockvisibility(lowcore, times, frequency, channel_bandwidth=channel_bandwidth,
weight=1.0, phasecentre=phasecentre, polarisation_frame=PolarisationFrame('stokesI'))
Read the venerable test image, constructing an image
In [ ]:
m31image = create_test_image(frequency=frequency, cellsize=0.0005)
nchan, npol, ny, nx = m31image.data.shape
m31image.wcs.wcs.crval[0] = bvt.phasecentre.ra.deg
m31image.wcs.wcs.crval[1] = bvt.phasecentre.dec.deg
m31image.wcs.wcs.crpix[0] = float(nx // 2)
m31image.wcs.wcs.crpix[1] = float(ny // 2)
fig=show_image(m31image)
Predict the visibility from this image
In [ ]:
vt = convert_blockvisibility_to_visibility(bvt)
vt = predict_list_serial_workflow(bvt, m31image, context='timeslice')
bvt = convert_visibility_to_blockvisibility(vt)
Create a gain table with modest amplitude and phase errors, smootheed over 16 channels
In [ ]:
gt = create_gaintable_from_blockvisibility(bvt)
gt = simulate_gaintable(gt, phase_error=1.0, amplitude_error=0.1, smooth_channels=16)
Plot the gains applied
In [ ]:
plt.clf()
for ant in range(4):
amp = numpy.abs(gt.gain[0,ant,:,0,0])
plt.plot(amp)
plt.title('Amplitude of bandpass')
plt.xlabel('channel')
plt.show()
plt.clf()
for ant in range(4):
phase = numpy.angle(gt.gain[0,ant,:,0,0])
plt.plot(phase)
plt.title('Phase of bandpass')
plt.xlabel('channel')
plt.show()
In [ ]:
cbvt = copy_visibility(bvt)
cbvt = apply_gaintable(cbvt, gt)
Solve for the gains
In [ ]:
gtsol=solve_gaintable(cbvt, bvt, phase_only=False)
Plot the solved relative to the applied. Declare antenna 0 to be the reference.
In [ ]:
plt.clf()
for ant in range(4):
amp = numpy.abs(gtsol.gain[0,ant,:,0,0]/gt.gain[0,ant,:,0,0])
plt.plot(amp)
plt.title('Relative amplitude of bandpass')
plt.xlabel('channel')
plt.show()
plt.clf()
for ant in range(4):
refphase = numpy.angle(gtsol.gain[0,0,:,0,0]/gt.gain[0,0,:,0,0])
phase = numpy.angle(gtsol.gain[0,ant,:,0,0]/gt.gain[0,ant,:,0,0])
plt.plot(phase-refphase)
plt.title('Relative phase of bandpass')
plt.xlabel('channel')
plt.show()