Basic Simulation¶

[1]:

import os
try:
    import sirius
    print('SiRIUS version',sirius.__version__,'already installed.')
except ImportError as e:
    print(e)
    print('Installing SiRIUS')
    os.system("pip install sirius")
    import sirius
    print('SiRIUS version',sirius.__version__,' installed.')

SiRIUS version 0.0.28 already installed.

Load Packages¶

[2]:

import pkg_resources
import xarray as xr
import numpy as np
from astropy.coordinates import SkyCoord
xr.set_options(display_style="html")
import os
try:
    from google.colab import output
    output.enable_custom_widget_manager()
    IN_COLAB = True
    %matplotlib widget
except:
    IN_COLAB = False
    %matplotlib inline

Load Telescope Layout¶

[3]:

########## Get telescope layout ##########
tel_dir = pkg_resources.resource_filename('sirius_data', 'telescope_layout/data/vla.d.tel.zarr')
tel_xds = xr.open_zarr(tel_dir,consolidated=False)
n_ant = tel_xds.dims['ant_name']
#tel_xds.attrs['telescope_name'] = 'EVLA'
tel_xds

[3]:

<xarray.Dataset>
Dimensions:        (ant_name: 27, pos_coord: 3)
Coordinates:
  * ant_name       (ant_name) <U3 'W01' 'W02' 'W03' 'W04' ... 'N07' 'N08' 'N09'
  * pos_coord      (pos_coord) int64 0 1 2
Data variables:
    ANT_POS        (ant_name, pos_coord) float64 dask.array<chunksize=(27, 3), meta=np.ndarray>
    DISH_DIAMETER  (ant_name) float64 dask.array<chunksize=(27,), meta=np.ndarray>
