Bandpass calibration demonstration


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()