Simulating ngVLA Data-CASA5.4.1

Introduction

The following tutorial shows how to create simulated data for the next generation Very Large Array (ngVLA). The ngVLA is composed of different subarrays that make up the current reference design. The configuration files for the different subarrays are found in the "ngVLA Configuration Tools" and can be used for simulations and calculations that investigate the scientific capabilities of the ngVLA. Each configuration (.cfg) file contains the name of the observatory, the antenna positions, the coordinate system of the antenna positions, and the diameter and the pad name of each antenna.

CASA provides two ways to create a simulation : the simobserve task and the sm toolkit. With these methods we can generate measurement sets, add thermal noise and predict model visibilities, from which we can explore the ngVLA’s imaging capabilities. The simobserve task is very user friendly but it is important to be aware that at the moment it has several limitations when making simulations for observatories other than ALMA and EVLA. For this reason, our suggested method for making simulations, specifically for advanced capabilities, is the sm toolkit. The sm toolkit provides much more flexibility in the observational parameters and is more compatible with observatories that are not recognized by CASA.

In this tutorial we present three example simulations: (i) simobserve using a model image, (ii) simobserve using a component list, and (iii) a sm toolkit simulation. For more information about component lists, see the CASA guides "here" and "here". The model image used in portions of this guide, a 93 GHz model of a protoplanetary disk, can be downloaded "here". In this guide we will create example continuum simulations at 93 GHz, consisting of a single channel observed for a total of 4 hours with a 60 second integration time (in order to keep the measurement set files small). We will then add to these simulations an amount of thermal noise that is representative of the ngVLA's continuum sensitivity. The array configuration used in this guide is the Main ngVLA subarray, which is composed of 214 18 m antennas and extends over a maximum baseline of 1005.4 km. The configuration file that we will use throughout this tutorial is called ngvla-main-revC.cfg and it can be found "here".

Estimating the scaling parameter for adding thermal noise

Simobserve's parameter to corrupt the simulated data is called thermalnoise and for interferometric data the allowed values are 'tsys-atm', 'tsys-manual' or False (see this "CASA guide" for more details). The option 'tsys-atm' is only applicable to ALMA since it uses site parameters which are specific to that observatory. For the option 'tsys-manual', it is necessary to supply several additional parameters which are required to construct an atmospheric model. In order to achieve the desired sensitivity without atmospheric modeling, we have chosen to corrupt the simulated data using the "sm.setnoise" function of the sm toolkit. In addition to the modes 'tsys-atm' and 'tsys-manual', this function also allows the option of 'simplenoise'. As the name indicates, 'simplenoise' adds random Gaussian noise to the visibilities based on a simple scaling factor. We will use this same technique for all three example simulations, i.e., those made with the simobserve task and those made with the sm toolkit.

In order to estimate the scaling factor for the 'simplenoise' parameter in sm.setnoise we use the following procedure:

The RMS noise ([math]\displaystyle{ \sigma_{NA} }[/math]) in an untapered, naturally-weighted Stokes I image will be approximately (see "setnoise function")

     [math]\displaystyle{ \sigma_{NA} \sim \frac{\sigma_{simple}}{ \sqrt{n_{ch}\,n_{pol}\,n_{baselines}\,n_{integrations} }} }[/math]   (1)

where [math]\displaystyle{ \sigma_{simple} }[/math] is the simplenoise parameter in sm.setnoise and corresponds to the noise per visibility, [math]\displaystyle{ n_{ch} }[/math] is the total number of channels across all spectral windows, [math]\displaystyle{ n_{pol} }[/math] is the number of polarizations in the measurement set and [math]\displaystyle{ n_{integrations} }[/math] is the number of correlator integration times in the measurement set (i.e., total on-source time / integration time). For the example simulations in this gude, the track time is 4 h and the integration time is 60 s, thus [math]\displaystyle{ n_{integrations}=240 }[/math]. Additionally, for these examples the total number of channels is 1 and the number of polarizations is 2. The number of baselines [math]\displaystyle{ n_{baselines} }[/math] is [math]\displaystyle{ N(N-1)/2 }[/math] where N is the number of antennas in the array. For the array configuration used in this guide (ngvla-main-revC.cfg), N=214 and therefore [math]\displaystyle{ n_{baselines}= 22791 }[/math].

If you already know the expected image noise ([math]\displaystyle{ \sigma_{NA} }[/math]) for your untapered, naturally-weighted image, you can solve for the scaling parameter [math]\displaystyle{ \sigma }[/math] in the above equation (1) and pass [math]\displaystyle{ \sigma_{simple} }[/math] to the simplenoise parameter in sm.setnoise.

If instead you want to calculate the expected sensitivity for an ngVLA image we suggest the below procedure:.

(i) Calculate the expected untapered, naturally weighted point source sensitivity ([math]\displaystyle{ \sigma_{NA} }[/math]) using one of the ngVLA performance tables. In Appendix D of "ngVLA memo #55 " there are key performance metrics for 6 subarrays which are tabulated as a function of frequency and resolution. For our example, we find in Table 10 of ngVLA memo #55 that the untapered, naturally weighted point source sensitivity of the Main interferometric array at 93 GHz is 0.83 uJy/beam for a 1 hour observation.

(ii) Scale that number to the desired observation length, in this case [math]\displaystyle{ t_{track}=4\,h }[/math]. Therefore, [math]\displaystyle{ \sigma_{NA} = 0.83/\sqrt{(t_{track}/1\,hour)} = 0.415\,\text{uJy/beam} }[/math].

(iii) Use the expected image noise ([math]\displaystyle{ \sigma_{NA} }[/math]) in the above equation (1) to solve for the scaling factor [math]\displaystyle{ \sigma_{simple} }[/math]. In this case, [math]\displaystyle{ \sigma_{simple}=0.415*\sqrt{1*2*22791*240} = 1.4\,\text{mJy} }[/math].

Once you have derived the scaling factor [math]\displaystyle{ \sigma_{simple} }[/math], run sm.setnoise and corrupt the visibilities of the noise-free measurement set. Your resulting untapered, naturally-weighted image will then have an RMS approximately equal to your desired image noise. Since it is not easy to undo this step, it is a good idea to make a copy of the noise-free measurement set before adding noise. The following mock example outlines this procedure:

# In CASA
os.system('cp -r noise_free.ms noisy.ms')
sm.openfromms('noisy.ms')
sm.setnoise(mode = 'simplenoise', simplenoise = sigma_simple)
sm.corrupt()
sm.done()

Note that this example will not execute without a measurement set named 'noise_free.ms' and a defined variable sigma_simple. See below for working examples of this procedure used in conjunction with the creation of simulated measurement sets and the prediction of model visibilities.

Simobserve using a model image

We will use the simobserve task to create our first noise-free measurement set (MS), using the configuration file and model image described in the Introduction.

# In CASA
## Using the configuration file obtained from the ngVLA's website
conf_file = 'ngvla-main-revC.cfg'
proj_name = 'ngVLA_214_ant_60s_noise_free'
model_file = 'ppmodel_image_93GHz'

simobserve(project = proj_name, 
                    skymodel = model_file+'.fits', 
                    inbright = '', 
                    incenter = '93GHz', 
                    inwidth = '1GHz' , 
                    setpointings = True, 
                    integration = '60s',  
                    direction = '',  
                    mapsize = '', 
                    obsmode = 'int', 
                    antennalist = conf_file, 
                    hourangle = 'transit', 
                    totaltime = '14400s',  
                    thermalnoise = '', 
                    graphics = 'none', 
                    overwrite = True)

project: Simobserve will create a folder with the project name in your current working directory, and this folder will contain all the resulting files including the noise-free ms.

skymodel: The input model image in Jy/pixel units. It could be a single image or a spectral cube. The simulated measurement set will inherit the number of channels, central frequency, and peak flux.

incenter: Central frequency of our simulation, in this case is 93 GHz. You only need to use this parameter if you want to override the reference frequency of the input model image.

inwidth: The bandwidth of the simulation.

%setpointings: allows simobserve to derive the pointing positions by its own algorithm. Given that the primary beam at Q-band is about 1 arcminutes (see the VLA observational status summary (OSS)), and the size of the model is less than an arcsecond, a single pointing will be adequate.

integration: We chose an integration time of 60 s in order to make the simulations faster. Note: An integration time of 60 s is not within the time-smearing limits for the ngVLA Main interferometric array, but for a source in the center of the map it should not be affected by time smearing. However, for future simulations and when considering sources at the edge of the primary beam this value needs to be reconsider and it should be << 0.24 s for this specific subarray.

direction: the center of the map. For a single pointing this is equivalent to the pointing center. Should we specify it?

obsmode: 'int' is used for interferometric data.

antennalist: the ngVLA configuration antenna position file in this case it corresponds to the ngVLA Main interferometric array. The files are available in the "ngVLA's website". Alternatively, the configuration files can be found in CASA in this path: casa.values()[0]['data']+'/alma/simmos/'.

hourangle: is used to simulate observations at a specific hour angle. We use 'transit' for culmination.

totaltime: This is the total track time, in this case the time on source which corresponds to 4 h.

thermalnoise: We leave this parameter empty for this noise-free simulation.

graphics: This will show graphics on the screen and/or save them as png files in the project directory. However, at the moment this is not working properly for baselines larger than a few hundred km. For this reason, we use graphics: 'none' in our example.

Now, to add thermal noise we do the following:

# In CASA
## Adding noise with 'simplenoise'
## set the noise level using the simplenoise parameter estimated in section 1
sigma_simple = '1.4mJy'

os.system('cp -r ngVLA_214_ant_60s_noise_free/ngVLA_214_ant_60s_noise_free.ngvla-main-revC.ms ngVLA_214_ant_60s_noisy.ms')
sm.openfromms('ngVLA_214_ant_60s_noisy.ms')
sm.setnoise(mode = 'simplenoise', simplenoise = sigma_simple)
## adds the noise: calculate random Gaussian numbers and add to visibilities
sm.corrupt()
sm.done()

Simobserve using a component list

Instead of using a model image we can generate a component list. Below is a simple example of a Gaussian source.

# In CASA
## Position of the source that we are observing. In this case  the source is in the same location where  the telescope is pointing
direction = me.direction(rf = 'J2000', v0 = qa.unit('00h0m0.0'), v1 = qa.unit('+24.0.0.000'))

## component list cl to make a model centered at the direction given above, and with a source of flux=3.2 Jy
cl.addcomponent(dir = direction, flux = 3.2, freq = '93GHz')

## name of the component list model 
cl.rename(filename = 'my_component.cl')

## close the component list 
cl.done()

Now, we can run simobserve using the generated component list:

# In CASA
## Using the configuration file obtained from the ngVLA's website
conf_file = 'ngvla-main-revC.cfg'
proj_name = 'ngVLA_214_ant_60s_noise_free'

simobserve(project = proj_name, 
                    skymodel = '', 
                    complist = 'my_component.cl' , 
                    compwidth = '1GHz',
                    setpointings = True, 
                    integration = '60s',  
                    direction = '',  
                    mapsize = '', 
                    obsmode = 'int', 
                    antennalist = conf_file, 
                    hourangle = 'transit', 
                    totaltime = '14400',  
                    thermalnoise = '', 
                    graphics = 'none', 
                    overwrite = True)

And we can add the thermal noise in the same way than in section 2.

sm toolkit using either a model image or a component list

# In CASA
## Using the configuration file obtained from the ngVLA's website
conf_file = 'ngvla-main-revC.cfg'

## If the configuration file is already included with the CASA distribution, use:
## configdir = casa.values()[0]['data']+'/alma/simmos/'

 
## Use simutil to read the .cfg file
from simutil import simutil
u = simutil()
xx,yy,zz,diam,padnames,telescope,posobs = u.readantenna(conf_file)
## or use this if the .cfg file is already part of the CASA distribution:
## xx,yy,zz,diam,padnames,telescope,posobs = u.readantenna(configdir+conf_file)



##############################################################
## Setting the observation framework  making the resources,
## similar to what we would do in the OPT when setting up 
## our observations.
##############################################################

## Simulate measurement set using the simulation utilities sm tool
ms_name = 'ngVLA_214_ant_60s_noise_free.ms'     ## Name of your measurement set
sm.open( ms_name )

## Get the position of the ngVLA using the measures utilities (me)
pos_ngVLA = me.observatory('ngvla')

## set the antenna configuration using the sm tool using the positions,
## diameter and names of the antennas as read from the configuration file
sm.setconfig(telescopename = telescope, x = xx, y = yy, z = zz,
                    dishdiameter = diam.tolist(), mount = 'alt-az',
                    antname = padnames,
                    coordsystem = 'global', referencelocation = pos_ngVLA)


## set the spectral windows, in this case is a single channel
## simulation with a channel resolution of 1 MHz and bandwidth of 1 MHz
sm.setspwindow(spwname = 'Band6', freq = '93GHz', deltafreq = '1MHz',
freqresolution = '1MHz', nchannels = 1, stokes = 'RR RL LR LL')

## set feed parameters for the antennas
sm.setfeed('perfect R L')

## set the field of observation that we are going to simulate
## (where the telescope is pointing), in this example we are using
## a Dec of +24deg
sm.setfield(sourcename = 'My source',
                  sourcedirection = ['J2000','00h0m0.0','+24.0.0.000'])

## set the limit of the observation for the antennas
sm.setlimits(shadowlimit = 0.001, elevationlimit = '8.0deg')

## weight to assign autocorrelation
sm.setauto(autocorrwt = 0.0)

## integration time or how often the array writes one visibility
integrationtime = '60s'
sm.settimes(integrationtime = integrationtime, usehourangle = True,
            referencetime = me.epoch('utc', 'today'))

## setting the observation time, which for our example is 4 h
starttime = '-2h'
stoptime = '2h'
sm.observe('My source', 'Band6', starttime = starttime, stoptime = stoptime)

Now we are going to write a scan of a source using the resource made above. Set the model either using a .fits model or a component list. After that we can predict the visibilities using sm.predict tool function. If using a component list follow the steps below:

## Using the same component list that we generated in section 2
## predicts the visibility of the source
sm.predict( complist = 'my_component.cl') 
sm.close()

However, if you want to use a model image instead please follow the steps below:

# In CASA
## To export the .fits file to .image 
model_file = 'ppmodel_image_93GHz'
importfits( fitsimage = model_file+'.fits', imagename = model_file+'.image')    

## Note: a warning is produced in CASA when running importfits about the image not having a beam or angular resolution. This is expected since the model is in units of Jy/pixel. Please ignore it.

## To predict the model visibilities
sm.predict( imagename = model_file+'.image') 

sm.close()

Note: the model image should have units of Jy/pixel and not Jy/beam which in that case will be the restored image.

Finally, in order to add thermal noise we do the following:

# In CASA
## Adding noise with 'simplenoise'
## set the noise level using the simplenoise parameter estimated in section 1
sigma_simple = '1.4mJy'
os.system('cp -r ngVLA_214_ant_60s_noise_free.ms ngVLA_214_ant_60s_noisy.ms')
sm.openfromms('ngVLA_214_ant_60s_noisy.ms')
sm.setnoise(mode = 'simplenoise', simplenoise = sigma_simple)
## adds the noise: calculate random Gaussian numbers and add to visibilities
sm.corrupt()
sm.done()

alt text

alt text

alt text

blah

Last checked on CASA Version 5.4.1.