Attributes:
    site_pos:        [{'m0': {'unit': 'm', 'value': -1601185.3650000016}, 'm1...
    telescope_name:  VLA

Create Time and Freq Xarrays¶

The chunking of time_xda and chan_xda determines the number of branches in the DAG (maximum parallelism = n_time_chunks x n_chan_chunks).

[4]:

from sirius.dio import make_time_xda
#time_xda = make_time_xda(time_start='2019-10-03T19:00:00.000',time_delta=3600,n_samples=10,n_chunks=4)
time_xda = make_time_xda(time_start='2019-10-03T19:00:00.000',time_delta=3600,n_samples=10,n_chunks=1)
time_xda

Number of chunks  1

[4]:

<xarray.DataArray 'array-51a16e79cab65fc8a0bc3efd397119f7' (time: 10)>
dask.array<array, shape=(10,), dtype=<U23, chunksize=(10,), chunktype=numpy.ndarray>
Dimensions without coordinates: time
Attributes:
    time_delta:  3600.0

[5]:

from sirius.dio import make_chan_xda
spw_name = 'SBand'
#chan_xda = make_chan_xda(freq_start = 3*10**9, freq_delta = 0.4*10**9, freq_resolution=0.01*10**9, n_channels=3, n_chunks=3)
chan_xda = make_chan_xda(freq_start = 3*10**9, freq_delta = 0.4*10**9, freq_resolution=0.01*10**9, n_channels=3, n_chunks=1)
chan_xda

Number of chunks  1

[5]:

<xarray.DataArray 'array-4bc1685db34e1d77709ed218dffb158a' (chan: 3)>
dask.array<array, shape=(3,), dtype=float64, chunksize=(3,), chunktype=numpy.ndarray>
Dimensions without coordinates: chan
Attributes:
    freq_resolution:  10000000.0
    spw_name:         sband
    freq_delta:       400000000.0

Beam Models¶

[6]:

from sirius_data.beam_1d_func_models.airy_disk import vla #Get airy disk parameters for VLA dish (obtained from https://open-bitbucket.nrao.edu/projects/CASA/repos/casa6/browse/casatools/src/code/synthesis/TransformMachines/PBMath.cc)
airy_disk_parms = vla
print(airy_disk_parms)

beam_models = [airy_disk_parms]
beam_model_map = np.zeros(n_ant,dtype=int) #Maps the antenna index to a model in beam_models.
beam_parms = {} #Use default beam parms.

#If no beam should be used:
#none_model = {'pb_func':'none','dish_diameter':0.0,'blockage_diameter':0.0}

'''
#If Zernike Polynomial should be used:
zpc_dir = pkg_resources.resource_filename('sirius_data', 'aperture_polynomial_coefficient_models/data/EVLA_avg_zcoeffs_SBand_lookup.apc.zarr')
zpc_xds = xr.open_zarr(zpc_dir,consolidated=False)

beam_models = [zpc_xds]
zpc_xds

beam_models = [zpc_xds,airy_disk_parms]

beam_model_map[0] = 1

zpc_xds
'''

{'func': 'casa_airy', 'dish_diam': 24.5, 'blockage_diam': 0.0, 'max_rad_1GHz': 0.014946999714079439}

[6]:

"\n#If Zernike Polynomial should be used:\nzpc_dir = pkg_resources.resource_filename('sirius_data', 'aperture_polynomial_coefficient_models/data/EVLA_avg_zcoeffs_SBand_lookup.apc.zarr')\nzpc_xds = xr.open_zarr(zpc_dir,consolidated=False)\n\nbeam_models = [zpc_xds]\nzpc_xds\n\nbeam_models = [zpc_xds,airy_disk_parms]\n\nbeam_model_map[0] = 1\n\nzpc_xds\n"

Polarization Setup¶

[7]:

#https://github.com/casacore/casacore/blob/dbf28794ef446bbf4e6150653dbe404379a3c429/measures/Measures/Stokes.h
# ['RR','RL','LR','LL'] => [5,6,7,8], ['XX','XY','YX','YY'] => [9,10,11,12]
pol = [5,8]

UVW Parameters¶

[8]:

# If using uvw_parms['calc_method'] = 'casa' .casarc must have directory of casadata.
import pkg_resources
casa_data_dir = pkg_resources.resource_filename('casadata', '__data__')
rc_file = open(os.path.expanduser("~/.casarc"), "a+")  # append mode
rc_file.write("\n measures.directory: " + casa_data_dir)
rc_file.close()

uvw_parms = {}
uvw_parms['calc_method'] = 'casa' #'astropy' or 'casa'
uvw_parms['auto_corr'] = False

Sources¶

point_source_ra_dec: [n_time, n_point_sources, 2] (singleton: n_time)¶

[9]:

point_source_skycoord = SkyCoord(ra='19h59m50.51793355s',dec='+40d48m11.3694551s',frame='fk5')
point_source_ra_dec = np.array([point_source_skycoord.ra.rad,point_source_skycoord.dec.rad])[None,None,:]

point_source_flux: [n_point_sources, n_time, n_chan, n_pol] (singleton: n_time, n_chan)¶

All 4 instramental pol values must be given even if only RR and LL are requested.

[10]:

point_source_flux = np.array([1.0, 0, 0, 1.0])[None,None,None,:]

Telescope Setup¶

phase_center: [n_time, 2] (singleton: n_time)¶

phase_center_names: [n_time] (singleton: n_time)¶

[11]:

phase_center = SkyCoord(ra='19h59m28.5s',dec='+40d44m01.5s',frame='fk5')
phase_center_ra_dec = np.array([phase_center.ra.rad,phase_center.dec.rad])[None,:]
phase_center_names = np.array(['field1'])

[12]:

%load_ext autoreload
%autoreload 2

[13]:

from sirius._sirius_utils._coord_transforms import _calc_rotation_mats, _directional_cosine,  _sin_project

print(phase_center_ra_dec[0,:],point_source_ra_dec[0,:,:])
lm_sin = _sin_project(phase_center_ra_dec[0,:],point_source_ra_dec[0,:,:])[0,:]
print('%.15f' % lm_sin[0],'%.15f' % lm_sin[1])

[5.23369701 0.71093805] [[5.2352982  0.71214946]]
0.001212034202758 0.001212034202764

pointing_ra_dec:[n_time, n_ant, 2] (singleton: n_time, n_ant) or None¶

[14]:

pointing_ra_dec = None #No pointing offsets

Noise¶

[15]:

noise_parms = None
#noise_parms ={}

Run Simulation¶

Write to disk as a measuremnt set and create DAGs.

[16]:

%load_ext autoreload
%autoreload 2

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload

[17]:

from sirius import simulation
save_parms = {'ms_name':'simple_sim.ms','write_to_ms':True,'DAG_name_vis_uvw_gen':'DAG_vis_uvw_gen.png','DAG_name_write':'DAG_write.png'}
ms_xds = simulation(point_source_flux,
                     point_source_ra_dec,
                     pointing_ra_dec,
                     phase_center_ra_dec,
                     phase_center_names,
                     beam_parms,beam_models,
                     beam_model_map,uvw_parms,
                     tel_xds, time_xda, chan_xda, pol, noise_parms, save_parms)
ms_xds

Setting default fov_scaling  to  4.0
Setting default mueller_selection  to  [ 0  5 10 15]
Setting default zernike_freq_interp  to  nearest
Setting default pa_radius  to  0.2
Setting default image_size  to  [1000 1000]
Setting default mode  to  dask_ms_and_sim_tool
10 351 3 2
<xarray.Dataset>
Dimensions:  (time: 10, pol: 2, chan: 3, baseline: 351, uvw: 3, time_chunk: 1,
              chan_chunk: 1, 4: 4)
Coordinates:
  * time     (time) <U23 '2019-10-03T19:00:00.000' ... '2019-10-04T04:00:00.000'
  * pol      (pol) int64 5 8
  * chan     (chan) float64 3e+09 3.4e+09 3.8e+09
Dimensions without coordinates: baseline, uvw, time_chunk, chan_chunk, 4
Data variables:
    DATA     (time, baseline, chan, pol) complex128 dask.array<chunksize=(10, 351, 3, 2), meta=np.ndarray>
    UVW      (time, baseline, uvw) complex128 dask.array<chunksize=(10, 351, 3), meta=np.ndarray>
    WEIGHT   (time, baseline, pol) float64 dask.array<chunksize=(10, 351, 2), meta=np.ndarray>
    SIGMA    (time, baseline, pol) float64 dask.array<chunksize=(10, 351, 2), meta=np.ndarray>
    TIMING   (time_chunk, chan_chunk, 4) float64 dask.array<chunksize=(1, 1, 4), meta=np.ndarray>
Meta data creation  0.4570884704589844
reshape time  0.0011935234069824219
*** Dask compute time 4.647282600402832

[17]:

<xarray.Dataset>
Dimensions:           (polarization_ids: 1, spw_ids: 1)
Coordinates:
  * polarization_ids  (polarization_ids) int64 0
  * spw_ids           (spw_ids) int64 0
Data variables:
    *empty*
Attributes:
    xds0:              <xarray.Dataset>\nDimensions:         (row: 3510, uvw_...
    SPECTRAL_WINDOW:   <xarray.Dataset>\nDimensions:          (row: 1, d1: 3)...
    POLARIZATION:      <xarray.Dataset>\nDimensions:       (row: 1, d1: 2, d2...
    DATA_DESCRIPTION:  <xarray.Dataset>\nDimensions:             (row: 1)\nCo...

xarray.Dataset

Dimensions:
- polarization_ids: 1
- spw_ids: 1
Coordinates: (2)
- polarization_ids
  (polarization_ids)
  int64
  0
```
array([0])
```
- spw_ids
  (spw_ids)
  int64
  0
```
array([0])
```
Data variables: (0)
Attributes: (4)
xds0 :
<xarray.Dataset> Dimensions: (row: 3510, uvw_index: 3, chan: 3, pol: 2, spw_id: 1, pol_id: 1) Coordinates: * chan (chan) float64 3e+09 3.4e+09 3.8e+09 * pol (pol) int32 5 8 * spw_id (spw_id) int64 0 * pol_id (pol_id) int64 0 Dimensions without coordinates: row, uvw_index Data variables: (12/24) UVW (row, uvw_index) float64 dask.array<chunksize=(3510, 3), meta=np.ndarray> FLAG (row, chan, pol) bool dask.array<chunksize=(3510, 3, 2), meta=np.ndarray> WEIGHT (row, pol) float32 dask.array<chunksize=(3510, 2), meta=np.ndarray> SIGMA (row, pol) float32 dask.array<chunksize=(3510, 2), meta=np.ndarray> ANTENNA1 (row) int32 dask.array<chunksize=(3510,), meta=np.ndarray> ANTENNA2 (row) int32 dask.array<chunksize=(3510,), meta=np.ndarray> ... ... TIME (row) datetime64[ns] dask.array<chunksize=(3510,), meta=np.ndarray> TIME_CENTROID (row) float64 dask.array<chunksize=(3510,), meta=np.ndarray> DATA (row, chan, pol) complex64 dask.array<chunksize=(3510, 3, 2), meta=np.ndarray> MODEL_DATA (row, chan, pol) complex64 dask.array<chunksize=(3510, 3, 2), meta=np.ndarray> CORRECTED_DATA (row, chan, pol) complex64 dask.array<chunksize=(3510, 3, 2), meta=np.ndarray> ROWID (row) int32 dask.array<chunksize=(3510,), meta=np.ndarray> Attributes: MS_VERSION: 2.0 column_descriptions: {'UVW': {'valueType': 'double', 'dataManagerType': ... info: {'type': 'Measurement Set', 'subType': 'simulator',... bad_cols: ['FLAG_CATEGORY']
SPECTRAL_WINDOW :
<xarray.Dataset> Dimensions: (row: 1, d1: 3) Dimensions without coordinates: row, d1 Data variables: (12/14) MEAS_FREQ_REF (row) int64 1 CHAN_FREQ (row, d1) float64 3e+09 3.4e+09 3.8e+09 REF_FREQUENCY (row) float64 3e+09 CHAN_WIDTH (row, d1) float64 4e+08 4e+08 4e+08 EFFECTIVE_BW (row, d1) float64 4e+08 4e+08 4e+08 RESOLUTION (row, d1) float64 4e+08 4e+08 4e+08 ... ... FREQ_GROUP_NAME (row) <U7 'Group 1' IF_CONV_CHAIN (row) int64 0 NAME (row) <U5 'sband' NET_SIDEBAND (row) int64 1 NUM_CHAN (row) int64 3 TOTAL_BANDWIDTH (row) float64 1.2e+09 Attributes: column_descriptions: {'MEAS_FREQ_REF': {'valueType': 'int', 'dataManager... info: {'type': '', 'subType': '', 'readme': ''} bad_cols: []
POLARIZATION :
<xarray.Dataset> Dimensions: (row: 1, d1: 2, d2: 2) Dimensions without coordinates: row, d1, d2 Data variables: CORR_TYPE (row, d1) int32 5 8 CORR_PRODUCT (row, d1, d2) int32 0 0 1 1 FLAG_ROW (row) bool False NUM_CORR (row) int64 2 Attributes: column_descriptions: {'CORR_TYPE': {'valueType': 'int', 'dataManagerType... info: {'type': '', 'subType': '', 'readme': ''} bad_cols: []
DATA_DESCRIPTION :
<xarray.Dataset> Dimensions: (row: 1) Coordinates: POLARIZATION_ID (row) int64 0 SPECTRAL_WINDOW_ID (row) int64 0 Dimensions without coordinates: row Data variables: FLAG_ROW (row) bool False Attributes: column_descriptions: {'FLAG_ROW': {'valueType': 'boolean', 'dataManagerT... info: {'type': '', 'subType': '', 'readme': ''} bad_cols: []

Image Simulated Dataset¶

[18]:

from casatasks import tclean
os.system('rm -rf simple.*')
tclean(vis=save_parms['ms_name'],imagename='simple',imsize=[400,400],cell=[5.0,5.0],specmode='cube',niter=1000,pblimit=0.1,pbmask=0.1,gridder='standard',stokes='RR')

[18]:

{}

[21]:

chan=0
from sirius.display_tools import display_image
print(np.abs(ms_xds.xds0.DATA[0,chan,0].values))
display_image(imname='simple.image',pbname='simple.pb',resname='simple.residual',ylim=[0.0,0.7],chan=chan)

0.6360993
Peak Intensity (chan0) : 0.6360901
PB at location of Intensity peak (chan0) : 0.6360983
max pixel location  (array([150]), array([250]), array([0]), array([0]))
Residual RMS : 0.0001055

[21]:

(0.636090099811554, 0.00010554427406769702)

[ ]: