Hacky coverage analysis for GBT/Argus

sparse_argus_map.py

#!/usr/bin/env python3

import itertools

import numpy as np

from matplotlib import patches
from matplotlib import pyplot as plt

from astropy import units as u
from astropy import (coordinates, convolution)


plt.rc('font', size=10, family='serif')
plt.rc('text', usetex=True)


cell_size   =  1  # arcsec
array_size  = 90
edge_size   = 90
pix_spacing = 30
#row_delta   = 22  # best for no subrow separation
#row_delta   = 16  # best for subrow sampling
row_delta   = 20  # 27 at -3deg, 5 at 30deg for 26.4 dither
field_size  =  2 * array_size + row_delta
pangles = np.linspace(0, 35, 36)

img_size = 2 * edge_size + 2 * array_size
n_pix = int(img_size / cell_size)
img = np.zeros((n_pix, n_pix))

i_x1 = i_y1 = int((edge_size + array_size / 2) / cell_size - row_delta / cell_size / 2)
i_x2 = i_y2 = i_x1 + int(array_size / cell_size + row_delta / cell_size)

argus_x_offsets = np.array([
    [-45, -15, 15, 45],
    [-45, -15, 15, 45],
    [-45, -15, 15, 45],
    [-45, -15, 15, 45],
])
argus_y_offsets = argus_x_offsets.T[::-1,:]
argus_offsets = np.array([argus_x_offsets.flatten(), argus_y_offsets.flatten()]).T


n_covg = 6
#subrow_offsets = itertools.cycle(np.array([-30, -15, 0, 15, 30]))
subrow_offsets_arr = np.linspace(-13.3, 13.3, n_covg).round().astype(int)
subrow_offsets = itertools.cycle(subrow_offsets_arr)
#test_pa_set = np.random.uniform(14, 40, size=2*n_covg).reshape(-1, 2) * u.deg
#test_pa_set = np.array([[17, 17], [23, 23], [30, 30], [34, 34], [36, 36]]) * u.deg
#test_pa_set = np.linspace(14, 40, 10).reshape(-1, 2) * u.deg
test_pa_start = 19
test_pa_set = np.linspace(test_pa_start, test_pa_start+6, 2*n_covg).reshape(-1, 2) * u.deg

for dix_subrow, (x_pa, y_pa) in zip(subrow_offsets, test_pa_set):
    for offset in argus_offsets:
        x_off, y_off = offset
        radius = np.sqrt(x_off**2 + y_off**2)
        theta = np.arctan(y_off / x_off) * u.rad
        i_xo_t = int(np.sign(x_off) * radius * np.cos(theta + x_pa) / cell_size)
        i_yo_t = int(np.sign(x_off) * radius * np.sin(theta + y_pa) / cell_size)
        # "ra scan"
        img[i_xo_t+i_x1+dix_subrow, :] += 1
        img[i_xo_t+i_x2+dix_subrow, :] += 1
        # "dec scan"
        img[:, i_yo_t+i_y1+dix_subrow] += 1
        img[:, i_yo_t+i_y2+dix_subrow] += 1

# convolution
hpbw = 8.0 / cell_size
std_pix = hpbw / (2 * np.sqrt(2 * np.log(2)))
kernel = convolution.Gaussian2DKernel(x_stddev=std_pix)
c_img = convolution.convolve(img, kernel, boundary='extend')


# compute image the integration time statistics over the field
rect_width = int(field_size / cell_size)
rect_corner = int(n_pix / 2 - rect_width / 2)
sub_img = c_img[rect_corner:rect_corner+rect_width+1,
                rect_corner:rect_corner+rect_width+1]
quants = np.quantile(sub_img.flatten(),
        [0.01, 0.10, 0.50, 0.90, 0.99])
q01, q10, q50, q90, q99 = quants
img_mad = np.median(np.abs(sub_img.flatten()/q50-1))
img_std = np.std(sub_img.flatten()/sub_img.mean()-1)
print(f'-- MAD       -> {img_mad}')
print(f'-- STD       -> {img_std}')
print(f'-- quantiles -> {quants}')
print(f'   med-ratio -> {q10/q50:2.3f}, {q90/q50:2.3f}')


fig, ax = plt.subplots(figsize=(4, 3.5))
imsh = ax.imshow(np.sqrt(c_img/q50), vmin=0.0, vmax=1.1, origin='lower', cmap='magma')
cbar = fig.colorbar(imsh, ax=ax, extend='max', pad=0.02, shrink=0.96)
cbar.minorticks_on()
mcol = '0.6'
ax.plot([i_x1, i_x1, i_x2, i_x2], [i_y1, i_y2, i_y1, i_y2], color=mcol,
        linestyle='none', marker='.', markersize=4)
#ax.plot([i_x1-30, i_x1-15, i_x1+15, i_x1+30, i_y1],
#        [i_y1, i_y1, i_y1, i_y1, i_y1], color=mcol,
#        linestyle='none', marker='|', markersize=4)
#ax.plot([i_x1, i_x1, i_x1, i_x1, i_x1],
#        [i_y1-30, i_y1-15, i_y1+15, i_y1+30, i_y1], color=mcol,
#        linestyle='none', marker='_', markersize=4)
ax.plot(len(subrow_offsets_arr)*[i_x1],
        i_y1+subrow_offsets_arr/cell_size, color='magenta',
        linewidth=0.8, linestyle='none', marker='_', markersize=4)
rect = patches.Rectangle((rect_corner, rect_corner), rect_width, rect_width,
        linewidth=0.8, linestyle='dashed', edgecolor=mcol, facecolor='none')
ax.add_patch(rect)
for offset in argus_offsets:
    b_pa = test_pa_start * u.deg
    x_off, y_off = offset
    radius = np.sqrt(x_off**2 + y_off**2)
    theta = np.arctan(y_off / x_off) * u.rad
    i_xo = int(np.sign(x_off) * radius * np.cos(theta + b_pa) / cell_size)
    i_yo = int(np.sign(x_off) * radius * np.sin(theta + b_pa) / cell_size)
    beam_pos = (i_x1 + i_xo, i_y1 + i_yo)
    beam = patches.Circle(beam_pos, radius=hpbw/2, edgecolor='black',
            facecolor='none', linewidth=0.7)
    ax.add_patch(beam)
ax.set_xticks(np.array([-120,-60,0,60,120])+n_pix/2)
ax.set_xticklabels([-2, -1, 0, 1, 2])
ax.set_yticks(np.array([-120,-60,0,60,120])+n_pix/2)
ax.set_yticklabels([-2, -1, 0, 1, 2])
ax.set_xlabel('RA Offset [arcmin]')
ax.set_ylabel('Dec Offset [arcmin]')
plt.tight_layout()
plt.savefig('test_stripes_dsub_drow.pdf', dpi=300)
plt.close('all')
del fig
del ax