Simulation Guide Component Lists (CASA 3.3)
↵ Simulating Observations in CASA
This guide is applicable to CASA version 3.3.
To create a script of the Python code on this page see Extracting scripts from these tutorials.
Explanation of the guide
When planning an interferometric observation it is often useful to simulate observations of very simple objects, like point sources, Gaussians, and disks. In CASA, observations can be simulated using task sim_observe and analyzed using task sim_analyze. This guide will demonstrate how to simulate an ALMA observations of a Gaussian and some point sources using these tasks as well as the CASA Toolkit.
We begin by using component lists in the Toolkit to create an image of a Gaussian flux distribution, and this will be saved as a FITS file. The fits file will them be "observed" using sim_observe and sim_analyze along with four point sources, added via the componentlist parameter. Finally, we show how the same observations could have been done without any skymodel in sim_observe, instead using only component lists.
Getting Started
To get started you need CASA version 3.3
To install CASA, follow the instructions given on the Obtaining CASA page.
CASA Basics
CASA is the post-processing package for ALMA and EVLA and can handle both interferometric and single dish data. To get a brief introduction to sim_observe and sim_analyze, the tasks within CASA that we will use here, go to Simulation Guide for New Users. To learn more about CASA in general, go to the CASA homepage. Walk-throughs of CASA data reduction for a variety of data sets can be found on the CASA Guides website.
Once you have installed CASA, you can launch it by typing "casapy" at the prompt or by double-clicking on the icon, depending on your system and preferences.
Making a Simple FITS Image
Here we show how to create a simple FITS image using the CASA tasks and the toolkit. The example here will be that of a Gaussian flux distribution. Enter the following lines at the CASA prompt:
# In CASA
direction = "J2000 10h00m00.0s -30d00m00.0s"
cl.done()
cl.addcomponent(dir=direction, flux=1.0, fluxunit='Jy', freq='230.0GHz', shape="Gaussian",
majoraxis="0.1arcmin", minoraxis='0.05arcmin', positionangle='45.0deg')
ia.fromshape("Gaussian.im",[256,256,1,1],overwrite=True)
cs=ia.coordsys()
cs.setunits(['rad','rad','','Hz'])
cell_rad=qa.convert(qa.quantity("0.1arcsec"),"rad")['value']
cs.setincrement([-cell_rad,cell_rad],'direction')
cs.setreferencevalue([qa.convert("10h",'rad')['value'],qa.convert("-30deg",'rad')['value']],type="direction")
cs.setreferencevalue("230GHz",'spectral')
cs.setincrement('1GHz','spectral')
ia.setcoordsys(cs.torecord())
ia.setbrightnessunit("Jy/pixel")
ia.modify(cl.torecord(),subtract=False)
exportfits(imagename='Gaussian.im',fitsimage='Gaussian.fits',overwrite=True)
The first line defines a string "direction" which will be the center of the Gaussian flux distribution.
cl.done closes any open component lists, if any.
cl.addcomponent creates a new component centered at "direction", with a flux of 1 Jy at a frequency of 230 GHz, a Gaussian shape of 0.1' by 0.05' with a position angle of 45 degrees.
ia.fromshape creates a new, empty CASA image with the name and dimensions given.
cs.coordsys gets the coordinate system of the image.
cs.setunits defines the units of the four axes of the new CASA image.
cell_rad will be the cell size and units in this CASA image, 0.1"
cs.setincrement tells CASA that RA increases to the right, Dec increases going up, and in a few lines that the one channel is 1 GHz wide.
cs.setreferencevalue sets the center of the image in RA, Dec, and frequency.
ia.setcoordsys puts the coordinates and frequencies into the image header.
ia.setbrightnessunit defines the brightness unit (Jy per pixel) of the CASA image.
ia.modify puts the Gaussian component into the image
exportfits writes the resultant CASA image as a FITS file (not strictly necessary for this guide, but useful to know in general).
As usual, more information can be found via the help in CASA
# In CASA
help(ia.modify) # syntax for help with toolkit or CASA tasks
help("exportfits") # syntax for help with CASA tasks, but not the toolkit
Simulating Observations with a FITS Image and a Component List
One use for component lists would to simulate the effect of having one or more point sources added to an input image, with the goal of finding out the effect on the simulated observations. For instance, one might want to know if a faint point source would be detectable if there is extended emission around it. Conversely, one might want to know if the artifacts from imaging a field with a bright point source would make a project tricky to carry out. In this example, we will component lists in the sim_observe task add four point sources to an input image. The input image will be the Gaussian flux distribution created above.
First we create the four point sources using the CASA toolkit.
# In CASA
cl.done()
cl.addcomponent(dir="J2000 10h00m00.08s -30d00m02.0s", flux=0.1, fluxunit='Jy', freq='230.0GHz', shape="point")
cl.addcomponent(dir="J2000 09h59m59.92s -29d59m58.0s", flux=0.1, fluxunit='Jy', freq='230.0GHz', shape="point")
cl.addcomponent(dir="J2000 10h00m00.40s -29d59m55.0s", flux=0.1, fluxunit='Jy', freq='230.0GHz', shape="point")
cl.addcomponent(dir="J2000 09h59m59.60s -30d00m05.0s", flux=0.1, fluxunit='Jy', freq='230.0GHz', shape="point")
cl.rename('point.cl')
cl.close()
We use cl.done to begin and end this sequence to close any open component list. The cl.addcomponent commands create point sources that are 0.1 Jy at 230 GHz at the coordinates given in each line.
We use sim_observe and sim_analyze to make the simulated observations of these point sources and the Gaussian flux distribution given in the FITS file we made previously in this guide.
default("sim_observe")
project = "cycle0_compact_8h"
skymodel = "Gaussian.fits"
inwidth = "1GHz"
complist = 'point.cl'
compwidth = '1GHz'
direction = "J2000 10h00m00.0s -30d00m00.0s"
repodir = os.getenv("CASAPATH").split(' ')[0]
antennalist = repodir+'/data/alma/simmos/alma.cycle0.compact.cfg'
totaltime = "28800s"
mapsize = ["20arcsec","20arcsec"]
sim_observe()
default("sim_analyze")
project = "cycle0_compact_8h"
vis="cycle0_compact_8h.alma.cycle0.compact.ms"
imsize = [100,100]
imdirection = "J2000 10h00m00.0s -30d00m00.0s"
cell = '0.2arcsec'
niter = 5000
treshold = '40.0mJy/beam'
analyze = True
sim_analyze()
