In [ ]:

    

import numpy as np
import os

import astropy.units as u
import astropy.constants as const
import fastimgproto.imager as imager
import fastimgproto.visibility as visibility
import fastimgproto.gridder.conv_funcs as kfuncs

from astropy.coordinates import Angle, SkyCoord, AltAz, EarthLocation
from astropy.time import Time

from fastimgproto.skymodel.helpers import SkyRegion, SkySource
from fastimgproto.sourcefind.image import SourceFindImage
from fastimgproto.telescope.readymade import Meerkat
from fastimgproto.bindings import cpp_image_visibilities, CppKernelFuncs


    

In [ ]:

    

import matplotlib.pyplot as plt
import matplotlib.image as mpimg
%matplotlib inline
# Plot image pixels in cartesian ordering (i.e. y-positive == upwards):
plt.rcParams['image.origin'] = 'lower'
# Make plots bigger
plt.rcParams['figure.figsize'] = 10, 10


    

In [ ]:

    

telescope = Meerkat()
print("Telescope with {} antennae == {} baselines".format(
    len(telescope.ant_local_xyz), len(telescope.baseline_local_xyz)))
print("Centre: {!r}, {!r}".format(telescope.lon, telescope.lat))


    

In [ ]:

    

pointing_centre = SkyCoord(0 * u.deg, -30 * u.deg)
obs_central_frequency = 3. * u.GHz
wavelength = const.c / obs_central_frequency
transit_time = telescope.next_transit(pointing_centre.ra,
                                      start_time=Time('2017-01-01'))


    

In [ ]:

    

altaz = pointing_centre.transform_to(
    AltAz(obstime=transit_time,
         location=telescope.centre))
altaz.alt.deg


    

In [ ]:

    

nstep=10
obs_times = transit_time + np.linspace(-1, 1, nstep) * u.hr
print("Generating UVW-baselines for {} timesteps".format(nstep))
uvw_m = telescope.uvw_tracking_skycoord(pointing_centre, obs_times)
# From here on we use UVW as multiples of wavelength, lambda:
uvw_lambda = (uvw_m / wavelength).to(u.dimensionless_unscaled).value



# Additional source to North-East of pointing centre
extra_src_position = SkyCoord(ra=pointing_centre.ra + 0.01 * u.deg,
                              dec=pointing_centre.dec + 0.01 * u.deg, )

steady_sources = [
    SkySource(pointing_centre, flux=1 * u.Jy),
    SkySource(extra_src_position, flux=0.4 * u.Jy),
]

# Simulate incoming data; includes transient sources, noise:
print("Simulating visibilities")
data_vis = visibility.visibilities_for_source_list(
    pointing_centre,
    source_list = steady_sources, 
    uvw = uvw_lambda)

vis_noise_level = 0.1 * u.Jy
data_vis = visibility.add_gaussian_noise(vis_noise_level, data_vis)


    

In [ ]:

    

image_size=1024 * u.pixel
cell_size=1 * u.arcsecond


    

In [ ]:

    

kernel_support = 3
trunc = kernel_support
kernel_func = kfuncs.GaussianSinc(trunc=trunc)
image, beam = imager.image_visibilities(
    data_vis, 
    uvw_lambda,
    image_size=image_size,
    cell_size=cell_size,
    kernel_func=kernel_func,
    kernel_support=kernel_support,
    kernel_exact=True,
    kernel_oversampling=None
    )
re_image = np.real(image)


    

In [ ]:

    

import stp_python
image_size = image_size.to(u.pix)
# Size of a UV-grid pixel, in multiples of wavelength (lambda):
grid_pixel_width_lambda = 1.0 / (cell_size.to(u.rad) * image_size)
uvw_in_pixels = (uvw_lambda / grid_pixel_width_lambda).value


    

In [ ]:

    

bind_image, bind_beam = cpp_image_visibilities(
    data_vis, 
    uvw_lambda,
    image_size=image_size,
    cell_size=cell_size,
    kernel_func_name=CppKernelFuncs.gauss_sinc,
    kernel_trunc_radius = trunc,
    kernel_support=kernel_support,
    kernel_exact=True,
    kernel_oversampling=None
    )


    

In [ ]:

    

bind_diffmax = np.max(np.abs(image-bind_image))
print(bind_diffmax)
bind_re_diffmax = np.max(np.abs(np.real(image-bind_image)))
print(bind_re_diffmax)
# assert cpp_diffmax < 2e-15


    

In [ ]:

    

vis_filepath = 'viz.npz'
with open(vis_filepath, 'wb') as f:
    np.savez(f,
         uvw_lambda=uvw_lambda,
         model=data_vis,
         vis=data_vis,
         )


    

In [ ]:

    

from fastimgproto.scripts.config import make_config
cpp_config_filepath = 'cpp_config.json'
with open(cpp_config_filepath, 'w') as f:
    make_config(f)


    

In [ ]:

    

import subprocess
cpp_image_filepath = 'cpp_img.npz'
result = subprocess.check_call([
    'fastimg_cpp_imagevis',
     cpp_config_filepath,
     vis_filepath,
     cpp_image_filepath
    ]
)
result


    

In [ ]:

    

with open(cpp_image_filepath, 'rb') as f:
    cpp_imgdata = np.load(f)
    cpp_image = cpp_imgdata['image']


    

In [ ]:

    

cpp_image = np.real(cpp_image)


    

In [ ]:

    

cpp_diffmax = np.max(np.abs(image-cpp_image))
print(cpp_diffmax)
cpp_re_diffmax = np.max(np.abs(np.real(image-cpp_image)))
print(cpp_re_diffmax)
# assert cpp_diffmax < 2e-15


    

In [ ]:

    

fig, ax = plt.subplots(1,2)
axlims = 400, 600
def plot_on_ax(ax,img):
    im = ax.imshow(np.real(img))
    ax.set_xlim(axlims)
    ax.set_ylim(axlims)
#     fig.colorbar(im)
plot_on_ax(ax[0], image)
plot_on_ax(ax[1], bind_image)


    

In [ ]:

    

fig, ax1 = plt.subplots(1,1)
im = ax1.imshow(np.real(bind_image-image))
axlims = 400, 600
ax1.set_xlim(axlims)
ax1.set_ylim(axlims)
fig.colorbar(im)


    

In [ ]: