Imaging simulations

This document serves as a brief user guide on performing simulations for dynamic range and noise estimation. We will be using standard packages and a Python environment. This assumes that the user has installed CASA and Meqtrees from the standard KERN-5 PPA on an Ubuntu (or derivative) 18.04 system. Please follow the instructions on the KERN 5 launchpad for details on adding the repository through the APT package management system.

We start by using the simms wrapper around CASA to simulate a dataset into the Measurement Set Standard v2.0 (CASA Memo 229). The MeerKAT antenna positions are already built-in. You can obtain a copy from (here we use release v1.2.0). We will be doing a transit-spanning 0.5 hour simulation with a fictitious catalogue centred at 04h00m00.0s, -33d00m00s (J2000), with correlator settings matching the 4k channel L-band MeerKAT correlation mode with 10 minute scans and 8s integration intervals (see the documentation for a discussion on all possible switches).

> simms -dir "J2000,04h00m00.0s,-33d00m00s" -T meerkat -dt 8 -st 0.5 -sl 0.16667 -nc 4096 -f0 856MHz -df 208.9843kHz -pl XX YY -n

Empty MS '' created

Validating ...

MS validated


You can verify that your empty database has the correct shape. The timespan may be different for your dataset.

> casa --nologger --log2term --nogui -c "listobs(vis='', verbose=False)"


2019-12-23 14:08:13 INFO listobs::::+ ##########################################

2019-12-23 14:08:13 INFO listobs::::+ ##### Begin Task: listobs        #####

2019-12-23 14:08:13 INFO listobs:::: listobs(vis="",selectdata=True,spw="",field="",antenna="",

2019-12-23 14:08:13 INFO listobs::::+        uvrange="",timerange="",correlation="",scan="",intent="",

2019-12-23 14:08:13 INFO listobs::::+        feed="",array="",observation="",verbose=True,listfile="",

2019-12-23 14:08:13 INFO listobs::::+        listunfl=False,cachesize=50,overwrite=False)

2019-12-23 14:08:13 INFO listobs::ms::summary ================================================================================

2019-12-23 14:08:13 INFO listobs::ms::summary+          MeasurementSet Name:  /home/hugo/userman_sims/  MS Version 2

2019-12-23 14:08:13 INFO listobs::ms::summary+   ================================================================================

2019-12-23 14:08:13 INFO listobs::ms::summary+  Observer: CASA simulator Project: CASA simulation  

2019-12-23 14:08:13 INFO listobs::ms::summary+   Observation: meerkat

2019-12-23 14:08:13 INFO listobs::MSMetaData::_computeScanAndSubScanProperties   Computing scan and subscan properties...

2019-12-23 14:08:13 INFO listobs::ms::summary Data records: 453600   Total elapsed time = 1800.02 seconds

2019-12-23 14:08:13 INFO listobs::ms::summary+  Observed from   22-Dec-2019/20:15:38.3   to   22-Dec-2019/20:45:38.3 (UTC)

2019-12-23 14:08:13 INFO listobs::ms::summary

2019-12-23 14:08:13 INFO listobs::ms::summary+  ObservationID = 0     ArrayID = 0

2019-12-23 14:08:13 INFO listobs::ms::summary+ Date    Timerange (UTC)      Scan  FldId FieldName         nRows SpwIds   Average Interval(s) ScanIntent

2019-12-23 14:08:13 INFO listobs::ms::summary+ 22-Dec-2019/20:15:38.3 - 20:25:38.3 1  0 00                  151200  [0]  [8] [OBSERVE_TARGET.ON_SOURCE]

2019-12-23 14:08:13 INFO listobs::ms::summary+             20:25:38.3 - 20:35:38.3 2  0 00                  151200  [0]  [8] [OBSERVE_TARGET.ON_SOURCE]

2019-12-23 14:08:13 INFO listobs::ms::summary+             20:35:38.3 - 20:45:38.3 3  0 00                  151200  [0]  [8] [OBSERVE_TARGET.ON_SOURCE]

2019-12-23 14:08:13 INFO listobs::ms::summary           (nRows = Total number of rows per scan)

2019-12-23 14:08:13 INFO listobs::ms::summary Fields: 1

2019-12-23 14:08:13 INFO listobs::ms::summary+ ID   Code Name            RA           Decl       Epoch   SrcId  nRows

2019-12-23 14:08:13 INFO listobs::ms::summary+ 0     00              04:00:00.000000 - J2000   0     453600

2019-12-23 14:08:13 INFO listobs::ms::summary Spectral Windows:  (1 unique spectral windows and 1 unique polarization setups)

2019-12-23 14:08:13 INFO listobs::ms::summary+ SpwID  Name   #Chans   Frame   Ch0(MHz)  ChanWid(kHz)  TotBW(kHz) CtrFreq(MHz)  Corrs  

2019-12-23 14:08:13 INFO listobs::ms::summary+ 0  00  4096   TOPO 856.000   208.984 855999.7   1283.8954   XX  YY

2019-12-23 14:08:13 INFO listobs::ms::summary Sources: 1

2019-12-23 14:08:13 INFO listobs::ms::summary+ ID   Name            SpwId RestFreq(MHz)  SysVel(km/s)

2019-12-23 14:08:13 INFO listobs::ms::summary+ 0 00              0 856        0        

2019-12-23 14:08:13 INFO listobs::ms::summary Antennas: 64:



Our elevation range automatically adjusts for optimal half-hour uv coverage and spans transit roughly in the middle of the observation. This can be adjusted according to your needs with the start and end time switches in simms:


Before proceeding to simulation we need to include the antenna far-field response. Here we use the Eidos package which includes Zernike model coefficients for the reconstruction of smoothed versions of L-band holographic measurements (K. Asad & M. de Villiers, private communication). The package can be obtained from (we will use version v1.0.4). The beam is evaluated every 5 MHz over the entire L-band bandwidth out to 4° (a maximum of 10° is possible at the time of writing). The package will dump 8 files - pairs of real and imaginary voltage responses for each of the correlation products in a 2x2 Jones matrix. 

> eidos -d 4 -r 0.015625 -f 856 1712 5 -P mk.holo.lband.4deg -o8

The I Mueller power response is shown below in logscale (the -10dB power point is roughly 2°15' in diameter at the lowest frequency and 1°10' in diameter at the highest frequency of the band):

With an empty database and antenna response in hand we can now initialize a basic fictitious catalog with two strong sources off-centre at beyond the primary beam of the antenna far-field response. Call this sources.txt with the following format:

#format: ra_h ra_m ra_s dec_d dec_m dec_s i freq0 spi

04 00 10.0 -33 00 00.0 1.0e-2 1.4e9 -0.75

04 10 10.0 -33 12 00.0 6.5e-2 1.4e9 -0.75

03 57 10.0 -32 55 00.0 4.0e-3 1.4e9 -0.75

03 58 10.0 -33 07 00.0 3.0e-4 1.4e9 -0.75

04 03 10.0 -33 20 00.0 5.0e-2 1.4e9 -0.75

04 07 10.0 -32 50 00.0 8.0e-4 1.4e9 -0.75

04 15 04.0 -32 00 00.0 8.0e-1 1.4e9 -0.75

04 00 00.0 -34 00 00.0 3.0 1.4e9 -0.75

We will use Meqtrees in interactive mode. Pipelining modes are available though we will not cover their use here. Start the browser with the following command and start the meqserver in multithreaded mode with the number of threads given via the -mt parameter:




From the TDL menu select Load TDL Script and navigate to the within the Siamese submodule of Cattery. The KERN installation will be in the Python 2.7 dist-packages as shown below:

The TDL compile time options will be displayed once the recipe is loaded. Notice that there is a capacity here to construct a Radio Interferometer Measurement Equation (RIME) that is both direction independent and direction dependent. We will demonstrate only the direction dependent antenna response options here. Start by loading in the empty database in the MS section, the sky catalogue into the Sky model -> TiggerSkyModel module and the 8 cube antenna response output by Eidos into the E Jones Siamese.OMS.pybeams_fits module. Ensure to change the CTYPE of L and M axis to ‘px’ and ‘py’ as per the antenna voltage response FITS file headers. We will also apply apparent sky rotation due to the ALT-AZ mounts and periodic pointing errors (here assumed to be 200 arcseconds). Lastly, we will add a realistic average noise of around 425 Jy SEFD across the band. Ensure the various parent options are ticked. Hit “Compile” once the options are filled:

Once the compute graph is constructed runtime parameters will be shown. Change the output column to the (empty) DATA column and increase the tile size to 128 time steps (or less if this does not fit into your memory constraints). Hit “Run simulation” to start simulating into the database:

Once the simulation has started, constituent interferometer and pointing error inspectors are available from the Bookmarks tab.

Once completed  the field can be imaged using your favourite imager. We will create an MFS map using WSClean to get an estimate of the instrument induced errors. 

> wsclean --join-channels -size 8192 8192 -scale 1.3asec -channelsout 8 -name briggs.-0.3 -niter 5000 -nmiter 10 -mgain 0.85 -weight briggs -0.3 -datacolumn DATA

The apparent source variation and spectrum steepening at the roll off points in the antenna response over such a wide bandwidth induces large unmodelled errors in the total system response. These errors manifest as direction dependent effects around the bright off-axis sources and are especially visible where the field of view is at its widest at the lower part of the instrumental bandwidth. Here we show the MFS map and the area around the bright source to the south of the pointing centre at the lowest and highest frequency respectively: 

MFS map showing pointing error induced direction dependent effects at the steepening edge of the primary beam of the antenna far-field response

Left: Direction Dependent Effects at the widest angle of view at the bottom of L-band and Right: direction dependent effects after suppression by the antenna response at the top of L-band. The image scales are identical.

Experience has shown that, although possible using a combination of packages like DDFacet, killMS and Cubical, effectively calibrating for direction dependent gains can be very challenging due to modelling inaccuracies over a large fractional bandwidth and SNR limitations at the high end of the band. A far easier mitigation approach is to preempt direction dependent calibration effects and use the antenna response to effectively suppress bright sources, or if necessary place bright sources more centrally within the primary beam of the antenna response to limit modulation from time-variable effects.