ACA Simulation (CASA 3.3): Difference between revisions

Revision as of 12:37, 25 October 2011

To create a script of the Python code on this page see Extracting scripts from these tutorials.

ALMA 12m + ACA + Total Power

Multiple configurations, or an observation in one 12m configuration plus one ACA observation, or other combinations, can be simulated by running sim_observe for each, and then combining the Measurement Sets to run sim_analyze. Total power observations can be simulated either in an independent run of sim_observe, or along with an interferometric simulation. Note that if you simulate total power and an interferometric observation simultaneously with sim_observe, they must have the same set of pointing centers and the same integration and total time, which is probably not realistic. (For example it is generally recommended to observe a larger area by 1/2 primary beam in total power mode to combine with a 12m ALMA mosaic).

Simdata version for CASA 3.3

Explanation of the script

Set sim_observe as current task

Reset all parameters to default, and then set the project name to m51c

# Set sim_observe to default parameters
default("sim_observe")
# Our project name will be m51c, and all simulation products will be placed in a subdirectory m51c/
project="m51c"

Specify sky model image

We'll use an Halpha image of M51 as an example model sky. Download from here and place in your working directory.

# Model sky = Halpha image of M51 
modelimage         =  "m51ha.fits"

Although the image has a world coordinate system, we want to override most of the parameters.

We'll place the source in the southern hemisphere with the indirection parameter,
set the pixel size to 0.1arcsec, effectively moving the galaxy further away (M51 itself would require a quite large mosaic, and in any case we need for the input model pixels to be significantly smaller than the synthesized beam that we'll be simulating, or else we won't be learning anything)
set the frequency to 330GHz, and since its a 2D image we'll set the single "channel" width to be 50MHz, and peak brightness of 0.004 Jy/pixel - parameters plausible for observing an emission line in a galaxy.

# Set model image parameters:
indirection="B1950 23h59m59.96s -34d59m59.50s"
incell="0.1arcsec"
inbright="0.004"
incenter="330.076GHz"
inwidth="50MHz"

We need to decide where to point the telescope. The image is 2/3 arcsec in size, so we only need one pointing. We could put that in a text file ourself, or let simdata create the ascii pointing file for us.

setpointings       =  True
direction          =  "J2000 18h00m00.031s -22d59m59.6s"
mapsize            =  "0.76arcsec"

The default pointingspacing is fine: we'll only fit one pointing in the small mapsize the default calculation maptype hexagonal is ok too since only one will fit anyway.

We do want to calculate visibilities in a measurement set: let's do a 20 min snapshot observation using the "out20" ALMA antenna configuration:

observe            =  True
totaltime          =  "1200s"

Use appropriate antenna configurations based on desired angular resolution (configuration 20 - alma.out20.cfg in this case - is the largest "compact" configuration)

# It might be helpful to confirm the alma.out20.cfg file exists in the path defined below
# If you have a problem, this might be the first thing to check, if you haven't already
repodir=os.getenv("CASAPATH").split(' ')[0]
antennalist        =  repodir+"/data/alma/simmos/alma.out20.cfg"

Deconvolve the visibilities back into an image

image              =  True
vis                =  "$project.ms"
imsize             =  [192, 192]

Specify number of iteration of cleaning task with proper threshold and weighting

niter              =  10000
threshold          =  "1e-7Jy"
weighting          =  "natural"

We'd like to calculate a difference and fidelity image, and see some diagnostics:

analyze            =  True

And see the array but not the UV coverage:

showarray          =  True
showuv             =  False

Plot both to the screen and the png files with lots of messages:

graphics           =  "both"
verbose            =  True
overwrite = True

Run simdata

# This commands CASA to execute simdata
simdata()

Output results:

Input:	Predict:
Image:	Analyze:

↵ Simulating Observations in CASA

@@ Line 15: / Line 15: @@
 ====Explanation of the script====
-=====Set sim_observe as current task and reset all parameters=====
+=====Set sim_observe as current task=====
+Reset all parameters to default, and then set the project name to m51c
 <source lang="python">
-# Setting everything in sim_observe to original defaults
+# Set sim_observe to default parameters
 default("sim_observe")
+# Our project name will be m51c, and all simulation products will be placed in a subdirectory m51c/
+project="m51c"
 </source>
 =====Specify sky model image=====
+We'll use an Halpha image of M51 as an example model sky.  Download from here and place in your working directory.
 <source lang="python">
-# We'll use an Halpha image of M51 as an illustrative model
+# Model sky = Halpha image of M51
-# Make sure you are running CASA in the same directory as this file
+modelimage         =  "m51ha.fits"
-modelimage         =  "input50pc_672GHz.fits"
 </source>
-=====Image coordinate system can be verified=====
+Although the image has a world coordinate system, we want to override most of the parameters.
+* We'll place the source in the southern hemisphere with the indirection parameter,
+* set the pixel size to 0.1arcsec, effectively moving the galaxy further away (M51 itself would require a quite large mosaic, and in any case we need for the input model pixels to be significantly smaller than the synthesized beam that we'll be simulating, or else we won't be learning anything)
+* set the frequency to 330GHz, and since its a 2D image we'll set the single "channel" width to be 50MHz, and peak brightness of 0.004 Jy/pixel - parameters plausible for observing an emission line in a galaxy.
 <source lang="python">
-# This reports image header parameters in the Log Messages window
+# Set model image parameters:
-imhead("input50pc_672GHz.fits")
+indirection="B1950 23h59m59.96s -34d59m59.50s"
+incell="0.1arcsec"
+inbright="0.004"
+incenter="330.076GHz"
+inwidth="50MHz"
 </source>
-=====Image center can be identified=====
-<source lang="python">
-# These are tools found in the CASA toolkit
-# They are very useful, but the interface is not as straightforward as the tasks
-# You can find the tool reference manual here: http://casa.nrao.edu/docs/CasaRef/CasaRef.html
-# The following command is used to open an image (this is part of the image analysis toolkit)
-# Whenever data are being viewed/manipulated by tools, what is being operated on needs to be explicitly opened
-# and closed (i.e. an image, a table, etc.)
-ia.open("input50pc_672GHz.fits")
-#  Out[9]: True
-#
-# Reports the length of each axis in the opened image
-ia.shape()
-#  Out[11]: [257L, 257L, 1L, 1L]
-#
-# This command converts from pixel (our source file) to world coordinates (something usable by simdata)
-ia.toworld([128.5,128.5])
-#  Out[12]:
-#{'numeric': array([  4.71239120e+00,  -4.01423802e-01,   1.00000000e+00,
-#         6.72000001e+11])}
-#
-# Formats the coordinate just converted into hms
-qa.formxxx("4.71239120rad",format='hms',prec=5)
-#  Out[13]: '18:00:00.03052'
-#
-# Formats one of the other coordinates into dms
-qa.formxxx("-0.401423802rad",format='dms',prec=5)
-#  Out[14]: '-022.59.59.602743'
-#
-# Final housekeeping by closing the image tool
-# The image tool will now be detached from the image
-ia.close()
-</source>
-=====Brightness scale can be viewed with 'imstat' task=====
-<source lang="python">
-# Default parameters are adequate for this
-imstat("input50pc_672GHz.fits")
-# ...
-#  'max': array([  6.52469971e-05]),
-# ...
-# that's 0.0652 mJy/pixel.
-</source>
-=====Let's call our project psim2=====
-<source lang="python">
-# This defines the root prefix for any output files from simdata
-project            =  "psim2"
-</source>
-=====We'll leave the sky model the way it is: simdata will create psim2.skymodel CASA image since this model is a fits file, and most but not all of CASA routines can operate directly on fits=====
-<source lang="python">
-skymodel           =  "input50pc_672GHz.fits"
-</source>
 =====We need to decide where to point the telescope.  The image is 2/3 arcsec in size, so we only need one pointing. We could put that in a text file ourself, or let simdata create the ascii pointing file for us.=====
 <source lang="python">