<jupyter_start><jupyter_text>Bandpass calibration demonstration<jupyter_code>%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))<jupyter_output><empty_output><jupyter_text>Construct LOW core configuration<jupyter_code>lowcore = create_named_configuration('LOWBD2-CORE')<jupyter_output><empty_output><jupyter_text>We create the visibility. This just makes the uvw, time, antenna1, antenna2, weight columns in a table<jupyter_code>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'))<jupyter_output><empty_output><jupyter_text>Read the venerable test image, constructing an image<jupyter_code>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)<jupyter_output><empty_output><jupyter_text>Predict the visibility from this image<jupyter_code>vt = convert_blockvisibility_to_visibility(bvt)
vt = predict_list_serial_workflow(bvt, m31image, context='timeslice')
bvt = convert_visibility_to_blockvisibility(vt)<jupyter_output><empty_output><jupyter_text>Create a gain table with modest amplitude and phase errors, smootheed over 16 channels<jupyter_code>gt = create_gaintable_from_blockvisibility(bvt)
gt = simulate_gaintable(gt, phase_error=1.0, amplitude_error=0.1, smooth_channels=16)<jupyter_output><empty_output><jupyter_text>Plot the gains applied<jupyter_code>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()


cbvt = copy_visibility(bvt)
cbvt = apply_gaintable(cbvt, gt)<jupyter_output><empty_output><jupyter_text>Solve for the gains<jupyter_code>gtsol=solve_gaintable(cbvt, bvt, phase_only=False)<jupyter_output><empty_output><jupyter_text>Plot the solved relative to the applied. Declare antenna 0 to be the reference.<jupyter_code>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()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1165155869997789, 'fraction_numerical': 0.021445942958213574, 'mean_word_length': 7.424581005586592, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 37, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '70957', 'license': 'Apache-2.0', 'n_tokens_mistral': 1682, 'n_tokens_neox': 1591, 'n_words': 282, 'path': 'workflows/notebooks/bandpass-calibration_serial.ipynb', 'repo_name': 'ska-telescope/algorithm-reference-library', 'size': '4523'}