In [ ]:
%matplotlib inline
import sys
sys.path.append('../..')
from matplotlib import pylab
from ipywidgets import interact
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
pylab.rcParams['figure.figsize'] = 12, 10
import functools
import numpy
import scipy
import scipy.special
from astropy.coordinates import SkyCoord
from astropy import units
from crocodile.clean import *
from crocodile.synthesis import *
from crocodile.simulate import *
from util.visualize import *
from arl.test_support import create_named_configuration
Generate baseline coordinates for an observation with the VLA over 10 minutes, with a visibility recorded every 10 seconds. The phase center is fixed at a declination of 30 degrees.
In [ ]:
vlas = create_named_configuration('VLAA')
ha0 = 10
tsnap = 30*60
tdump = 20
ha_range = numpy.arange(numpy.radians(ha0),
numpy.radians(ha0 + 360 * tsnap / 3600 / 24),
numpy.radians(360 * tdump / 3600 / 24))
dec = numpy.radians(30)
vobs = xyz_to_baselines(vlas.data['xyz'], ha_range, dec)
# Wavelength: 5 metres
wvl=5
uvw = vobs / wvl
ax = plt.figure().add_subplot(121, projection='3d')
ax.scatter(uvw[:,0], uvw[:,1] , uvw[:,2])
max_uvw = numpy.amax(uvw)
ax.set_xlabel('U [$\lambda$]'); ax.set_xlim((-max_uvw, max_uvw))
ax.set_ylabel('V [$\lambda$]'); ax.set_ylim((-max_uvw, max_uvw))
ax.set_zlabel('W [$\lambda$]'); ax.set_zlim((-max_uvw, max_uvw))
ax.view_init(0, 20)
pylab.show()
In [ ]:
import itertools
vis = numpy.zeros(len(uvw), dtype=complex)
for u,v in itertools.product(range(-3, 4), range(-3, 4)):
vis += 1.0*simulate_point(uvw, 0.010*u, 0.010*v)
plt.clf()
uvdist=numpy.sqrt(uvw[:,0]**2+uvw[:,1]**2)
plt.plot(uvdist, numpy.abs(vis), '.', color='r')
In [ ]:
mean_ha = numpy.mean(ha_range)
pc = SkyCoord(ra=-mean_ha, dec=dec, unit=units.rad)
new_pc = SkyCoord(ra=0, dec=vlas.location.latitude, unit=units.deg)
uvw_r, vis_r = visibility_rotate(uvw, vis, pc, new_pc)
print("Max w before rotation:", numpy.max(numpy.abs(uvw[:,2])))
print("Max w after rotation: ", numpy.max(numpy.abs(uvw_r[:,2])))
Now we get the old phase centre back into the centre of the picture by doing an FFT shift. We will have to replicate this shift in the kernels laters later on.
In [ ]:
l,m,n = skycoord_to_lmn(pc, new_pc)
print("Shifting image by (%f,%f)" % (-l,-m))
vis_s = visibility_shift(uvw_r, vis_r, -l,-m)
Even thought the image revolves around the old phase centre, it is still transformed. Reason is that when changing phase centre we change the projection plane. This can both cause skews as well as rotations.
We can determine the skew matrix relatively easily: The $(l,m,n)$ and $(u,v,w)$ coordinate system are parallel. This means we can easily learn the inverse transformation matrix by feeding the unit vectors to visibility_rotate. The upper 2x2 matrix is the approximate transformation matrix of the old to the new $(l,m)$ coordinate system.
Note that this is not perfect, as the new $(l,m)$ also depends a bit on $n$. This cannot be corrected without adjusting $w$, so we will have to deal with it using image reprojection after the fact.
In [ ]:
T3, _ = visibility_rotate(numpy.eye(3), [0,0,0], new_pc, pc)
T = T3[0:2,0:2]
Ti = numpy.linalg.inv(T)
print("Image transformation:\n%s\nDropped n transformation:\n%s" % (T, T3[0:2,2]))
uvw_t = uvw_transform(uvw_r, T)
uvw = uvw_t
vis = vis_s
We can visualise the new $u,v,w$ distribution. Note how every baseline is roughly centered around $w=0$ now.
In [ ]:
plt.scatter(uvw[:,0], uvw[:,1] , c=uvw[:,2], s=4,lw=0)
plt.scatter(-uvw[:,0], -uvw[:,1] , c=-uvw[:,2], s=4,lw=0)
max_uvw = numpy.amax(uvw)*1.1
plt.xlabel('U [$\lambda$]'); plt.xlim((-max_uvw, max_uvw))
plt.ylabel('V [$\lambda$]'); plt.ylim((-max_uvw, max_uvw))
plt.colorbar(shrink=.92);
In [ ]:
theta = 2*0.05
lam = 19000
wmax = numpy.max(numpy.abs(uvw[:,2]))
Nccvf = 2*theta*numpy.sqrt((wmax*theta/2)**2 + (wmax**1.5 * theta / 2 / numpy.pi / 0.01))
Naa = 30
NpixKern = int(numpy.ceil((numpy.sqrt(Naa**2 + Nccvf**2)-1) / 2)) * 2 + 1
print("Kernel size: %dx%d (%dx%d * %dx%d)" % (NpixKern, NpixKern, Nccvf,Nccvf, Naa,Naa))
Grid resolution in $uv$ is $1/\theta$, which is then oversampled. We choose the resolution in $w$ accordingly. This tells us how many kernels we need for imaging.
In [ ]:
Qpx=2
wstep=1/theta/Qpx
wcachesize=2*int(numpy.ceil(wmax/wstep))
print("Making w-kernel cache of %d kernels" % wcachesize)
Now we can generate our kernel cache and cosntruct the imaging worker. Note that we need to account for the transformations we did above: The kernel image needs to be shifted by $(l,m)$ and transformed with $T^{-1}$.
In [ ]:
import functools
wcache=pylru.FunctionCacheManager(functools.partial(w_kernel, T=numpy.linalg.inv(T)), wcachesize)
imgfn = functools.partial(w_cache_imaging, kernel_cache=wcache,
wstep=wstep, Qpx=Qpx, NpixFF=256, NpixKern=NpixKern, dl=l, dm=m)
After everything is set up, we can start imaging:
In [ ]:
d_w,p_w,_=do_imaging(theta, lam, uvw, None, vis, imgfn)
show_image(d_w, "image", theta)
In [ ]:
step=int(theta*lam/10)
def zoom_w(x=theta*lam/2,y=step): pylab.matshow(d_w[y:y+2*step,x:x+2*step]); pylab.colorbar(shrink=.4,pad=0.025); pylab.show()
interact(zoom_w, x=(0,d_w.shape[0]-2*step,step), y=(0,d_w.shape[1]-2*step,step))
We can have a look at the generated kernel cache to confirm what kernels were used:
In [ ]:
cached_gcfs = list(sorted(wcache.cache.items()))
def inspect_cache(i, ox=0,oy=0):
(pars, kwargs), gcf = cached_gcfs[i]
print("theta=%f, w=%f, %s" % (pars[0], pars[1], ", ".join(["%s=%f" % kv for kv in kwargs ]) ))
pylab.matshow(gcf[oy,ox].real); pylab.colorbar(shrink=.4,pad=0.025);
pylab.matshow(gcf[oy,ox].imag); pylab.colorbar(shrink=.4,pad=0.025);pylab.show()
interact(inspect_cache, i=(0,len(cached_gcfs)-1), ox=(0,Qpx-1), oy=(0,Qpx-1));
In [ ]: