Skip to content

Instantly share code, notes, and snippets.

@rstofi

rstofi/considerations_of_PB_correction_for_MeerKLASS_OTFM.py

Last active December 7, 2023 17:06
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save rstofi/1d73ac5a785147bd887856050022a3c0 to your computer and use it in GitHub Desktop.
Save rstofi/1d73ac5a785147bd887856050022a3c0 to your computer and use it in GitHub Desktop.
Download ZIP
Code used for a MeerKLASS memo on the importance of OTFM PB correction and the amplitude errors introduced
Raw
considerations_of_PB_correction_for_MeerKLASS_OTFM.py
"""Simple code to re-create OTFM beam smearing plots. For the default plots, and
mathematical background see:
https://iopscience.iop.org/article/10.3847/1538-4357/aaef7c/pdf
Also, the effect is computed and results are re-computed for MeerKAT. Both for
constant RA and constant elevation (MeerKLASS)
The code is based on a circular GAussian beam input for now.
I try to vectorise things rather than parallelize and hope the underlying C code
is fast enough for these things
MeerKLASS general parameters:
https://skaafrica.atlassian.net/wiki/spaces/ESDKB/pages/277315585/MeerKAT+specifications
Author: Kristof Rozgonyi
e-mail: rstofi@gmail.com
License: MIT
"""
#=== Imports ===
import os, sys
import numpy as np
from scipy.special import erf
import matplotlib
import matplotlib.pyplot as plt
from matplotlib.offsetbox import AnchoredText
from matplotlib.patheffects import withStroke
#Astropy
from astropy.coordinates import Longitude, Latitude, EarthLocation
from astropy.coordinates import SkyCoord, AltAz
from astropy.time import Time
from astropy import units as u
from astropy.coordinates import AltAz
#=== Globals ===
#RCparams for plotting
matplotlib.rcParams['xtick.direction'] = 'in'
matplotlib.rcParams['ytick.direction'] = 'in'
matplotlib.rcParams['xtick.major.size'] = 9
matplotlib.rcParams['ytick.major.size'] = 9
matplotlib.rcParams['xtick.major.width'] = 3
matplotlib.rcParams['ytick.major.width'] = 3
matplotlib.rcParams['xtick.minor.size'] = 6
matplotlib.rcParams['ytick.minor.size'] = 6
matplotlib.rcParams['xtick.minor.width'] = 2
matplotlib.rcParams['ytick.minor.width'] = 2
matplotlib.rcParams['axes.linewidth'] = 2
plt.rcParams['xtick.labelsize']=16
plt.rcParams['ytick.labelsize']=16
#4 sampled colors from viridis
c0 = '#440154';#Purple
c1 = '#30678D';#Blue
c2 = '#35B778';#Greenish
c3 = '#FDE724';#Yellow
# === Telescope related
# The wideband coarse (4k)
N_CHAN = 410 #I average every 10 channels together not to simulate all 4096chans
L_band_PIXELSIZE = 12 # Arcseconds
L_BAND_NU_CENTRAL = 1.285
L_BAND_CW = 0.00209 # 209kHz channelwidth (208.984 kHz)
UHF_band_PIXELSIZE = 14 # Arcseconds
UHF_BAND_NU_CENTRAL = 0.7975
UHF_BAND_CW = 0.000133 # 132.812 kHz
#=== Functions ===
def deg_to_rad(deg):
return deg * (np.pi / 180)
def arcsecond_to_arcmin(arcsec):
return float(arcsec) / 60.
def arcmin_to_deg(arcmin):
return float(arcmin) / 60.
def arcsecond_to_deg(arcsec):
return float(arcsec) / 360.
def FWHM_to_std(FWHM):
"""Convert Gaussian FWHM to standard deviation
https://en.wikipedia.org/wiki/Full_width_at_half_maximum
"""
return FWHM / 2.355
def std_to_FWHM(std):
"""Convert Gaussian FWHM to standard deviation
https://en.wikipedia.org/wiki/Full_width_at_half_maximum
"""
return std * 2.355
def get_pixel_coordinate_array(N,dn):
"""Quick subroutine to compute the grid cell midpoints around zero (one point
is centered at zero) with N points with dn width. If N is even more negative
pints will be computed as this is preffered fro numpy FFT
"""
return np.linspace(-dn*(np.floor(N/2)),dn*np.ceil((N/2-1)),N,endpoint=True)
def get_VLA_PB_FWHM_from_freq(nu):
"""Compute the VLA PB FWHM based on Pereley VLA memo 195 and the
Mooley et al. 2018 paper
Parameters
==========
nu: float
The given frequency in GHz
Returns
=======
The PB FWHM in arcseconds
"""
return (42 * 60) / nu
def get_MeerKAT_PB_FWHM_from_freq(nu):
"""Compute the MeerKAT PB FWHM based on eq. 4 and 5 from Mauch et al. 2019
(deep2 paper). I use the Horizontal PB model scaling
Parameters
==========
nu: float
The given frequency in GHz
Returns
=======
The PB FWHM in arcseconds
"""
return (86.25 * 60) / nu #i.e. (57.5 * 1.5) / nu
def circular_gaussian_beam_value(x, y, sigma=0.5):
"""A simple 2D Gaussian beam value computed on a given image point
Parameters
==========
x: <numpy.ndarray>
The x (RA) directional coordinates relative to the pointing centre
y: <numpy.ndarray>
The y (Dec) directional coordinates relative to the pointing centre
sigma: float, optional
The standard deviation of the circular Gaussian
Returns
=======
z: float
The primary beam value at [x,y]
"""
if np.shape(x) != np.shape(y):
raise ValueError('x and y have different shapes!')
z = np.exp(-(np.power(x,2) + np.power(y,2)) / (2 * np.power(sigma,2)))
return z
def RA_smeared_gaussian_beam_value(x, y, dx, sigma=0.5):
"""Compute the smeared beam (eq. 10) from
https://iopscience.iop.org/article/10.3847/1538-4357/aaef7c/pdf
Parameters
==========
x: <numpy.ndarray>
The x (RA) directional coordinates relative to the pointing centre
y: <numpy.ndarray>
The y (Dec) directional coordinates relative to the pointing centre
dx: float
The max RA (x) coordinate offset of the beam for a given phase centre
sigma: float, optional
The standard deviation of the circular Gaussian
Returns
=======
z: float
The primary beam value at [x,y]
"""
if np.shape(x) != np.shape(y):
raise ValueError('x and y have different shapes!')
if dx == 0:
return circular_gaussian_beam_value(x, y, sigma=sigma)
else:
b0 = np.exp(-(np.power(x,2) + np.power(y,2)) / (2 * np.power(sigma,2)))
A = (sigma / dx) * b0 * np.exp(np.power(x,2) / (2 * np.power(sigma,2))) \
* np.sqrt(np.pi / 2)
z = A * (erf((x + (dx / 2)) / (np.sqrt(2) * sigma)) \
- erf((x - (dx / 2)) / (np.sqrt(2) * sigma)))
return z
def general_smeared_gaussian_beam_value(x, y, dx, dy, sigma=0.5):
"""Compute the general smeared beam computed using a linear approximation
See the memo for the formula and derivation.
Parameters
==========
x: <numpy.ndarray>
The x (RA) directional coordinates relative to the pointing centre
y: <numpy.ndarray>
The y (Dec) directional coordinates relative to the pointing centre
dx: float
The max RA (x) coordinate offset of the beam for a given phase centre
in arcseconds
dx: float
The max Dec (y) coordinate offset of the beam for a given phase centre
in arcseconds
.
sigma: float, optional
The standard deviation of the circular Gaussian
Returns
=======
z: float
The primary beam value at [x,y]
"""
if np.shape(x) != np.shape(y):
raise ValueError('x and y have different shapes!')
if dx == 0 and dy == 0:
return circular_gaussian_beam_value(x, y, sigma=sigma)
else:
ds = np.sqrt(np.power(dx,2) + np.power(dy,2))
A = (sigma / ds) * np.sqrt(np.pi / 2)
s_exp = np.exp(-(np.power(((x * (dy / ds)) + (y * (dx / ds))),2)) \
/ (2 * np.power(sigma,2)))
s_smear = (erf(((x * (dx / ds)) - (y * (dy / ds)) + (ds / 2)) \
/ (np.sqrt(2) * sigma)) \
- erf(((x * (dx / ds)) - (y * (dy / ds)) - (ds / 2)) \
/ (np.sqrt(2) * sigma)))
z = A * s_exp * s_smear
return z
def compute_beam_image(N, dn,
beam_model='circular_gaussian',
sigma=0.5,
dx=None,
dy=None):
"""Generate a regular image of NxN with dn pixel size and generate the given
beam model.
The beam values are computed at the pixel centre points
Parameters
==========
N: int
The size of the sky image
dn: int
The image pixel size
beam_model: str, optional
Name the primary beam model from the following list:
['circular_gaussian', 'RA_smeared_gaussian', 'general_smeared_gaussian']
sigma: float, optional
The standard deviation of the circular Gaussian
dx: float, optional
The max RA (x) coordinate offset of the beam for a given phase centre (scan)
dy: float, optional
The max RA (x) coordinate offset of the beam for a given phase centre (scan)
Returns
=======
x: <numpy.ndarray>
The pixel x coordinates (RA) as a grid (NxN matrix with constant columns)
y: <numpy.ndarray>
The pixel y coordinates(Dec) as a grid (NxN matrix with constant rows)
z: <numpy.ndarray>
The primary beam values as a grid (NxN matrix)
"""
#This is a dirty type solution
N = int(np.ceil(N))
#Extra end at negative values => This is good for numpy FFT
x = get_pixel_coordinate_array(N,dn)
y = get_pixel_coordinate_array(N,dn)
x, y = np.meshgrid(x,y)
if beam_model == 'circular_gaussian':
z = circular_gaussian_beam_value(x,y,sigma=sigma)
elif beam_model == 'RA_smeared_gaussian':
if dx == None:
raise ValueError('no dx is defined!')
z = RA_smeared_gaussian_beam_value(x,y,sigma=sigma,dx=dx)
elif beam_model == 'general_smeared_beam':
if dx == None:
raise ValueError('no dx is defined!')
if dy == None:
raise ValueError('no dy is defined!')
z = general_smeared_gaussian_beam_value(x,y,sigma=sigma,dx=dx,dy=dy)
else:
raise NotImplementedError('Beam model is wrong or not implemented yet')
#else:
# z = circular_gaussian_beam_value(x,y,sigma=sigma)
# z = z - RA_smeared_gaussian_beam_value(x,y,sigma=sigma,dx=dx)
return x, y, z
def get_beam_slice(x, y, z, slice_dir='x', slice_at=0.):
"""Get a slice from a PB image at constant RA or Dec
Parameters
==========
x: <numpy.ndarray>
The pixel x coordinates (RA) as a grid (NxN matrix with constant columns)
y: <numpy.ndarray>
The pixel y coordinates(Dec) as a grid (NxN matrix with constant rows)
z: <numpy.ndarray>
The primary beam values as a grid (NxN matrix)
slice_dir: string, optional
Slice direction i.e. the direction which the slice is returned
ca be 'RA', 'x' or 'Dec', 'y'
e.g if slice_dir = x a slice at constant y is returned
slice_at: float, optional
The constant RA or Dec value in which the slice is taken
Returns
=======
xy_slice: <numpy.ndarray>
The slice pixel (x or y) coordinates (1D array)
z_slice: <numpy.ndarray>
The beam values at the slice (1D array)
"""
#I skip checking size matches
if slice_dir == 'RA' or slice_dir == 'x':
slice_dir_array = x[0,:]
perpendicular_dir_array = y
elif slice_dir == 'Dec' or slice_dir == 'y':
slice_dir_array = y[:,0]
perpendicular_dir_array = x
else:
raise ValueError('Invalid slice direction is given!')
#Get the closest value in array to the given slice at
#I use the fact that the bem is computed on an NxN regular image
slice_at = slice_dir_array[np.argmin(np.abs(
slice_dir_array - slice_at))]
xy_slice = np.copy(slice_dir_array)
z_slice = np.copy(z[perpendicular_dir_array == slice_at])
return xy_slice, z_slice
def fractional_PB_cange_of_single_slice(pixelsize=1,
PB_FWHM=60,
Ndx = 100,
fdx = 0.01,
slice_dir='x',
slice_at=0.,
N_image=None,
beam_model = 'RA_smeared_gaussian'
):
"""Generates fig 2 data matrix from Mooley et al. 2018
Currently only working with the RA smeared beam model
TO DO: make this function more general to work with an arbritary PB model
Parameters
==========
pixelsize: float, opt
The image pixelsize in arcseconds
PB_FWHM: float, opt
The Primary Beam FWHM in arcseconds. Used to calculate the std of the PB
Ndx: int, optional
The number of dx coordinate offset values to test
fdx: float, optional
The fractional change of the x coordinate for a single subscan (one phase
centre) from 0 to Ndx * fdx *PB_FWHM
slice_dir: string, optional
Slice direction i.e. the direction which the slice is returned
ca be 'RA', 'x' or 'Dec', 'y'
e.g if slice_dir = x a slice at constant y is returned
slice_at: float, optional
The constant RA or Dec value in which the slice is taken
N_image: int, optional
If None the image size is 2 * PB_FWHM, otherwise this parameter controls
the image sixe (x axis of the plot)
beam_model: str, optional
Name the (smeared) primary beam model from the following list:
['RA_smeared_gaussian']
Returns
=======
dPB_array: <numpy.ndarray>
The pixel slice coordinates as a grid (N_image x N_dx matrix with
constant columns)
dx_array: <numpy.ndarray>
The generated dx value coordinates as a grid (N_image x N_dx matrix with
constant rows)
PB_err: <numpy.ndarray>
The fractional bema change values as a grid (N_image x N_dx matrix)
"""
#Set up the parameters
PB_sigma = FWHM_to_std(PB_FWHM)
if N_image == None:
N_image = np.ceil((PB_FWHM * 2) / pixelsize)
N_image = int(N_image)
ddx = fdx * PB_FWHM
#Get the output arrays
dx_array = np.arange(ddx,(Ndx * ddx) + 1,ddx)
dPB_array = get_pixel_coordinate_array(N_image,pixelsize)
PB_err = np.zeros((Ndx, N_image))
#Get the reference beam
x0, y0, z0 = compute_beam_image(N = N_image,
dn = pixelsize,
beam_model = 'circular_gaussian',
sigma = PB_sigma)
ref_xy_slice, ref_z_slice = get_beam_slice(x0,y0,z0,
slice_dir=slice_dir, slice_at=slice_at)
#Fill up the PB_err matrix
for i in range(0,Ndx):
x, y, z = compute_beam_image(N = N_image,
dn = pixelsize,
beam_model = beam_model,
sigma = PB_sigma,
dx = dx_array[i])
xy_slice, z_slice = get_beam_slice(x,y,z,
slice_dir=slice_dir, slice_at=slice_at)
PB_err[i,:] = np.divide(z_slice,ref_z_slice)
#Compute the meshgrid coordinates
dPB_array, dx_array = np.meshgrid(dPB_array,dx_array)
del x0, y0, z0, x, y, z, ref_xy_slice, ref_z_slice, xy_slice, z_slice
return dPB_array, dx_array, PB_err
def fractional_PB_cange_with_frequency_of_single_slice(pixelsize=1,
N_image = (14 * 60 * 1.6) + 1,
nu_central = 3,
Nnu = 201,
dnu = 0.01,
PB_scaling_model = 'VLA',
dx = 240,
slice_dir = 'x',
slice_at = 0.,
beam_model = 'RA_smeared_gaussian'
):
"""Generates fig 2 data matrix from Mooley et al. 2018
Currently only working with the RA smeared beam model
TO DO: make this function more general to work with an arbritary PB model
Parameters
==========
pixelsize: float, opt
The image pixelsize in arcseconds
N_image: int, optional
The image size
nu_central: float, opt
The central frequency in GHz
Nnu: int, optional
The number of frequency channels to test
dnu: float, opt
The width of a single frequency channel in GHz
PB_scaling_model: str, opt
The model of how the PB FWHM depends on the observing frequency. Can be
any from the following list: ['VLA', 'MeerKAT']
dx: float, optional
The absoolute change of the x coordinate for a single subscan used for
all frequencies (in arcseconds)
slice_dir: string, optional
Slice direction i.e. the direction which the slice is returned
ca be 'RA', 'x' or 'Dec', 'y'
e.g if slice_dir = x a slice at constant y is returned
slice_at: float, optional
The constant RA or Dec value in which the slice is taken
beam_model: str, optional
Name the (smeared) primary beam model from the following list:
['RA_smeared_gaussian']
Returns
=======
dPB_array: <numpy.ndarray>
The pixel slice coordinates as a grid (N_image x N_dx matrix with
constant columns)
dx_array: <numpy.ndarray>
The generated dx value coordinates as a grid (N_image x N_dx matrix with
constant rows)
PB_err: <numpy.ndarray>
The fractional bema change values as a grid (N_image x N_dx matrix)
"""
N_image = int(np.ceil(N_image))
#Get the output arrays
nu_array = np.add(get_pixel_coordinate_array(Nnu,dnu), nu_central)
dPB_array = get_pixel_coordinate_array(N_image,pixelsize)
PB_err = np.zeros((Nnu, N_image))
#Get the frequency deendent PB FWHM array
PB_FWHM_array = np.zeros(Nnu)
if PB_scaling_model == 'VLA':
for i in range(0,Nnu):
PB_FWHM_array[i] = get_VLA_PB_FWHM_from_freq(nu_array[i])
elif PB_scaling_model == 'MeerKAT':
for i in range(0,Nnu):
PB_FWHM_array[i] = get_MeerKAT_PB_FWHM_from_freq(nu_array[i])
#Fill up the PB_err matrix
for i in range(0,Nnu):
#Get the reference beam
x0, y0, z0 = compute_beam_image(N = N_image,
dn = pixelsize,
beam_model = 'circular_gaussian',
sigma = FWHM_to_std(PB_FWHM_array[i]))
ref_xy_slice, ref_z_slice = get_beam_slice(x0,y0,z0,
slice_dir=slice_dir, slice_at=slice_at)
#Compute the smeared beam
x, y, z = compute_beam_image(N = N_image,
dn = pixelsize,
beam_model = beam_model,
sigma = FWHM_to_std(PB_FWHM_array[i]),
dx = dx)
xy_slice, z_slice = get_beam_slice(x,y,z,
slice_dir=slice_dir, slice_at=slice_at)
PB_err[i,:] = np.divide(z_slice,ref_z_slice)
#Compute the meshgrid coordinates
dPB_array, nu_array = np.meshgrid(dPB_array,nu_array)
del x0, y0, z0, x, y, z, ref_xy_slice, ref_z_slice, xy_slice, z_slice
return dPB_array, nu_array, PB_err
def get_smearing_from_altaz_subscan(t0 = Time('2019-02-24T19:48:28'),
alt0 = 42.,
az0 = 59.6,
slew_velocity = 5.,
slew_direction = -1,
dt = 2.,
telescope = 'MeerKAT',
separation_only=False):
"""Compute the smeared beam from the altaz parameters pelescope pointing
position and sub-scan integration time
==========
t0: <astropy.Time>, optional
The the observation time in UTC with y-m-dTh:m:s. format at the time
centroid of the sub-scan
alt0: float, optional
The altitude of the observation at the phase centre
az0: float, optional
The azimuth at the observation at the phase centre
slew_velocity: float, optional
The slew velocity in arcseconds (!). Defined in the azimutal direction
slew_direction: int, optional
The direction of slewing. Can be [+1,-1]
dt: float, optional
The sub-scan integration time
telescope: str, optional
The telescope used for the subscan simulations. Possible values are:
['MeerKAT']
separetion_only: bool, optional
If true, onlt the Delta s angular separation is returned
Returns
=======
dRA: float
The RA separation for the sub-scan (corrected for Dec) in arcseconds
dDec: float
The Dec separation for the sub-scan in arcseconds
"""
if telescope == 'MeerKAT':
location = EarthLocation.from_geodetic(
Longitude('21:26:38.0', u.degree, wrap_angle=180. * u.degree, copy=False),
Latitude('-30:42:39.8', u.degree, copy=False),
height=u.Quantity(1086.6, u.m, copy=False))
else:
raise NotImplementedError('Only MeerKAT telescope is implemented yet')
#Convert slew velocity to degrees
slew_velocity_deg = arcmin_to_deg(slew_velocity)
dt_half = (dt / 2)
t1 = t0 - (dt_half * u.second)
t2 = t0 + (dt_half * u.second)
az1 = az0 - (slew_direction * np.multiply(slew_velocity_deg, dt_half))
az2 = az0 + (slew_direction * np.multiply(slew_velocity_deg, dt_half))
centre_icrs_coordinates = SkyCoord(alt = alt0 * u.deg,
az = az0 * u.deg,
obstime = t0,
frame = 'altaz',
location = location).transform_to('icrs')
start_icrs_coordinates = SkyCoord(alt = alt0 * u.deg,
az = az1 * u.deg,
obstime = t1,
frame = 'altaz',
location = location).transform_to('icrs')
end_icrs_coordinates = SkyCoord(alt = alt0 * u.deg,
az = az2 * u.deg,
obstime = t2,
frame = 'altaz',
location = location).transform_to('icrs')
dra = (end_icrs_coordinates.ra-start_icrs_coordinates.ra).arcsec \
* np.cos(centre_icrs_coordinates.dec.rad)
dDec = (end_icrs_coordinates.dec-start_icrs_coordinates.dec).arcsec
sep = end_icrs_coordinates.separation(start_icrs_coordinates).arcsec
if separation_only:
return sep
else:
return dra, dDec
#=== Survey and scanning pattern defining functions ===
def get_stripe_coordinate_values_constant_az_slew(t0 = Time('2019-02-24T19:48:28'),
alt0 = 42.,
az0 = 59.6,
slew_velocity = 5.,
slew_direction = -1,
dt = 2.,
T_stripe = 216.,
t0_rel = 0.,
telescope = 'MeerKAT',
frame = 'icrs'):
"""This function returns a relative time array and the corresponding RA and
Dec arrays for a given sstripe.
The stripe slew is performed with a constant altitude (elevation) pointing and so
slewing only in the azimuthal direction.
Parameters
==========
t0: <astropy.Time>, optional
The start of the observation time in UTC with y-m-dTh:m:s. format
alt0: float, optional
The altitude of the observation at t0
az0: float, optional
The azimuth at the observation at t0
slew_velocity: float, optional
The slew velocity in arcseconds (!). Defined in the azimutal direction
slew_direction: int, optional
The direction of slewing. Can be [+1,-1]
dt: float, optional
The visibility sampling time in seconds
T_stripe: float, optional
The time which under a single stripe is performed
t0_rel: float, optional
The relative start time of the observation (used for several stripe)
telescope: str, optional
The telescope used for the stripe simulations. Possible values are:
['MeerKAT']
frame: stra, optional
Available frames: ['icrs','altaz']
Returns
=======
t_stripe: <numpy.ndarray>
The time values of the stripe relative to the observation's starting
time in seconds.
x_stripe: <numpy.ndarray>
The RA (icrs) or Alt (altaz) values of the stripe in degrees
y_stripe: <numpy.ndarray>
The Dec (icrs) or Az (altaz) values of the stripe in degrees
ti: <astropy.Time>
The last time sample's UTC time stamp (to be used as an input for the next
stripe)
azi: float
The last ime sample's azimuth in degrees (to be used as an input for
the next stripe)
"""
if telescope == 'MeerKAT':
location = EarthLocation.from_geodetic(
Longitude('21:26:38.0', u.degree, wrap_angle=180. * u.degree, copy=False),
Latitude('-30:42:39.8', u.degree, copy=False),
height=u.Quantity(1086.6, u.m, copy=False))
else:
raise NotImplementedError('Only MeerKAT telescope is implemented yet')
#Convert slew velocity to degrees
slew_velocity_deg = arcmin_to_deg(slew_velocity)
#Get the relative time array
t_stripe = np.arange(t0_rel,T_stripe,dt)
#Get alt az vectors
alt_stipe = alt0 * np.ones(np.size(t_stripe))
az_stipe = az0 + (slew_direction * np.multiply(slew_velocity_deg, t_stripe))
t_obs_stripe = t0 + t_stripe * u.second
#Convert to RA Dec
stripe_altaz_coordinates = SkyCoord(alt = alt_stipe * u.deg,
az = az_stipe * u.deg,
obstime = t_obs_stripe,
frame = 'altaz',
location = location)
if frame == 'altaz':
x_stripe = stripe_altaz_coordinates.alt.deg
y_stripe = stripe_altaz_coordinates.az.deg
elif frame == 'icrs':
stripe_icrs_coordinates = stripe_altaz_coordinates.transform_to('icrs')
x_stripe = stripe_icrs_coordinates.ra.deg
y_stripe = stripe_icrs_coordinates.dec.deg
else:
raise NotImplementedError('Only icrs and altaz frames are implemented')
ti = t_obs_stripe[-1]
#alti = stripe_altaz_coordinates[-1].alt.deg
azi = stripe_altaz_coordinates[-1].az.deg
return t_stripe, x_stripe, y_stripe, ti, azi
def get_field_scan_coordinate_values_constant_az_slew(N_scan = 26,
t0 = Time('2019-02-24T19:48:28'),
alt0 = 42.,
az0 = 59.6,
slew_velocity = 5.,
slew_direction = -1,
dt = 2.,
T_stripe = 216.,
t0_rel = 0.,
telescope = 'MeerKAT',
frame = 'icrs'):
"""Create the full scan of the target field, defined via a starting stripe
and the number of consequent stripe.
Parameters
==========
N_scan: int, optional,
The number of stripes used to cover the field simulated
t0: <astropy.Time>, optional
The start of the observation time in UTC with y-m-dTh:m:s. format
alt0: float, optional
The altitude of the observation at t0
az0: float, optional
The azimuth at the observation at t0
slew_velocity: float, optional
The slew velocity in arcseconds (!). Defined in the azimutal direction
slew_direction: int, optional
The direction of slewing. Can be [+1,-1]
dt: float, optional
The visibility sampling time in seconds
T_stripe: float, optional
The time which under a single stripe is performed
t0_rel: float, optional
The relative start time of the observation (used for several stripes)
telescope: str, optional
The telescope used for the stripe simulations. Possible values are:
['MeerKAT']
frame: stra, optional
Available frames: ['icrs','altaz']
Returns
=======
t_field: <numpy.ndarray>
The time values of the full field scans relative to the observation's
starting time in seconds.
x_field: <numpy.ndarray>
The RA (icrs) or Alt (altaz) values of the full field scans in degrees
y_field: <numpy.ndarray>
The Dec (icrs) or Az (altaz) values of the full field scans in degrees
"""
#Define the output arrays
t_field = np.arange(t0_rel,T_stripe * N_scan,dt)
x_field = np.zeros(np.size(t_field))
y_field = np.zeros(np.size(t_field))
#Loop through the stripes
ti = t0
azi = az0
for i in range(0,N_scan):
t_stripe, x_stripe, y_stripe, ti, azi = \
get_stripe_coordinate_values_constant_az_slew(
t0 = ti,
alt0 = alt0,
az0 = azi,
slew_velocity = slew_velocity,
slew_direction = slew_direction,
dt = dt,
T_stripe = T_stripe,
t0_rel = t0_rel,
telescope = telescope,
frame = frame)
#ti_rel += T_stripe * dt
slew_direction *= -1
stripe_size = np.size(t_stripe) #This is ugly but I am mad atm for somethiong else so I dont care...
counter = i * stripe_size
t_field[counter:counter + stripe_size] = t_stripe + (i * T_stripe)
x_field[counter:counter + stripe_size] = x_stripe
y_field[counter:counter + stripe_size] = y_stripe
return t_field, x_field, y_field
#=== Plotting functions ===
def add_inner_title(ax, title, loc, prop=None, white_border=True, **kwargs):
"""Create a fancy inner title for plots
Parameters
==========
ax: `matplotlib.axex`
The axis which the inner title will be created
title: str
The title created
loc: int
Location. Can be verbose as well I think...
prop: dict, optional
A dictionary containing the properties of the title
eg. dict(size=18, color='green')
white_border: bool, optional
If True, no white borders are drawn around the text
Return
======
at: artist
The inner title artist that is added to the plot
"""
if prop is None:
prop = dict(size=plt.rcParams['legend.fontsize'])
at = AnchoredText(title, loc=loc, prop=prop,
pad=0., borderpad=0.75,
frameon=False, **kwargs)
ax.add_artist(at)
if white_border:
at.txt._text.set_path_effects([withStroke(foreground="w", linewidth=3)])
return at
def plot_beam(x,y,z,ptitle=None):
"""Create a simple beam plot
Parameters
==========
x, y, z: <numpy.ndarray>
Output positional and beam value vectors from `compute_beam_image`
ptitle: str
The title of the plot
Returns
=======
Create an image
"""
fig = plt.figure(1, figsize=(12,12))
if ptitle != None:
fig.suptitle('Primary beam model: {0:s}'.format(ptitle), fontsize = 18)
plt.gca().set_aspect('equal')
beam = plt.pcolor(x, y, z, rasterized = True, shading='auto')
cb = plt.colorbar(beam, aspect=30, fraction=0.04975, pad=0)
plt.rcParams['contour.negative_linestyle'] = 'solid'
cb.ax.yaxis.get_offset_text().set_fontsize(18)
cb.ax.tick_params(labelsize=18)
cb.ax.tick_params(direction='in', length=6, width=2)
cb.ax.set_ylabel(r'B', fontsize = 18)
CS = plt.contour(x, y, z,
levels = [np.amax(z)/2,np.amax(z)*2/3],
colors='w', alpha=1.)
plt.clabel(CS, fontsize=14, inline=True)
plt.xlabel(r'x', fontsize = 18)
plt.ylabel(r'y', fontsize = 18)
plt.show()
plt.close()
def plot_fractional_PB_change(dPB_array, dx_array, PB_err,
fname,
ptitle=None,
PB_FWHM = 60,
PB_normalisation = True,
cont_levels = [],
half_beam = True):
"""Generates fig 2 plot from Mooley et al. 2018
Currently only working with the RA smeared beam model
TO DO: make this function more general to work with an arbritary PB model
for more info see `fractional_PB_cange_of_single_slice()`
Parameters
==========
dPB_array: <numpy.ndarray>
The pixel slice coordinates as a grid (N_image x N_dx matrix with
constant columns)
dx_array: <numpy.ndarray>
The generated dx value coordinates as a grid (N_image x N_dx matrix with
constant rows)
PB_err: <numpy.ndarray>
The fractional bema change values as a grid (N_image x N_dx matrix)
fname: str
The full path and name and extension of the output image
ptitle: str, opt
The title of the plot
PB_FWHM: float, opt
The Primary Beam FWHM in arcseconds. Used to calculate the std of the PB
PB_normalisation: bool, opt
If True the plots x and y axis will be normalised by the PB FWHM
cont_levels: list, optional
The contour level values (default values if not given)
half_beam: bool, optional
If True, only half of the beam is being polotted due to symmetry reasons
Returns
=======
Create an image
"""
fig = plt.figure(1, figsize=(9,4.5))
if ptitle != None:
plt.title(ptitle, fontsize = 18)
if PB_normalisation:
beam = plt.pcolor(dPB_array / PB_FWHM, dx_array / PB_FWHM, PB_err,
rasterized=True, shading='auto')
else:
beam = plt.pcolor(dPB_array, dx_array, PB_err,
rasterized=True, shading='auto')
cb = plt.colorbar(beam, aspect=30, fraction=0.04975, pad=0)
cb.ax.yaxis.get_offset_text().set_fontsize(18)
cb.ax.tick_params(labelsize=18)
cb.ax.tick_params(direction='in', length=6, width=2)
cb.ax.set_ylabel(r'B$_{eff, \nu_0}$/B$_{0, \nu_0}$', fontsize = 18)
#Contours
if len(cont_levels) == 0:
CS = plt.contour(dPB_array / PB_FWHM, dx_array / PB_FWHM, PB_err,
colors='w', alpha=1.)
else:
CS = plt.contour(dPB_array / PB_FWHM, dx_array / PB_FWHM, PB_err,
levels = cont_levels,
colors='w', alpha=1.)
plt.clabel(CS, fontsize=14, inline=True)
if PB_normalisation:
plt.xlabel(r's$_0$/$\Theta_{B_{0, \nu_0}}$', fontsize = 18)
plt.ylabel(r'$\Delta$s/$\Theta_{B_{0, \nu_0}}$', fontsize = 18)
else:
plt.xlabel(r's$_0$', fontsize = 18)
plt.ylabel(r'$\Delta$s', fontsize = 18)
if half_beam:
plt.xlim(0)
#plt.show()
plt.savefig(fname, bbox_inches='tight')
plt.close()
def plot_fractional_PB_change_with_frequency(dPB_array, nu_array, PB_err,
fname,
ptitle=None,
PB_FWHM = None,
PB_normalisation = True,
cont_levels = [],
half_beam = True,
inner_title = None):
"""Generates fig 2 plot from Mooley et al. 2018
Currently only working with the RA smeared beam model
TO DO: make this function more general to work with an arbritary PB model
for more info see `fractional_PB_cange_of_single_slice()`
Parameters
==========
dPB_array: <numpy.ndarray>
The pixel slice coordinates as a grid (N_image x N_nu matrix with
constant columns)
nu_array: <numpy.ndarray>
The generated frequency value coordinates as a grid (N_image x N_nu
matrix with constant rows) given in GHz
PB_err: <numpy.ndarray>
The fractional bema change values as a grid (N_image x N_nu matrix)
fname: str
The full path and name and extension of the output image
ptitle: str, opt
The title of the plot
PB_FWHM_reference: float, opt
The Primary Beam FWHM for a given ref frequency in arcseconds.
PB_normalisation: bool, opt
If True the plots x and y axis will be normalised by the PB FWHM
cont_levels: list, optional
The contour level values (default values if not given)
half_beam: bool, optional
If True, only half of the beam is being polotted due to symmetry reasons
Returns
=======
Create an image
"""
#get the PB at the lowest frequency
if PB_FWHM == None:
PB_FWHM = get_VLA_PB_FWHM_from_freq(np.amin(nu_array))
fig = plt.figure(1, figsize=(9,4.5))
if ptitle != None:
plt.title(ptitle, fontsize = 18)
if PB_normalisation:
beam = plt.pcolor(dPB_array / PB_FWHM, nu_array, PB_err,
rasterized=True, shading='auto')
else:
beam = plt.pcolor(dPB_array, nu_array, PB_err,
rasterized=True, shading='auto')
cb = plt.colorbar(beam, aspect=30, fraction=0.04975, pad=0)
cb.ax.yaxis.get_offset_text().set_fontsize(18)
cb.ax.tick_params(labelsize=18)
cb.ax.tick_params(direction='in', length=6, width=2)
cb.ax.set_ylabel(r'B$_{eff}$/B$_0$', fontsize = 18)
#Contours
if len(cont_levels) == 0:
CS = plt.contour(dPB_array / PB_FWHM, nu_array, PB_err,
colors='w', alpha=1.)
else:
CS = plt.contour(dPB_array / PB_FWHM, nu_array, PB_err,
levels = cont_levels,
colors='w', alpha=1.)
plt.clabel(CS, fontsize=14, inline=True)
if PB_normalisation:
plt.xlabel(r's$_0$/$\Theta_{B,\nu_{0}}$', fontsize = 18)
else:
plt.xlabel(r's$_0$', fontsize = 18)
plt.ylabel(r'$\nu$ [GHz]', fontsize = 18)
if half_beam:
plt.xlim(0)
if inner_title != None:
ax = plt.gca()
add_inner_title(ax, inner_title, loc=3, prop=dict(size=16))
#plt.show()
plt.savefig(fname, bbox_inches='tight')
plt.close()
def plot_field_scan(RA_field,
Dec_field,
fname,
ptitle = '',
setting = False,
RA_field_setting = None,
Dec_field_setting = None):
"""Plot a scan pattern of a field, both the rising and setting scans if the
latter is provided.
Basicaly creates similar figures as Fig. 3 in Wang et aql. 2021
Parameters
==========
RA_field: <numpy.ndarray>
The RA values of the full field scans in degrees
Dec_field: <numpy.ndarray>
The Dec values of the full field scans in degrees
fname: str
The full path and name and extension of the output image
ptitle: str, opt
The title of the plot
setting: bool, optional
If True, the setting scanning pattern is also drawn
RA_field_setting: <numpy.ndarray>, optional
The RA values of the full field setting scans in degrees
Dec_field_setting: <numpy.ndarray>, optional
The Dec values of the full field setting scans in degrees
Returns
=======
Create an image
"""
if setting:
if np.shape(RA_field) != np.shape(RA_field_setting) or \
np.shape(Dec_field) != np.shape(Dec_field_setting):
raise ValueError("Input scan coordinate lengths are different!")
fig = plt.figure(1, figsize=(9,4.5))
if ptitle != None:
plt.title(ptitle, fontsize = 18)
if setting:
plt.plot(RA_field, Dec_field, '-', c=c0, lw=2., alpha=1., label='rising')
plt.plot(RA_field_setting, Dec_field_setting, '-', c=c2, lw=2., alpha=1.,
label='setting')
legend = plt.legend(loc='upper right', fontsize=18, framealpha=1., edgecolor='black')
legend.get_frame().set_linewidth(2.)
else:
plt.plot(RA_field, Dec_field, '-', c=c0, lw=2., alpha=1.)
ax = plt.gca()
ax.invert_xaxis()
plt.ylabel(r'Dec (J2000) [deg]', fontsize = 18)
plt.xlabel(r'RA (J2000) [deg]', fontsize = 18)
plt.savefig(fname, bbox_inches='tight')
plt.close()
#=== FUNCTIONS GENERATING THE ANALYTIC PLOTS ===
def VLA_OTFM_plots():
"""Re-create the Mooley et al. 2018 Figs 2. & 3.
More precisely, fig 3. is not re-created but a similar figure with more
physical mening/relevance is shown. Where the x axis uses a fixed PB
FWHM (the smallest one), as the pointings are frequency dependent.
This is basically a check that our code works as intended
"""
telescope = 'VLA'
pixelsize = 2.5 #This equals the synthesised beam size for now
dx = 4 * 60 #The absolute RA offset within each subscan
nu_central = 3.
N_chan = 201
chan_width = 0.01
#PB_FWHM = 60 * 10 #The PB FWHM in arcseconds at 3 GHz
PB_FWHM = get_VLA_PB_FWHM_from_freq(nu_central)
#Get PB values
nu_array = get_pixel_coordinate_array(N_chan,chan_width) + nu_central
#=== Fractional beam change ===
contours = [0.825,0.85,0.875,0.9,0.925,0.95,0.97,0.98,0.99,0.995,1.,\
1.005,1.01,1.02,1.03,1.05,1.1,1.2,1.3]
#RA
dPB_array, dx_array, PB_err = fractional_PB_cange_of_single_slice(
pixelsize = pixelsize,
PB_FWHM = PB_FWHM,
Ndx = 99, #So it goes up to dx=0.1
fdx = 0.01,
slice_dir='x',
slice_at=0.,
N_image=((PB_FWHM * 1.6) / pixelsize) + 1) #So I plot up to the +/- 0.8
#of the PB FWHM
parameters_string = r'{0:s} OTFM at {1:0.1f} GHz with' \
'$\Theta_B$={2:.1f}\' & $\Theta_s$={3:.1f}\" RA slice'.format(
telescope,nu_central,PB_FWHM / 60,pixelsize)
fname = './figures/{0:s}_fractional_beam_change_RA_slice.pdf'.format(
telescope)
plot_fractional_PB_change(dPB_array, dx_array, PB_err,
fname = fname,
ptitle = parameters_string,
PB_FWHM = PB_FWHM,
cont_levels=contours)
#=== Frequency dependence ===
contours = [0.97,0.98,0.99,0.995,1.,1.005,1.01,1.02,1.03,1.04,1.06]
#Normalising with max freq
nu_ext = np.amax(nu_array)
nu_ext_PB_FWHM = get_VLA_PB_FWHM_from_freq(nu_ext)
dPB_array, dnu_array, PB_err = fractional_PB_cange_with_frequency_of_single_slice(
pixelsize = pixelsize,
N_image = ((nu_ext_PB_FWHM * 1.6) / pixelsize) + 1,
nu_central = nu_central,
Nnu = N_chan,
dnu = chan_width,
dx = dx,
PB_scaling_model = telescope,
slice_dir='x')
parameters_string = r'{0:s} OTFM with ' \
r'$\Delta$x={1:.1f}'"'"r' & $\nu$$_0$={2:.1f} GHz RA slice'.format(
telescope, dx/60, nu_ext)
fname = './figures/\
freq_dependent_{0:s}_fractional_beam_change.pdf'.format(telescope)
plot_fractional_PB_change_with_frequency(dPB_array, dnu_array, PB_err,
fname = fname,
ptitle = parameters_string,
PB_FWHM = nu_ext_PB_FWHM,
cont_levels = contours)
def plot_example_MeerKLASS_field_scans():
"""Create a similar figure for the simulated scans as in Wang et al. 2021
Fig. 3.
"""
t_field_rising, x_field_rising, y_field_rising = \
get_field_scan_coordinate_values_constant_az_slew()
t_field_setting, x_field_setting, y_field_setting = \
get_field_scan_coordinate_values_constant_az_slew(
t0 = Time('2019-02-25T00:40:11'),
alt0 = 41.,
az0=-43.)
plot_field_scan(RA_field = x_field_rising,
Dec_field = y_field_rising,
fname = './figures/example_scan_coverage.pdf',
ptitle = '',
setting = True,
RA_field_setting = x_field_setting,
Dec_field_setting = y_field_setting)
def MeerKLASS_OTFM_plots(band='L'):
"""Create the OTFM smearing plots for MeerKLASS
"""
telescope = 'MeerKAT'
#From the MeerKLass setup i.e Wang et al. 2021
if band == 'L':
pixelsize = L_band_PIXELSIZE
nu_central = L_BAND_NU_CENTRAL
N_chan = N_CHAN
chan_width = L_BAND_CW
elif band == 'UHF':
pixelsize = UHF_band_PIXELSIZE
nu_central = UHF_BAND_NU_CENTRAL
N_chan = N_CHAN
chan_width = UHF_BAND_CW
else:
raise ValueError('Invalid band name (S is not supported)!')
#I use the top frequency of the band as discussed in the memo
#PB_FWHM = get_MeerKAT_PB_FWHM_from_freq(nu_central)
#Get PB values
nu_array = get_pixel_coordinate_array(N_chan,chan_width) + nu_central
PB_FWHM = get_MeerKAT_PB_FWHM_from_freq(np.amax(nu_array))
#=== Fractional beam change ===
contours = [0.85,0.9,0.95,0.99,1.,1.01,1.05,1.1,1.2,1.3,1.4,1.5]
dPB_array, dx_array, PB_err = fractional_PB_cange_of_single_slice(
pixelsize = pixelsize,
PB_FWHM = PB_FWHM,
Ndx = 99, #So it goes up to dx=0.1
fdx = 0.01,
slice_dir='x',
slice_at=0.,
N_image=((PB_FWHM * 1.6) / pixelsize) + 1) #So I plot up to the +/- 0.8
#of the PB FWHM
#parameters_string = r'Fractional PB cange for {0:s} OTFM at {1:0.1f} GHz with' \
# '$\Theta_B$={2:.1f}\' & $\Theta_s$={3:.1f}\" x slice'.format(
# telescope,nu_central,PB_FWHM / 60,pixelsize)
fname = './figures/{0:s}_fractional_beam_change_{1:s}.pdf'.format(
telescope, band)
plot_fractional_PB_change(dPB_array, dx_array, PB_err,
fname = fname,
ptitle = None,
PB_FWHM = PB_FWHM,
cont_levels=contours)
#=== Frequency dependence ===
contours = [0.9,0.925,0.95,0.975,0.99,1.,1.01,1.05,1.1,1.2]
#Normalising with max freq
nu_ext = np.amax(nu_array)
nu_ext_PB_FWHM = get_MeerKAT_PB_FWHM_from_freq(nu_ext)
for dt in [2.,4.,6.,8.,10.,12.]:
#Get ds for the frequncy dependence plots
#Rising field parameters
dRA, dDec = get_smearing_from_altaz_subscan(t0 = Time('2019-02-24T19:48:28'),
alt0 = 42.,
az0 = 59.6,
slew_velocity = 5.,
slew_direction = -1,
dt = dt)
ds = np.sqrt(np.power(dRA,2) + np.power(dDec,2))
dPB_array, dnu_array, PB_err = fractional_PB_cange_with_frequency_of_single_slice(
pixelsize = pixelsize,
N_image = ((nu_ext_PB_FWHM * 1.6) / pixelsize) + 1,
nu_central = nu_central,
Nnu = N_chan,
dnu = chan_width,
dx = ds,
PB_scaling_model = telescope,
slice_dir='x')
#parameters_string = r'{0:s} OTFM with ' \
#r'$\Delta$x={1:.1f}'"'"r' & $\nu$$_0$={2:.1f} GHz RA slice'.format(
# telescope, dx/60, nu_ext)
fname = './figures/\
{0:s}_freq_dependent_fractional_beam_change_ds_{1:d}_{2:s}.pdf'.format(telescope,
int(dt), band)
plot_fractional_PB_change_with_frequency(dPB_array, dnu_array, PB_err,
fname = fname,
ptitle = '',
PB_FWHM = nu_ext_PB_FWHM,
cont_levels = contours,
inner_title = \
r"$\Delta$T = {0:d}s & $\Delta$s = {1:.2f}".format(
int(dt),ds / nu_ext_PB_FWHM) \
+ r" $\Theta_{B,\nu_{0}}$")
def print_PB_parameters_and_smearing_difference(band='L'):
"""Just print out various useful parameters for the discussion
"""
telescope = 'MeerKAT'
#From the MeerKLass setup i.e Wang et al. 2021
if band == 'L':
pixelsize = L_band_PIXELSIZE
#nu_central = 1.284
nu_central = L_BAND_NU_CENTRAL
N_chan = N_CHAN
chan_width = L_BAND_CW
elif band == 'UHF':
pixelsize = UHF_band_PIXELSIZE
nu_central = UHF_BAND_NU_CENTRAL
N_chan = N_CHAN
chan_width = UHF_BAND_CW
else:
raise ValueError('Invalid band name (S is not supported)!')
print('Selected observing band: {0:s}'.format(band))
#Get PB values
nu_array = get_pixel_coordinate_array(N_chan,chan_width) + nu_central
print('The PB FWHM at min freq: {0:.2f} arcminutes'.format(
arcsecond_to_arcmin(get_MeerKAT_PB_FWHM_from_freq(np.amin(nu_array)))))
print('The PB FWHM at max freq: {0:.2f} arcminutes'.format(
arcsecond_to_arcmin(get_MeerKAT_PB_FWHM_from_freq(np.amax(nu_array)))))
#Compute the difference for the rising and setting skyes
for dt in [2.,4.,6.,8.,10.,12.]:
#Get ds for the frequncy dependence plots
#Rising field parameters
dRA, dDec = get_smearing_from_altaz_subscan(
t0 = Time('2019-02-24T19:48:28'),
alt0 = 42.,
az0 = 59.6,
slew_velocity = 5.,
slew_direction = -1,
dt = dt)
ds = np.sqrt(np.power(dRA,2) + np.power(dDec,2))
#Setting field parameters
dRA_setting, dDec_setting = get_smearing_from_altaz_subscan(
t0 = Time('2019-02-25T00:40:11'),
alt0 = 41.,
az0=-43.,
dt = dt)
ds_setting = np.sqrt(np.power(dRA_setting,2) + np.power(dDec_setting,2))
print('For dT={0:d}s the rising and setting smearings are {1:.2f} and \
{2:.2f} arcminutes, with a difference of {3:.2f} arcminutes'.format(
int(dt),
arcsecond_to_arcmin(ds),
arcsecond_to_arcmin(ds_setting),
arcsecond_to_arcmin(np.fabs(np.subtract(ds,ds_setting)))))
#=== MAIN ===
if __name__ == "__main__":
# Need to have a figures/ subdir
#VLA_OTFM_plots()
#plot_example_MeerKLASS_field_scans()
#MeerKLASS_OTFM_plots(band='L')
#print_PB_parameters_and_smearing_difference(band='L')
MeerKLASS_OTFM_plots(band='UHF')
print_PB_parameters_and_smearing_difference(band='UHF')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment