M51 at z = 0.1 and z = 0.3 (CASA 3.3): Difference between revisions

From CASA Guides
Jump to navigationJump to search
No edit summary
 
(55 intermediate revisions by 2 users not shown)
Line 13: Line 13:
This tutorial presents a simulation of ALMA observations of the well-known galaxy M51 in the CO 1-0 transition. This galaxy is located relatively nearby (luminosity distance = 9 Mpc), but, in this simulation, we will model how it would appear at redshift z = 0.1 (luminosity distance = 460 Mpc).  
This tutorial presents a simulation of ALMA observations of the well-known galaxy M51 in the CO 1-0 transition. This galaxy is located relatively nearby (luminosity distance = 9 Mpc), but, in this simulation, we will model how it would appear at redshift z = 0.1 (luminosity distance = 460 Mpc).  


The goal of this tutorial is to provide a ''complete'' run-through of a relatively simple simulation. Included in this simulation are the effects of (''u'', ''v'') sampling of a 50-antenna ALMA, the primary beam of the ALMA, and thermal noise. Neither calibration overheads nor errors are included, and so this simulation should be viewed as optimistic.
The goal of this tutorial is to provide a ''complete'' run-through of a relatively simple simulation. The simulation will include effects of (''u'', ''v'') sampling of a 50-antenna ALMA, the primary beam of the ALMA, and thermal noise. Neither calibration overheads nor errors are included, and so this simulation should be viewed as optimistic.  The CASA task <tt>sim_observe</tt> will generate simulated visibility measurements from an input model.  The task <tt>sim_analyze</tt> will image the visibilities and generate diagnostic plots.


For this tutorial, we'll use the [http://nedwww.ipac.caltech.edu/level5/March02/SONG/SONG.html BIMA SONG] ([http://adsabs.harvard.edu//abs/2003ApJS..145..259H Helfer et al. 2003]) observations of M51 as the basis for the model. Grab the file [http://nedwww.ipac.caltech.edu/level5/March02/SONG/NGC5194.bima12m.cm.fits.gz NGC5194.bima12m.cm.fits.gz] and uncompress it in a working directory.
For this tutorial, we'll use the [http://nedwww.ipac.caltech.edu/level5/March02/SONG/SONG.html BIMA SONG] ([http://adsabs.harvard.edu//abs/2003ApJS..145..259H Helfer et al. 2003]) observations of M51 as the basis for the model. Grab the file [http://nedwww.ipac.caltech.edu/level5/March02/SONG/NGC5194.bima12m.cm.fits.gz NGC5194.bima12m.cm.fits.gz] and uncompress it in a working directory.
Line 22: Line 22:
</source>
</source>


=== Load data and Initialize <tt>sim_observe</tt> ===  
=== Initializing the Data and Tasks ===  


For convenience, we store the name of the resulting measurement set into the Python global <tt>cubeName</tt>. The task <tt>importfits</tt> is needed to convert the image FITS file into a CASA image. Here we invoke <tt>importfits</tt> as a function, where the task inputs are set as function input parameters. The task will produce a new directory with the name designated by imagename (in this case <tt>m51-song.image</tt>) in the current working directory.
Task <tt>sim_observe</tt> takes an input model in the form of a CASA image. So, we begin by converting the BIMA SONG FITS cube into a CASA image using task <tt>importfits</tt>. Here we invoke <tt>importfits</tt> as a function, where the task inputs are set as function input parameters. For convenience, we store the name of the resulting image into the Python variable <tt>cubeName</tt>. The task will produce the new image (which is actually a directory) in the current working directory.


<source lang="python">
<source lang="python">
# in CASA
# in CASA
cubeName = 'm51-song.image'
cubeName = 'm51-song.image'
# If m51-song.image exists, it needs to be removed for importfits to work as set
os.system('rm -rf '+cubeName)
importfits(fitsimage='NGC5194.bima12m.cm.fits', imagename=cubeName)
importfits(fitsimage='NGC5194.bima12m.cm.fits', imagename=cubeName)
</source>
</source>


The task <tt>sim_observe</tt> will generate the simulated measurement set using the input modelWe initialize the task like so.
We initialize the simulation tasks by setting their parameters to the default values.  We then make the new CASA image the input model to <tt>sim_observe</tt> by setting parameter '''skymodel'''.
 
   
<source lang="python">
<source lang="python">
# in CASA
# in CASA
default("sim_observe")
default("sim_observe")
skymodel=cubeName
default("sim_analyze")
skymodel = cubeName
</source>
</source>


Line 48: Line 51:
=== Cosmology Calculations ===
=== Cosmology Calculations ===


Next we'll set up some python globals to handle the scaling of the model coordinates and flux densities appropriate for new redshift. We'll primarily need the angular size and luminosity distances for a given cosmology. To keep things simple, we'll use [http://www.astro.ucla.edu/~wright/CosmoCalc.html Ned Wright's CosmoCalc] ([http://adsabs.harvard.edu/abs/2006PASP..118.1711W 2006]) with the default cosmology; redshifts were collected from [http://nedwww.ipac.caltech.edu/cgi-bin/nph-objsearch?objname=NGC5194&extend=no&hconst=73&omegam=0.27&omegav=0.73&corr_z=1&out_csys=Equatorial&out_equinox=J2000.0&obj_sort=RA+or+Longitude&of=pre_text&zv_breaker=30000.0&list_limit=5&img_stamp=YES NED].  
Next we'll set up some python globals to handle the scaling of the model coordinates and flux densities appropriate for new redshift. We'll primarily need the angular size and luminosity distances for a given cosmology. To keep things simple, we'll use [http://www.astro.ucla.edu/~wright/CosmoCalc.html Ned Wright's CosmoCalc] ([http://adsabs.harvard.edu/abs/2006PASP..118.1711W 2006]) with the default cosmology; redshifts were collected from [http://nedwww.ipac.caltech.edu/cgi-bin/nph-objsearch?objname=NGC5194&extend=no&hconst=73&omegam=0.27&omegav=0.73&corr_z=1&out_csys=Equatorial&out_equinox=J2000.0&obj_sort=RA+or+Longitude&of=pre_text&zv_breaker=30000.0&list_limit=5&img_stamp=YES NED].


<source lang="python">
<source lang="python">
#z's
# in CASA
# Distinguish between z_lsrk, which sets the observed frequency for scaling, from z_cmb,  
# Distinguish between z_lsrk, which sets the observed frequency for scaling,  
#      which is needed to get cosmological distances.
# and z_cmb, which is needed to get cosmological distances.
z_old_cmb = 0.002122 # CMB-referenced z for cosmological distances from NED
z_old_cmb = 0.002122 # CMB-referenced z for cosmological distances from NED
z_old_lsrk = 0.001544 # from NED
z_old_lsrk = 0.001544 # from NED
Line 67: Line 70:
</source>
</source>


The convention is ''old'' refers to M51 as observed at its proper redshift, and ''new'' refers to the new, higher redshift for our model.
We will use the convention where ''old'' refers to M51 as observed at its proper redshift, and ''new'' refers to the new, higher redshift for our model.


=== Preparing the Model ===
=== Preparing the Model ===


The next step is to scale the M51 data cube into a model cube appropriate for <tt>sim_observe</tt>. First, we establish the file naming convention by setting the task parameter <tt>project</tt>.
The next step is to scale the M51 data cube into a model cube appropriate for <tt>sim_observe</tt>. First, we establish the file naming convention by setting the task parameter '''project'''.  Thinking ahead, we might want to rerun this simulation for other redshifts.  So, to make this easy we will include the redshift in the '''project''' name: we'll write 0.1 as <tt>p1</tt> (or "point one") and assign it to Python variable <tt>suffix</tt>.  Then we'll use <tt>suffix</tt> to define '''project'''.


<source lang="python">
<source lang="python">
# in CASA
# in CASA
# z = 0.1, or point-1; useful to distinguish from repeated simulations at different z's
suffix = "-p1"
suffix = "-p1"
# project ID to assign output filenames
project = "M51-AtZ" + suffix   
project = "M51-ATZ" + suffix   
</source>
</source>


==== Flux Density Scaling ====
==== Flux Density Scaling ====


Our goal in this step is to set the <tt>inbright</tt> parameter of <tt>sim_observe</tt>, which will scale the data cube appropriately for its new luminosity distance.  The task wants models in units of Jy / pixel, but the BIMA SONG cube has units of Jy / beam. That's an easy conversion.
Our goal in this step is to set the '''inbright''' parameter of <tt>sim_observe</tt>, which will scale the data cube appropriately for its new luminosity distance.  The task wants models in units of Jy / pixel, but the BIMA SONG cube has units of Jy / beam. We'll use the task <tt>imhead</tt> to get the BIMA SONG synthesized beam dimensions and we'll use the <tt>qa</tt> tool to convert the dimensions from radians to arcseconds.


<source lang="python">
<source lang="python">
Line 98: Line 98:
</source>
</source>


Next, scale the flux for the new luminosity distance. Make the approximation that each pixel is a point source, and use the inverse square law to scale. <tt>sim_observe</tt> can scale the peak surface brightness using the parameter <tt>inbright</tt>.
Next, scale the flux for the new luminosity distance. Make the approximation that each pixel is a point source, and use the inverse square law to scale. We then set '''inbright''' to the newly scaled peak flux.


<source lang="python">
<source lang="python">
Line 114: Line 114:
==== Angular Size Scaling ====
==== Angular Size Scaling ====


The sky coordinates axes of the model need to be adjusted (1) to place M51 in the southern hemisphere and (2) to correct for the new angular size distance. These tasks can be accomplished by the <tt>sim_observe</tt> parameters <tt>indirection</tt> (to change the location of M51 on the sky) and <tt>incell</tt> (to change the angular scale of an input pixel).  
The sky coordinates axes of the model need to be adjusted (1) to place M51 in the southern hemisphere and (2) to correct for the new angular size distance. These tasks can be accomplished by the <tt>sim_observe</tt> parameters '''indirection''' (to change the location of M51 on the sky) and '''incell''' (to change the angular scale of an input pixel).  
To perform task (1), we'll just flip the sign of the declination using [[imhead]]. We'll use the [http://casa.nrao.edu/docs/casaref/quanta-Module.html#x409-4100001.4 qa] tool to convert radians to sexagesimal units.
To perform task (1), we'll just flip the sign of the declination using [[imhead]]. We'll use the [http://casa.nrao.edu/docs/casaref/quanta-Module.html#x409-4100001.4 qa] tool to convert radians to sexagesimal units.


Line 129: Line 129:
</source>
</source>


Next, adjust the pixel scale for the new angular size distance. To perform this adjustment, we'll use [[imhead]] with <tt>mode = "get"</tt> to read in the original pixel scale, then use the <tt>sim_observe</tt> input parameter <tt>incell</tt> to designate the new pixel size.
Next, adjust the pixel scale for the new angular size distance. To perform this adjustment, we'll use [[imhead]] with '''mode = "get"''' to read in the original pixel scale, then use the <tt>sim_observe</tt> input parameter '''incell''' to designate the new pixel size.


<source lang="python">
<source lang="python">
Line 162: Line 162:
</source>
</source>


=== <tt>sim_observe</tt> ===
=== Map Size and Mosaicking ===


The original BIMA SONG image is about 480 arcseconds across; scale this image size to the new redshift.  We will use this value later.
The original BIMA SONG image is about 480 arcseconds across; scale this image size to the new redshift and assign to parameter '''mapsize'''.  We will also use this value later.


<source lang="python">
<source lang="python">
# in CASA
# in CASA
imSize = 480.0 * (da_old / da_new) # arcseconds
imSize = 480.0 * (da_old / da_new) # arcseconds
mapsize = "%farcsec" % imSize
</source>
</source>


Line 177: Line 178:
primaryBeam = 17.0 * (300e9 / newFreq) # arcseconds
primaryBeam = 17.0 * (300e9 / newFreq) # arcseconds
pointingspacing = "%farcsec" % (primaryBeam / 2.0)  
pointingspacing = "%farcsec" % (primaryBeam / 2.0)  
mapsize = "%farcsec" % imSize
</source>
</source>


We also need to estimate the desired synthesized beam size. We don't want the new beam to be so large so that we cannot resolve a rotation curve, but we also don't want it to be so small that we effectively resolve out the BIMA SONG data. The BIMA SONG beam was 5 arcsec, and so as a reasonable guess we'll adopt the equivalent of a 15" beam at its true distance(3 times coarser than the BIMA SONG beam) and then scale appropriately to z = 0.1.
=== Array Configuration ===


<source lang="python">
We want to set '''antennalist''' such that <tt>sim_observe</tt> will decide on an appropriate ALMA configuration based on the desired beam size.  Thus, we first must estimate the desired synthesized beam size.  We don't want the new beam to be so large so that we cannot resolve a rotation curve, but we also don't want it to be so small that we effectively resolve out the BIMA SONG data. The BIMA SONG beam was 5 arcsec, and so as a reasonable guess we'll adopt the equivalent of a 15" beam at its true distance (3 times coarser than the BIMA SONG beam) and then scale appropriately to z = 0.1. 
# in CASA
beamNew = 15.0 * (da_old / da_new)
</source>


We want pixels that sample the beam at least 3 times for stable deconvolution; we'll use 4 times sampling, rounded off to the nearest milliarcsec.
Looking ahead to imaging our simulated data, it's important to sample the beam at least 3 times for stable deconvolution.  We'll sample the beam by 4 times, rounded off to the nearest milliacsecond.


<source lang="python">
<source lang="python">
# in CASA
# in CASA
beamNew = 15.0 * (da_old / da_new)
pixelSize = round(beamNew * 1000.0 / 4.0) / 1000.0
pixelSize = round(beamNew * 1000.0 / 4.0) / 1000.0
</source>
</source>
Line 201: Line 199:
</source>
</source>


Now we know both the image size and pixel size in arcseconds, but <tt>sim_analyze</tt> wants the ratio: the number of pixels along the RA or Dec axis. To keep the image from becoming too small, set the minimum image size to be 256 pixels.
Now we can set parameter '''antennalist''' using the beam size (but see [[#Aside: Antenna Configurations|Aside: Antenna Configurations]] below).


<source lang="python">
<source lang="python">
# in CASA
# in CASA
imSizePix = int(round(imSize / pixelSize))
antennalist = "alma;%farcsec" % beamNew
if imSizePix < 256: imSizePix = 256
</source>
</source>


Let <tt>sim_observe</tt> decide on an appropriate ALMA configuration based on the desired beam size. Set the parameter '''antennalist''' as follows (but see [[#Other Antenna Configurations|Other Antenna Configurations]] below).
=== Integration and Total Time ===


<source lang="python">
We will set the integration time to 500 seconds and the total integration time to 8 hours.  For a real observation, 10 seconds will be a more appropriate integration time.  However, we can decrease the simulation CPU time by increasing the integration interval (See [[Etime study]] for details).  Whenever simulating a new data set, we recommend beginning with a large integration time for initial testing and then lowering the integration time to generate a more realistic measurement set.   
# in CASA
antennalist = "alma;%farcsec" % beamNew
</source>
 
We will set the integration time to 100 seconds and the total integration time to 8 hours.  For a real observation, 10 seconds may be a more appropriate integration time.  However, we can decrease the simulation CPU time by increasing the integration interval (See [[Etime study]] for details]]).  Whenever simulating a new data set, we recommend beginning with a large integration time for initial testing and then lowering the integration time to generate a more realistic measurement set.   


<source lang="python">
<source lang="python">
Line 225: Line 217:
</source>
</source>


We will add thermal noise to the simulated data by setting parameter '''thermalnoise''' to <tt>tsys-atm</tt>.  We will keep the default values of precipitable water vapor (1 mm) and ambient temperature (269 K).
=== Noise ===
 
We will add thermal noise to the simulated data by setting parameter '''thermalnoise''' to <tt>tsys-atm</tt>.  We keep the default values of precipitable water vapor (1 mm) and ambient temperature (269 K).


<source lang="python">
<source lang="python">
Line 232: Line 226:
</source>
</source>


Finally, we set '''graphics''' to <tt>"both"</tt> to send plots to the screen and disk; we set '''verbose''' to <tt>True</tt> to print additional information to the logger and terminal; and we set '''overwrite''' to <tt>True</tt> to overwrite existing files in the '''project''' subdirectory.   
=== <tt>sim_observe</tt> Execution ===
Now we are ready to run <tt>sim_observe</tt>.
 
Finally, we set '''graphics''' to <tt>"both"</tt> to send plots to the screen and disk; we set '''verbose''' to <tt>True</tt> to print additional information to the logger and terminal; and we set '''overwrite''' to <tt>True</tt> to overwrite existing files in the '''project''' subdirectory.  Now we are ready to run <tt>sim_observe</tt>.


<source lang="python">
<source lang="python">
Line 242: Line 237:
sim_observe()
sim_observe()
</source>
</source>
=== Aside: Antenna Configurations ===
[[file:Beamsummary.png|thumb|ALMA synthetic beam size as a function of array configuration number]]
In this tutorial the selection of the {{ALMA}} antenna configuration is automated to produce a decent cube of M51 as viewed at z = 0.1, but, for other simulations (or this one, for that matter), it will be worth playing with the configurations, or perhaps evaluating the possibility of detections in CSV or early science.  CASA comes with stock antenna configurations in the directory <tt>$CASAPATH/data/alma/simmos/</tt>.  Use the Python <tt>os</tt> module to resolve the CASA path.  For example:
CASA <37>: os.getenv('CASAPATH').split(' ')[0] + '/data/alma/simmos/'
  Out[37]: '/usr/lib64/casapy/stable/data/alma/simmos/'
CASA <38>: ls /usr/lib64/casapy/stable/data/alma/simmos/
aca.i.cfg                alma.out08.cfg  alma.out20.cfg  carma.d.cfg        vla.a.cfg
aca.ns.cfg                alma.out09.cfg  alma.out21.cfg  carma.e.cfg        vla.b.cfg
aca.tp.cfg                alma.out10.cfg  alma.out22.cfg  meerkat.cfg        vla.bna.cfg
alma.cycle0.compact.cfg  alma.out11.cfg  alma.out23.cfg  pdbi-a.cfg          vla.c.cfg
alma.cycle0.extended.cfg  alma.out12.cfg  alma.out24.cfg  pdbi-b.cfg          vla.cnb.cfg
alma.out01.cfg            alma.out13.cfg  alma.out25.cfg  pdbi-c.cfg          vla.d.cfg
alma.out02.cfg            alma.out14.cfg  alma.out26.cfg  pdbi-d.cfg          vla.dnc.cfg
alma.out03.cfg            alma.out15.cfg  alma.out27.cfg  sma.compact.cfg    WSRT.cfg
alma.out04.cfg            alma.out16.cfg  alma.out28.cfg  sma.compact.n.cfg
alma.out05.cfg            alma.out17.cfg  carma.a.cfg    sma.extended.cfg
alma.out06.cfg            alma.out18.cfg  carma.b.cfg    sma.subcompact.cfg
alma.out07.cfg            alma.out19.cfg  carma.c.cfg    sma.vextended.cfg
The version of CASA used to generate this example contains 28 50-element array configurations for ALMA, compact and extended configurations for ALMA cycle 0, as well as configuration files for [http://almascience.nao.ac.jp/ ACA], {{EVLA}}, [http://www.cfa.harvard.edu/sma/ SMA], and other interferometers.


=== <tt>sim_observe</tt> Output ===
=== <tt>sim_observe</tt> Output ===


Outputs from <tt>sim_observe</tt> will be written to a directory specified by the '''project''' parameter, in this case "M51-ATZ-p1".  <tt>sim_observe</tt> will generate these files:
[[File:M51-ATZ-p1.alma_0.360000arcsec.observe.png|thumb|right|Summary plot output by sim_observe (M51-ATZ-p1.alma_0.360000arcsec.observe.png).]]
 
Outputs from <tt>sim_observe</tt> will be written to a directory specified by the '''project''' parameter, in this case <tt>M51-AtZ-p1</tt>.  <tt>sim_observe</tt> will generate these files:


{| border="1"
{| border="1"
Line 251: Line 273:
! Description
! Description
|-
|-
| M51-ATZ-p1.alma_0.360000arcsec.ms
| M51-AtZ-p1.alma_0.360000arcsec.ms
| Model measurement set ''sans'' noise
| Model measurement set ''sans'' noise
|-  
|-  
| M51-ATZ-p1.alma_0.360000arcsec.noisy.ms
| M51-AtZ-p1.alma_0.360000arcsec.noisy.ms
| Model measurement set with thermal noise
| Model measurement set with thermal noise
|-
|-
| M51-ATZ-p1.alma_0.360000arcsec.quick.psf
| M51-AtZ-p1.alma_0.360000arcsec.quick.psf
| CASA image of the point spread function
| CASA image of the point spread function
|-
|-
| M51-ATZ-p1.alma_0.360000arcsec.skymodel
| M51-AtZ-p1.alma_0.360000arcsec.skymodel
| CASA image of the input sky model rescaled according to the skymodel sub-parameters
| CASA image of the input sky model rescaled according to the skymodel sub-parameters
|-
|-
| M51-ATZ-p1.alma_0.360000arcsec.skymodel.flat
| M51-AtZ-p1.alma_0.360000arcsec.skymodel.flat
| flattened CASA image of the input sky model rescaled
| flattened CASA image of the input sky model rescaled
|-  
|-  
| M51-ATZ-p1.alma_0.360000arcsec.skymodel.png
| M51-AtZ-p1.alma_0.360000arcsec.skymodel.png
| flattened PNG image of the input sky model rescaled and overlaid with the mosaic pattern specified by the setpointings sub-parameters
| flattened PNG image of the input sky model rescaled and overlaid with the mosaic pattern specified by the setpointings sub-parameters
|-
|-
| M51-ATZ-p1.alma_0.360000arcsec.observe.png
| M51-AtZ-p1.alma_0.360000arcsec.observe.png
| 2x2 PNG summary plot showing: source elevation vs. time, antenna position, uv coverage, and the point spread function
| 2x2 PNG summary plot showing: source elevation vs. time, antenna position, uv coverage, and the point spread function
|-
|-
| M51-ATZ-p1.alma_0.360000arcsec.ptg.txt
| M51-AtZ-p1.alma_0.360000arcsec.ptg.txt
| ASCII text listing of mosaic pointings
| ASCII text listing of mosaic pointings
|}
|}


And there are plenty of other auxiliary files. The image at the right can be created by running the task viewer and loading the appropriate image file (e.g. viewer(infile = 'M51-ATZ-p1/M51-ATZ-p1.image', channel = 16)).
=== <tt>sim_analyze</tt> Execution ===
 
After running <tt>sim_observe</tt> and producing a simulated measurement set, the task <tt>sim_analyze</tt> will quickly and easily image the simulated data and produce diagnostic plots and images. (Remember that <tt>sim_analyze</tt> is a convenience task. You can also use other CASA tasks, like <tt>clean</tt>, directly on the simulated data.)


=== <tt>sim_analyze</tt> ===
To setup the <tt>sim_analyze</tt> parameters we will keep the values of '''project''' and '''verbose''' that we set for <tt>sim_observe</tt>.  We set '''vis''' to image the noisy measurement set by concatenating '''project''', '''antennalist''' with ";" replaced by "_", and ".noisy.ms".  We set '''modelimage = cubeName''' so <tt>sim_analyze</tt> will use the model image when cleaning the simulated data.  We then set the clean threshold to 1 mJy.


<source lang="python">
<source lang="python">
# in CASA
# in CASA
modelimage = cubeName
modelimage = cubeName
# make simulated images/cubes using noisy data
vis = project +'.'+ antennalist.replace(';','_') +'.noisy.ms'
image=True # SIM_ANALYZE
vis = '$project.noisy.ms' # SIM_ANALYZE
cell = "%farcsec" % pixelSize
imsize = [imSizePix,imSizePix]
threshold = "1.0mJy"
threshold = "1.0mJy"
weighting = "natural"
stokes = 'I'
verbose = True
graphics="both"
overwrite = True
</source>
</source>


==== Other Antenna Configurations ====
We now set '''cell''' to the pixel size calculated above.


<source lang="python">
cell = "%farcsec" % pixelSize 
</source>


[[file:Beamsummary.png|thumb|ALMA synthetic beam size as a function of array configuration number]]
Now we have both the image size and pixel size in arcseconds, but the <tt>sim_analyze</tt> parameter '''imsize''' wants the ratio: the number of pixels along the RA or Dec axis. To keep the image from becoming too small, set the minimum image size to be 256 pixels.
Finally, we need to know which ALMA configuration number based on the desired angular resolution. [[simdata CASA 3.2|Simdata]] makes this easy by allowing users to specify the desired angular resolution in the parameter '''antennalist'''.


This tutorial is somewhat automated to produce a decent cube of M51 as viewed at z = 0.1. The selection of the {{ALMA}} antenna configuration is automated, but, for other simulations (or this one, for that matter), it will be worth playing with the configurations, or perhaps evaluating the possibility of detections in CSV or early science.
<source lang="python">
# in CASA
imSizePix = int(round(imSize / pixelSize))
if imSizePix < 256: imSizePix = 256
imsize = [imSizePix,imSizePix]
</source>


Notice from the [[simdata CASA 3.2|simdata]] inputs that CASA comes with stock antenna configurations in the directory $CASAPATH/data/alma/simmos/ (the python task '''os.getenv''' is used to look up CASAPATH automatically). For CASA 3.2, here is a list of included ALMA antenna configurations.
Now we are ready to execute the task.


{| border = 1
<source lang="python">
|+ '''ALMA Configuration Files'''
sim_analyze()
|-
</source>
| alma.cycle0.compact.cfg
|-
| alma.cycle0.extended.cfg
|-
| alma.out01.cfg
|-
| alma.out02.cfg
|-
| alma.out03.cfg
|-
| ...
|-
| alma.out27.cfg
|-
| alma.out28.cfg
|}
 
There are also configuration data for the [http://www.narrabri.atnf.csiro.au/ ACA], {{EVLA}}, and [http://www.cfa.harvard.edu/sma/ SMA].
 
=== Take a Break ===
 
If you have got this far, you've earned it. [[simdata CASA 3.2|Simdata]] will be running for a while, and coffee sure sounds good right now.


=== Results ===
=== <tt>sim_analyze</tt> Output ===


[[File:M51sim-09.png|thumb|Channel 16 of the M51 at z=0.1 simulation. The rms noise is about 0.25 mJy/beam, and the peak flux density on this channel is about 2 mJy/beam.]]
[[File:M51-ATZ-p1.alma 0.360000arcsec.analysis.png|thumb|<tt>sim_analyze</tt> summary plot displayed in the CASA plotter.]]


Here is an inventory of some of the [[simdata]] products. They will appear in the directory "M51-ATZ-p1", a sub-directory of your working directory.
The <tt>sim_analyze</tt> output files will appear in the CASA plotter and in subdirectory <tt>M51-AtZ-p1</tt>.


{| border="1"
{| border="1"
Line 342: Line 342:
! Description
! Description
|-
|-
| M51-ATZ-p1.ms
| M51-AtZ-p1.alma_0.360000arcsec.noisy.image
| Model measurement set ''sans'' noise
| CLEAN-deconvolved cube of M51-AtZ-p1.alma_0.360000arcsec.noisy.ms
|-
| M51-AtZ-p1.alma_0.360000arcsec.noisy.residual
| Residual image after cleaning
|-  
|-  
| M51-ATZ-p1.noisy.ms
| M51-AtZ-p1.alma_0.360000arcsec.skymodel.flat.regrid
| Model measurement set with thermal noise
| Input model image regridded to match the output image
|-
|-
| M51-ATZ-p1.image
| M51-AtZ-p1.alma_0.360000arcsec.skymodel.flat.regrid.conv
| CLEAN-deconvolved cube of M51-ATZ-p1.noisy.ms
| Input model image regridded to match the output image and convolved with the output synthesized beam
|}
|}


And there are plenty of other auxiliary files. The image at the right can be created by running the task viewer and loading the appropriate image file (e.g. viewer(infile = 'M51-ATZ-p1/M51-ATZ-p1.image', channel = 16)).
[[File:M51sim-09.png|thumb|Channel 16 of the M51 at z=0.1 simulation. The rms noise is about 0.25 mJy/beam, and the peak flux density on this channel is about 2 mJy/beam.]]
 
Additional (or fewer) outputs can be produced by <tt>sim_analyze</tt> by setting the '''analyze''' subparameters. The image at the right can be displayed by running the task viewer and loading the appropriate image file, e.g. <tt>viewer(infile = 'M51-AtZ-p1/M51-AtZ-p1.image', channel = 16)</tt>.


The rest frequency will have been lost in the simulation, and it's worth restoring to the header.
The rest frequency will have been lost in the simulation, and it's worth restoring to the header.


<source lang="python">
<source lang="python">
imhead(imagename="M51-ATZ-p1/M51-ATZ-p1.image", mode="put", hdkey="restfreq",  
# in CASA
      hdvalue="115.27120180GHz")
imhead(imagename="M51-AtZ-p1/M51-AtZ-p1.alma_0.360000arcsec.noisy.image",  
      mode="put", hdkey="restfreq", hdvalue="115.27120180GHz")
</source>
</source>


The simulated data cube can be analyzed just like any other CASA image -- examples are given in the [[Imaging a Mosaicked Spectral Line Dataset|CARMA tutorial]] and the [[NGC 5921: red-shifted HI emission#Cube Moments|VLA 21cm tutorial]].
The simulated data cube can be analyzed just like any other CASA image -- examples are given in the [[Imaging a Mosaicked Spectral Line Dataset|CARMA tutorial]] and the [[NGC 5921: red-shifted HI emission#Cube Moments|VLA 21cm tutorial]].


==== Moment Maps ====
=== Moment Maps ===
 
[[File:M51simnew.png|thumb|Moment maps of the M51 CO 1-0 at z=0.1 simulation. ''Left'': Moment 0 (integrated intensity) map. ''Right'': Moment 1 (velocity field) map. Note that the rotation sense has been flipped, exactly [[#Adjusting the Frequency Axis | as expected]].]]
[[File:M51simnew.png|thumb|Moment maps of the M51 CO 1-0 at z=0.1 simulation. ''Left'': Moment 0 (integrated intensity) map. ''Right'': Moment 1 (velocity field) map. Note that the rotation sense has been flipped, exactly [[#Adjusting the Frequency Axis | as expected]].]]


Line 369: Line 376:


<source lang="python">
<source lang="python">
immoments(imagename='M51-ATZ-p1/M51-ATZ-p1.image',moments=[0,1],axis='spectral',
# Must remove this file before immoments is run, if this has been run previously in the directory
          excludepix=[-100,0.0006],outfile='M51-ATZ-p1/M51-ATZ-p1.moments')
os.system('rm -rf M51-AtZ-p1/M51-AtZ-p1.alma_0.360000arcsec.noisy.moments.*')
immoments(imagename='M51-AtZ-p1/M51-AtZ-p1.alma_0.360000arcsec.noisy.image',
          moments=[0,1],axis='spectral', excludepix=[-100,0.0006],
          outfile='M51-AtZ-p1/M51-AtZ-p1.alma_0.360000arcsec.noisy.moments')
</source>
</source>


The results are shown at right.
The results are shown at right.


You can view these results by running the task viewer, then selecting the moments maps found in the directory M51-ATZ-p1 (M51-ATZ-p1.moments.integrated and M51-ATZ-p1.moments.weighted_coord)
You can view these results by running the task viewer, then selecting the moments maps found in the directory <tt>M51-AtZ-p1</tt>.


=== Pushing M51 out to z = 0.3 ===
=== Pushing M51 out to z = 0.3 ===


We'll leave it as an exercise how to tune this simulation to push M51 all the way out to z = 0.3 (luminosity distance = 1540 Mpc). [[simdata CASA 3.2|Simdata]] produces a reasonable detection in 8 hrs integration, but you'll need to consider more carefully the antenna configuration needed to produce the requisite sensitivity. To conclude, here are the moment maps for the z = 0.3 simulation.
We'll leave as an exercise for the reader the tuning of this simulation to push M51 all the way out to z = 0.3 (luminosity distance = 1540 Mpc). <tt>sim_observe</tt> produces a reasonable detection in 8 hrs integration, but you'll need to consider more carefully the antenna configuration needed to produce the requisite sensitivity. To conclude, here are the moment maps for the z = 0.3 simulation.


[[File:M51sim78.png|400px|Simulation of M51 CO 1-0 at z = 0.3. ''Left'': integrated intensity map. ''Right'': (marginally resolved) velocity field.]]
[[File:M51sim78.png|400px|Simulation of M51 CO 1-0 at z = 0.3. ''Left'': integrated intensity map. ''Right'': (marginally resolved) velocity field.]]


{{Simulations Intro}}
{{Simulations Intro}}
 
{{Checked 3.3.0}}
{{Checked 3.2.0}}
 
--[[User:Jgallimo|Jack Gallimore]] 15:55, 30 April 2010 (UTC), updated for CASA 3.2 by A. Kimball, 28th April 2011

Latest revision as of 11:15, 22 November 2011


Simulating Observations in CASA

This guide is applicable to CASA version 3.3. For older versions of CASA, see M51 at z = 0.1 and z = 0.3 (CASA 3.2).

Follow this process, Extracting scripts from these tutorials, to extract a casapy script from this web page. Also notice that this page uses global Python variables (e.g. fluxScale), not to be confused with task inputs (e.g. inbright).

Overview

CO 1-0 Moment maps of the original BIMA SONG measurements of M51. Left: Moment 0 (integrated intensity). Right: Moment 1 (velocity field).

This tutorial presents a simulation of ALMA observations of the well-known galaxy M51 in the CO 1-0 transition. This galaxy is located relatively nearby (luminosity distance = 9 Mpc), but, in this simulation, we will model how it would appear at redshift z = 0.1 (luminosity distance = 460 Mpc).

The goal of this tutorial is to provide a complete run-through of a relatively simple simulation. The simulation will include effects of (u, v) sampling of a 50-antenna ALMA, the primary beam of the ALMA, and thermal noise. Neither calibration overheads nor errors are included, and so this simulation should be viewed as optimistic. The CASA task sim_observe will generate simulated visibility measurements from an input model. The task sim_analyze will image the visibilities and generate diagnostic plots.

For this tutorial, we'll use the BIMA SONG (Helfer et al. 2003) observations of M51 as the basis for the model. Grab the file NGC5194.bima12m.cm.fits.gz and uncompress it in a working directory.

# in bash (or other unix shell)
gunzip NGC5194.bima.12m.cm.fits.gz

Initializing the Data and Tasks

Task sim_observe takes an input model in the form of a CASA image. So, we begin by converting the BIMA SONG FITS cube into a CASA image using task importfits. Here we invoke importfits as a function, where the task inputs are set as function input parameters. For convenience, we store the name of the resulting image into the Python variable cubeName. The task will produce the new image (which is actually a directory) in the current working directory.

# in CASA
cubeName = 'm51-song.image'
# If m51-song.image exists, it needs to be removed for importfits to work as set
os.system('rm -rf '+cubeName)
importfits(fitsimage='NGC5194.bima12m.cm.fits', imagename=cubeName)

We initialize the simulation tasks by setting their parameters to the default values. We then make the new CASA image the input model to sim_observe by setting parameter skymodel.

# in CASA
default("sim_observe")
default("sim_analyze")
skymodel = cubeName

Noise in the Input Model

The input model is actually an observation and, as such, certainly contains noise. We're OK as long as the noise in the model falls below the expected thermal noise of the ALMA observation. Here's a quick, back-of-the-envelope calculation.

The noise on the BIMA SONG channels measures roughly 0.1 Jy/beam, which scales to 0.04 mJy/beam at z = 0.1 (to within factors of powers of 1+z for cosmology and root-small-integer for the to-be-degraded beam; see below). For comparison, the expected thermal noise of the ALMA observation is 0.1 mJy / beam (using the ALMA Sensitivity Calculator, assuming 8 hrs integration, and, to match roughly the BIMA SONG observation, a 4 MHz channel width), a factor of 2.5 greater than the anticipated contribution of noise from the model. Added in quadrature, the noise components total 0.108 mJy, and so we can expect that the input model noise degrades the simulation noise by about 8%.

Cosmology Calculations

Next we'll set up some python globals to handle the scaling of the model coordinates and flux densities appropriate for new redshift. We'll primarily need the angular size and luminosity distances for a given cosmology. To keep things simple, we'll use Ned Wright's CosmoCalc (2006) with the default cosmology; redshifts were collected from NED.

# in CASA
# Distinguish between z_lsrk, which sets the observed frequency for scaling, 
# and z_cmb, which is needed to get cosmological distances.
z_old_cmb = 0.002122 # CMB-referenced z for cosmological distances from NED
z_old_lsrk = 0.001544 # from NED
z_new = 0.1
#
# angular size distances from CosmoCalc
da_old = 9.0
da_new = 375.9
#
# luminosity distances from CosmoCalc
dl_old = 8.937
dl_new = 454.8

We will use the convention where old refers to M51 as observed at its proper redshift, and new refers to the new, higher redshift for our model.

Preparing the Model

The next step is to scale the M51 data cube into a model cube appropriate for sim_observe. First, we establish the file naming convention by setting the task parameter project. Thinking ahead, we might want to rerun this simulation for other redshifts. So, to make this easy we will include the redshift in the project name: we'll write 0.1 as p1 (or "point one") and assign it to Python variable suffix. Then we'll use suffix to define project.

# in CASA
suffix = "-p1"
project = "M51-AtZ" + suffix

Flux Density Scaling

Our goal in this step is to set the inbright parameter of sim_observe, which will scale the data cube appropriately for its new luminosity distance. The task wants models in units of Jy / pixel, but the BIMA SONG cube has units of Jy / beam. We'll use the task imhead to get the BIMA SONG synthesized beam dimensions and we'll use the qa tool to convert the dimensions from radians to arcseconds.

# in CASA
# Get the major and minor axes of the model clean beam
bmaj = imhead(imagename=cubeName,mode='get',hdkey='beammajor')
bmin = imhead(imagename=cubeName,mode='get',hdkey='beamminor')
# Convert radians to 1-arcsecond pixels using the qa tool
bmaj = qa.convert(bmaj,'arcsec')['value']
bmin = qa.convert(bmin,'arcsec')['value']
# Gaussian beam conversion = beams / pixel
toJyPerPix = 1.0 / (1.1331 * bmaj * bmin)

Next, scale the flux for the new luminosity distance. Make the approximation that each pixel is a point source, and use the inverse square law to scale. We then set inbright to the newly scaled peak flux.

# in CASA
# Correct flux density for luminosity distance
fluxScale = (dl_old/dl_new)**2 * (1.0 + z_new) / (1.0 + z_old_cmb)
# Current peak flux (Jy / pixel)
peak = imstat(cubeName)['max'][0] * toJyPerPix
# Desired peak flux (using Python formatting convention)
inbright = "%fJy/pixel" % (peak*fluxScale)

Notice that there is an additional (1+z) correction because we are scaling flux densities rather than bolometric fluxes.

Angular Size Scaling

The sky coordinates axes of the model need to be adjusted (1) to place M51 in the southern hemisphere and (2) to correct for the new angular size distance. These tasks can be accomplished by the sim_observe parameters indirection (to change the location of M51 on the sky) and incell (to change the angular scale of an input pixel). To perform task (1), we'll just flip the sign of the declination using imhead. We'll use the qa tool to convert radians to sexagesimal units.

# in CASA
# For clarity, build up "indirection" string one term at a time.
indirection = "J2000 " # Epoch
# RA
crval1 = imhead(imagename=cubeName,mode='get',hdkey='crval1')['value']
indirection += qa.formxxx(str(crval1)+'rad',format='hms') + " " 
# Dec * -1
crval2 = imhead(imagename=cubeName,mode='get',hdkey='crval2')['value']
indirection += qa.formxxx('%frad' % (-1*float(crval2)),format='dms')

Next, adjust the pixel scale for the new angular size distance. To perform this adjustment, we'll use imhead with mode = "get" to read in the original pixel scale, then use the sim_observe input parameter incell to designate the new pixel size.

# in CASA
# Scale pixel size: imhead returns radians; convert to arcsec
cdelt2 = imhead(imagename=cubeName,mode='get',hdkey='cdelt2')['value']
oldCell = cdelt2 * 206265
# Scale for new angular size distance
newCell = oldCell * da_old / da_new 
# Format the new pixel size for input to sim_observe
incell = "%farcsec" % (newCell)

Adjusting the Frequency Axis

Changing the frequency axis of the model header is just a straightforward (1+z) correction. sim_observe adjusts the channelwidth using the inwidth parameter and the observing frequency using incenter.

Notice that the absolute value of the BIMA SONG channel width is used. The noise calculation in sim_observe needs positive channel widths, but the input cube is ordered in increasing velocity rather than increasing frequency. We could transpose the cube; on the other hand, the only penalty of changing the sign is that the sense of rotation will be flipped. Since this simulation is a simple detection experiment, the rotation sense is irrelevant, and so we'll take the easy way out.

# in CASA
# Move freq to z_new
oldFreq = imhead(imagename=cubeName,mode='get',hdkey='crval3')['value']
newFreq = oldFreq * (1.0 + z_old_lsrk) / (1.0 + z_new)
nchan = imstat("NGC5194.bima12m.cm.fits")['trc'][2]
# Adjust frequency channelwidth for new z
oldDnu = imhead(imagename=cubeName,mode='get',hdkey='cdelt3')['value']
newDnu = abs((1.0+z_old_lsrk) /(1.0+z_new)*oldDnu) # make channel widths positive
inwidth = "%fHz" % newDnu
# Specify the observing frequency at the center of the observing band:
incenter = "%fHz" % (newFreq + 0.5*nchan*newDnu)

Map Size and Mosaicking

The original BIMA SONG image is about 480 arcseconds across; scale this image size to the new redshift and assign to parameter mapsize. We will also use this value later.

# in CASA
imSize = 480.0 * (da_old / da_new) # arcseconds
mapsize = "%farcsec" % imSize

For relatively high redshifts, there should be no need to mosaic the observations. We'll nevertheless allow for mosaicking in case we want to repeat the simulation for lower redshift. sim_observe needs the spacing between pointings in the mosaic; we'll require pointings spaced by half of the primary beam. The ALMA primary beam is 17 arcseconds at 300 GHz.

# in CASA
primaryBeam = 17.0 * (300e9 / newFreq) # arcseconds
pointingspacing = "%farcsec" % (primaryBeam / 2.0)

Array Configuration

We want to set antennalist such that sim_observe will decide on an appropriate ALMA configuration based on the desired beam size. Thus, we first must estimate the desired synthesized beam size. We don't want the new beam to be so large so that we cannot resolve a rotation curve, but we also don't want it to be so small that we effectively resolve out the BIMA SONG data. The BIMA SONG beam was 5 arcsec, and so as a reasonable guess we'll adopt the equivalent of a 15" beam at its true distance (3 times coarser than the BIMA SONG beam) and then scale appropriately to z = 0.1.

Looking ahead to imaging our simulated data, it's important to sample the beam at least 3 times for stable deconvolution. We'll sample the beam by 4 times, rounded off to the nearest milliacsecond.

# in CASA
beamNew = 15.0 * (da_old / da_new)
pixelSize = round(beamNew * 1000.0 / 4.0) / 1000.0

Now, to guard against undersampling the beam as a result of rounding error, reset the desired beam to 4 times the pixel size.

# in CASA
beamNew = 4.0 * pixelSize

Now we can set parameter antennalist using the beam size (but see Aside: Antenna Configurations below).

# in CASA
antennalist = "alma;%farcsec" % beamNew

Integration and Total Time

We will set the integration time to 500 seconds and the total integration time to 8 hours. For a real observation, 10 seconds will be a more appropriate integration time. However, we can decrease the simulation CPU time by increasing the integration interval (See Etime study for details). Whenever simulating a new data set, we recommend beginning with a large integration time for initial testing and then lowering the integration time to generate a more realistic measurement set.

# in CASA
integration = '100s'
# 8 hours of total integration
totaltime = '28800s'

Noise

We will add thermal noise to the simulated data by setting parameter thermalnoise to tsys-atm. We keep the default values of precipitable water vapor (1 mm) and ambient temperature (269 K).

# in CASA
thermalnoise="tsys-atm"

sim_observe Execution

Finally, we set graphics to "both" to send plots to the screen and disk; we set verbose to True to print additional information to the logger and terminal; and we set overwrite to True to overwrite existing files in the project subdirectory. Now we are ready to run sim_observe.

# in CASA
graphics="both"
verbose = True
overwrite = True
sim_observe()

Aside: Antenna Configurations

ALMA synthetic beam size as a function of array configuration number

In this tutorial the selection of the ALMA antenna configuration is automated to produce a decent cube of M51 as viewed at z = 0.1, but, for other simulations (or this one, for that matter), it will be worth playing with the configurations, or perhaps evaluating the possibility of detections in CSV or early science. CASA comes with stock antenna configurations in the directory $CASAPATH/data/alma/simmos/. Use the Python os module to resolve the CASA path. For example:

CASA <37>: os.getenv('CASAPATH').split(' ')[0] + '/data/alma/simmos/'
  Out[37]: '/usr/lib64/casapy/stable/data/alma/simmos/'
CASA <38>: ls /usr/lib64/casapy/stable/data/alma/simmos/
aca.i.cfg                 alma.out08.cfg  alma.out20.cfg  carma.d.cfg         vla.a.cfg
aca.ns.cfg                alma.out09.cfg  alma.out21.cfg  carma.e.cfg         vla.b.cfg
aca.tp.cfg                alma.out10.cfg  alma.out22.cfg  meerkat.cfg         vla.bna.cfg
alma.cycle0.compact.cfg   alma.out11.cfg  alma.out23.cfg  pdbi-a.cfg          vla.c.cfg
alma.cycle0.extended.cfg  alma.out12.cfg  alma.out24.cfg  pdbi-b.cfg          vla.cnb.cfg
alma.out01.cfg            alma.out13.cfg  alma.out25.cfg  pdbi-c.cfg          vla.d.cfg
alma.out02.cfg            alma.out14.cfg  alma.out26.cfg  pdbi-d.cfg          vla.dnc.cfg
alma.out03.cfg            alma.out15.cfg  alma.out27.cfg  sma.compact.cfg     WSRT.cfg
alma.out04.cfg            alma.out16.cfg  alma.out28.cfg  sma.compact.n.cfg
alma.out05.cfg            alma.out17.cfg  carma.a.cfg     sma.extended.cfg
alma.out06.cfg            alma.out18.cfg  carma.b.cfg     sma.subcompact.cfg
alma.out07.cfg            alma.out19.cfg  carma.c.cfg     sma.vextended.cfg 

The version of CASA used to generate this example contains 28 50-element array configurations for ALMA, compact and extended configurations for ALMA cycle 0, as well as configuration files for ACA, EVLA, SMA, and other interferometers.

sim_observe Output

Summary plot output by sim_observe (M51-ATZ-p1.alma_0.360000arcsec.observe.png).

Outputs from sim_observe will be written to a directory specified by the project parameter, in this case M51-AtZ-p1. sim_observe will generate these files:

Filename Description
M51-AtZ-p1.alma_0.360000arcsec.ms Model measurement set sans noise
M51-AtZ-p1.alma_0.360000arcsec.noisy.ms Model measurement set with thermal noise
M51-AtZ-p1.alma_0.360000arcsec.quick.psf CASA image of the point spread function
M51-AtZ-p1.alma_0.360000arcsec.skymodel CASA image of the input sky model rescaled according to the skymodel sub-parameters
M51-AtZ-p1.alma_0.360000arcsec.skymodel.flat flattened CASA image of the input sky model rescaled
M51-AtZ-p1.alma_0.360000arcsec.skymodel.png flattened PNG image of the input sky model rescaled and overlaid with the mosaic pattern specified by the setpointings sub-parameters
M51-AtZ-p1.alma_0.360000arcsec.observe.png 2x2 PNG summary plot showing: source elevation vs. time, antenna position, uv coverage, and the point spread function
M51-AtZ-p1.alma_0.360000arcsec.ptg.txt ASCII text listing of mosaic pointings

sim_analyze Execution

After running sim_observe and producing a simulated measurement set, the task sim_analyze will quickly and easily image the simulated data and produce diagnostic plots and images. (Remember that sim_analyze is a convenience task. You can also use other CASA tasks, like clean, directly on the simulated data.)

To setup the sim_analyze parameters we will keep the values of project and verbose that we set for sim_observe. We set vis to image the noisy measurement set by concatenating project, antennalist with ";" replaced by "_", and ".noisy.ms". We set modelimage = cubeName so sim_analyze will use the model image when cleaning the simulated data. We then set the clean threshold to 1 mJy.

# in CASA
modelimage = cubeName
vis = project +'.'+ antennalist.replace(';','_') +'.noisy.ms'
threshold = "1.0mJy"

We now set cell to the pixel size calculated above.

cell = "%farcsec" % pixelSize

Now we have both the image size and pixel size in arcseconds, but the sim_analyze parameter imsize wants the ratio: the number of pixels along the RA or Dec axis. To keep the image from becoming too small, set the minimum image size to be 256 pixels.

# in CASA
imSizePix = int(round(imSize / pixelSize))
if imSizePix < 256: imSizePix = 256
imsize = [imSizePix,imSizePix]

Now we are ready to execute the task.

sim_analyze()

sim_analyze Output

sim_analyze summary plot displayed in the CASA plotter.

The sim_analyze output files will appear in the CASA plotter and in subdirectory M51-AtZ-p1.

Filename Description
M51-AtZ-p1.alma_0.360000arcsec.noisy.image CLEAN-deconvolved cube of M51-AtZ-p1.alma_0.360000arcsec.noisy.ms
M51-AtZ-p1.alma_0.360000arcsec.noisy.residual Residual image after cleaning
M51-AtZ-p1.alma_0.360000arcsec.skymodel.flat.regrid Input model image regridded to match the output image
M51-AtZ-p1.alma_0.360000arcsec.skymodel.flat.regrid.conv Input model image regridded to match the output image and convolved with the output synthesized beam
Channel 16 of the M51 at z=0.1 simulation. The rms noise is about 0.25 mJy/beam, and the peak flux density on this channel is about 2 mJy/beam.

Additional (or fewer) outputs can be produced by sim_analyze by setting the analyze subparameters. The image at the right can be displayed by running the task viewer and loading the appropriate image file, e.g. viewer(infile = 'M51-AtZ-p1/M51-AtZ-p1.image', channel = 16).

The rest frequency will have been lost in the simulation, and it's worth restoring to the header.

# in CASA
imhead(imagename="M51-AtZ-p1/M51-AtZ-p1.alma_0.360000arcsec.noisy.image", 
       mode="put", hdkey="restfreq", hdvalue="115.27120180GHz")

The simulated data cube can be analyzed just like any other CASA image -- examples are given in the CARMA tutorial and the VLA 21cm tutorial.

Moment Maps

Moment maps of the M51 CO 1-0 at z=0.1 simulation. Left: Moment 0 (integrated intensity) map. Right: Moment 1 (velocity field) map. Note that the rotation sense has been flipped, exactly as expected.

Use immoments to calculate the integrated intensity and velocity field maps from the simulated cube. The excludepix option applies a 3σ cut.

# Must remove this file before immoments is run, if this has been run previously in the directory
os.system('rm -rf M51-AtZ-p1/M51-AtZ-p1.alma_0.360000arcsec.noisy.moments.*')
immoments(imagename='M51-AtZ-p1/M51-AtZ-p1.alma_0.360000arcsec.noisy.image',
          moments=[0,1],axis='spectral', excludepix=[-100,0.0006],
          outfile='M51-AtZ-p1/M51-AtZ-p1.alma_0.360000arcsec.noisy.moments')

The results are shown at right.

You can view these results by running the task viewer, then selecting the moments maps found in the directory M51-AtZ-p1.

Pushing M51 out to z = 0.3

We'll leave as an exercise for the reader the tuning of this simulation to push M51 all the way out to z = 0.3 (luminosity distance = 1540 Mpc). sim_observe produces a reasonable detection in 8 hrs integration, but you'll need to consider more carefully the antenna configuration needed to produce the requisite sensitivity. To conclude, here are the moment maps for the z = 0.3 simulation.

Simulation of M51 CO 1-0 at z = 0.3. Left: integrated intensity map. Right: (marginally resolved) velocity field.

Simulating Observations in CASA

Last checked on CASA Version 3.3.0.