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

From CASA Guides
Jump to navigationJump to search
(Undo revision 4879 by Akimball (talk))
 
(48 intermediate revisions by 6 users not shown)
Line 1: Line 1:
[[Category: Simulations]] [[Category:ALMA]]
[[Category: Simulations]] [[Category:ALMA]]
{{Simulations Intro 3.1}}
Remember that you need to follow this process: [[Extracting scripts from these tutorials]] to extract a casapy script from this web page.


=== Overview ===
=== Overview ===
Line 21: Line 25:
cubeName = 'm51-song'
cubeName = 'm51-song'
importfits(fitsimage='NGC5194.bima12m.cm.fits', imagename=cubeName)
importfits(fitsimage='NGC5194.bima12m.cm.fits', imagename=cubeName)
 
# initialize simdata
# initialize simdata2
default("simdata")
default("simdata2")
skymodel=cubeName
skymodel=cubeName
modifymodel=True
modifymodel=True
Line 32: Line 35:
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 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 [http://www.eso.org/sci/facilities/alma/observing/tools/etc/index.html 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%.
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 [https://almascience.nrao.edu/call-for-proposals/sensitivity-calculator 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 ===
=== Cosmology Calculations ===
Line 45: Line 48:
z_old_lsrk = 0.001544 # from NED
z_old_lsrk = 0.001544 # from NED
z_new = 0.1
z_new = 0.1
 
#
# angular size distances from CosmoCalc
# angular size distances from CosmoCalc
da_old = 9.0
da_old = 9.0
da_new = 375.9
da_new = 375.9
 
#
# luminosity distances from CosmoCalc
# luminosity distances from CosmoCalc
dl_old = 8.937
dl_old = 8.937
Line 59: Line 62:
=== Preparing the Model ===
=== Preparing the Model ===


The next step is to scale the M51 data cube into a model cube appropriate for [[simdata2]]. First, we'll set up some globals to establish some file naming conventions.
The next step is to scale the M51 data cube into a model cube appropriate for [[simdata CASA 3.1| simdata]]. First, we'll set up some globals to establish some file naming conventions.


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


Line 69: Line 72:
==== Flux Density Scaling ====
==== Flux Density Scaling ====


Our goal in this step is to set the '''inbright''' parameter of [[simdata2]], which will scale the data cube appropriately for its new luminosity distance.
Our goal in this step is to set the '''inbright''' parameter of [[simdata CASA 3.1 | simdata]], which will scale the data cube appropriately for its new luminosity distance.


[[simdata2|Simdata2]] wants models in units of Jy / pixel, but the BIMA SONG cube is in units of Jy / beam. That's an easy conversion.
[[simdata CASA 3.1|Simdata]] wants models in units of Jy / pixel, but the BIMA SONG cube is in units of Jy / beam. That's an easy conversion.


<source lang="python">
<source lang="python">
Line 77: Line 80:
bmaj = imhead(imagename=cubeName,mode='get',hdkey='beammajor')
bmaj = imhead(imagename=cubeName,mode='get',hdkey='beammajor')
bmin = imhead(imagename=cubeName,mode='get',hdkey='beamminor')
bmin = imhead(imagename=cubeName,mode='get',hdkey='beamminor')
 
#
# use [[qa]] tool to convert beam in radians to (1 arcsec) pixels:
# use qa tool to convert beam in radians to (1 arcsec) pixels:
bmaj = qa.convert(bmaj,'arcsec')['value']
bmaj = qa.convert(bmaj,'arcsec')['value']
bmin = qa.convert(bmin,'arcsec')['value']
bmin = qa.convert(bmin,'arcsec')['value']
 
#
toJyPerPix = 1.0 / (1.1331 * bmaj * bmin) # gaussian beam conversion = beams / pixel
toJyPerPix = 1.0 / (1.1331 * bmaj * bmin) # gaussian beam conversion = beams / pixel
</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. [[simdata2|Simdata2]] can scale the peak surface brightness using the parameter '''inbright'''.
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. [[simdata CASA 3.1|Simdata]] can scale the peak surface brightness using the parameter '''inbright'''.


<source lang="python">
<source lang="python">
Line 94: Line 97:
# desired peak flux in Jy / pixel:
# desired peak flux in Jy / pixel:
inbright = "%fJy/pixel" % (peak*fluxScale) # use python formatting convention
inbright = "%fJy/pixel" % (peak*fluxScale) # use python formatting convention
# inbright is a simdata2 parameter
# inbright is a simdata parameter
</source>
</source>


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


==== 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 [[simdata2]] parameters '''indirection''' (change the location of M51 on the sky) and '''incell''' (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 [[simdata CASA 3.1| simdata]] parameters '''indirection''' (change the location of M51 on the sky) and '''incell''' (change the angular scale of an input pixel).
 
a
To perform task (1), we'll just flip the sign of the declination using [[imhead]]. Notice the use of the [[qa]] tool to convert radians to sexagesimal.
To perform task (1), we'll just flip the sign of the declination using [[imhead]]. Notice the use of the [http://casa.nrao.edu/docs/casaref/quanta-Module.html#x409-4100001.4 qa] tool to convert radians to sexagesimal.


<source lang="python">
<source lang="python">
Line 127: Line 130:
newCell = oldCell * da_old / da_new  
newCell = oldCell * da_old / da_new  
# Format the new pixel size for input to simdata
# Format the new pixel size for input to simdata
incell = "%farcsec" % (newCell) # parameter for simdata2
incell = "%farcsec" % (newCell) # parameter for simdata
</source>
</source>


==== Adjusting the Frequency Axis ====
==== Adjusting the Frequency Axis ====


Changing the frequency axis of the model header is just a straightforward (1+''z'') correction. [[simdata2|Simdata2]] adjusts the channelwidth using the '''inwidth''' parameter and the observing frequency using '''incenter'''.  
Changing the frequency axis of the model header is just a straightforward (1+''z'') correction. [[simdata CASA 3.1|Simdata]] 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 [[simdata CASA 3.1|Simdata]] 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.


<source lang="python">
<source lang="python">
Line 139: Line 144:
newFreq = oldFreq * (1.0 + z_old_lsrk) / (1.0 + z_new)
newFreq = oldFreq * (1.0 + z_old_lsrk) / (1.0 + z_new)
nchan = imstat("NGC5194.bima12m.cm.fits")['trc'][2]
nchan = imstat("NGC5194.bima12m.cm.fits")['trc'][2]
 
#
# Adjust frequency channelwidth for new z
# Adjust frequency channelwidth for new z
oldDnu = float(imhead(imagename=cubeName,mode='get',hdkey='cdelt3')['value'])  # Hz
oldDnu = float(imhead(imagename=cubeName,mode='get',hdkey='cdelt3')['value'])  # Hz
newDnu = abs((1.0+z_old_lsrk) /(1.0+z_new)*oldDnu)
newDnu = abs((1.0+z_old_lsrk) /(1.0+z_new)*oldDnu) # need positive channel widths to ensure noise calc goes well
inwidth = "%fHz" % newDnu # parameter for simdata
inwidth = "%fHz" % newDnu # parameter for simdata
 
#
# Specify the observing frequency at the center of the observing band:
# Specify the observing frequency at the center of the observing band:
incenter = "%fHz" % (newFreq + 0.5*nchan*newDnu)
incenter = "%fHz" % (newFreq + 0.5*nchan*newDnu)
</source>
</source>


=== Simdata2 ===
=== Simdata ===


The CASA task [[simdata2]] will monolithically simulate an ALMA observation, produce measurement sets with and without thermal noise, and finally produce a CLEANed image cube based on the simulated observation. Details can be found in the [[Simulating Observations in CASA]] tutorial.
The CASA task [[simdata]] will monolithically simulate an ALMA observation, produce measurement sets with and without thermal noise, and finally produce a CLEANed image cube based on the simulated observation. Details can be found in the [[Simulating Observations in CASA]] tutorial.


The original BIMA SONG image is about 480 arcseconds across; scale this image size to the new redshift.
The original BIMA SONG image is about 480 arcseconds across; scale this image size to the new redshift.
Line 160: Line 165:
</source>
</source>


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. [[simdata2|Simdata2]] needs the spacing between pointings in the mosaic; we'll require pointings spaced by half of the primary beam.  
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. [[simdata|Simdata]] needs the spacing between pointings in the mosaic; we'll require pointings spaced by half of the primary beam.  


<source lang="python">
<source lang="python">
# mosaicking info
# mosaicking info
primaryBeam = 17.0 * (300. / thisFreq1) # in arcseconds; ALMA primary beam = 17 arcsec at 300 GHz
primaryBeam = 17.0 * (300e9 / newFreq) # in arcseconds; ALMA primary beam = 17 arcsec at 300 GHz
pointingSpacing = primaryBeam / 2.0 # in arcseconds
pointingspacing = "%farcsec" % (primaryBeam / 2.0)
mosaicSize = "%farcsec" % imSize # how big to make the mosaic
mapsize = "%farcsec" % imSize # how big to make the mosaic
</source>
</source>


Line 183: Line 188:


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


<source lang="python">
<source lang="python">
Line 189: Line 193:
</source>
</source>


Now we know both the image size and pixel size in arcseconds, but [[simdata2]] 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 know both the image size and pixel size in arcseconds, but [[simdata CASA 3.1| simdata]] 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.


<source lang="python">
<source lang="python">
Line 196: Line 200:
</source>
</source>


[[file:Beamsummary.png|thumb|ALMA synthetic beam size as a function of array configuration number]]
Let [[simdata]] 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).
Finally, we need to know which ALMA configuration number based on the desired angular resolution. [[simdata2|Simdata2]] makes this easy by allowing users to specify the desired angular resolution in the parameter '''antennalist'''.


<source lang="python">
<source lang="python">
Line 204: Line 207:
</source>
</source>


Now we have enough information to run [[simdata2]], and hopefully some of the python global variables that were defined above will start to make sense. The following CASA and python commands set up the remaining parameters for the [[simdata2]] task.  
Now we have enough information to run [[simdata CASA 3.1| simdata]], and hopefully some of the python global variables that were defined above will start to make sense. The following CASA and python commands set up the remaining parameters for the [[simdata]] task.  


<source lang="python">
<source lang="python">
Line 211: Line 214:
# 10s is usually safe for large mosaics, but will require more memory
# 10s is usually safe for large mosaics, but will require more memory
totaltime = '28800s' # 8 hr integration
totaltime = '28800s' # 8 hr integration
</source>
<font color="red">If you are doing this as a tutorial, we suggest that you consider increasing the integration time so that your learning is not slowed by processing time</font> See [[Etime study]] for details.


<source lang="python">
# make simulated images/cubes
# make simulated images/cubes
image=True
image=True
Line 221: Line 228:
weighting = "natural"
weighting = "natural"
stokes = 'I'
stokes = 'I'
 
#
verbose = True
verbose = True
graphics="both"
graphics="both"
overwrite = True
overwrite = True
 
#
inp("simdata2")
inp("simdata")
simdata2()
simdata()
</source>
</source>




==== Other Antenna Configurations ====
==== Other Antenna Configurations ====
[[file:Beamsummary.png|thumb|ALMA synthetic beam size as a function of array configuration number]]
Finally, we need to know which ALMA configuration number based on the desired angular resolution. [[simdata CASA 3.1|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.
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.


Notice from the [[simdata2]] 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.0.1, here is a list of included antenna configurations.
Notice from the [[simdata CASA 3.1|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.0.1 and 3.1, here is a list of included antenna configurations.


{| border = 1
{| border = 1
Line 265: Line 276:
=== Take a Break ===
=== Take a Break ===


If you have got this far, you've earned it. [[simdata2|Simdata2]] will be running for a while, and coffee sure sounds good right now.
If you have got this far, you've earned it. [[simdata CASA 3.1|Simdata]] will be running for a while, and coffee sure sounds good right now.


=== Results ===
=== Results ===


[[File:M51sim-01.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: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.]]


Here is an inventory of some of the [[simdata2]] products.  
Here is an inventory of some of the [[simdata]] products.  


{| border="1"
{| border="1"
Line 292: Line 303:


<source lang="python">
<source lang="python">
imhead(imagename="M51-ATZ-p1.clean.image", mode="put", hdkey="restfreq", hdvalue="115.27120180GHz")
imhead(imagename="M51-ATZ-p1.image", mode="put", hdkey="restfreq", hdvalue="115.27120180GHz")
</source>
</source>


Line 298: Line 309:


==== Moment Maps ====
==== Moment Maps ====
[[File:M51sim34.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.]]
[[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]].]]


Use [[immoments]] to calculate the integrated intensity and velocity field maps from the simulated cube. The '''excludepix''' option applies a 3&sigma; cut.
Use [[immoments]] to calculate the integrated intensity and velocity field maps from the simulated cube. The '''excludepix''' option applies a 3&sigma; cut.


<source lang="python">
<source lang="python">
immoments(imagename='M51-ATZ-p1.clean.image',moments=[0,1],axis='spectral',
immoments(imagename='M51-ATZ-p1.image',moments=[0,1],axis='spectral',
           excludepix=[-100,0.0006],outfile='M51-ATZ-p1.moments')
           excludepix=[-100,0.0006],outfile='M51-ATZ-p1.moments')
</source>
</source>
Line 311: Line 322:
=== 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). [[simdata2|Simdata2]] 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 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.1|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.


[[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 3.1}}
{{Checked 3.0.2}}


--[[User:Jgallimo|Jack Gallimore]] 15:55, 30 April 2010 (UTC)
--[[User:Jgallimo|Jack Gallimore]] 15:55, 30 April 2010 (UTC), updated for CASA 3.1 M. Lacy 1st December 2010

Latest revision as of 19:00, 28 April 2011


Simulating Observations in CASA 3.1

Remember that you need to follow this process: Extracting scripts from these tutorials to extract a casapy script from this web page.

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. 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.

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

Load these data into CASA. For convenience, store the name of the resulting measurement set into the python global cubeName.

cubeName = 'm51-song'
importfits(fitsimage='NGC5194.bima12m.cm.fits', imagename=cubeName)
# initialize simdata
default("simdata")
skymodel=cubeName
modifymodel=True

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.

#z's 
# Distinguish between z_lsrk, which sets the observed frequency for scaling, from 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

The convention is 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 simdata. First, we'll set up some globals to establish some file naming conventions.

suffix = "-p1" # z = 0.1, or point-1; useful to distinguish from repeated simulations at different z's
project = "M51-ATZ" + suffix  # project ID to assign output filenames (simdata parameter)


Flux Density Scaling

Our goal in this step is to set the inbright parameter of simdata, which will scale the data cube appropriately for its new luminosity distance.

Simdata wants models in units of Jy / pixel, but the BIMA SONG cube is in units of Jy / beam. That's an easy conversion.

# BIMA SONG beam 
bmaj = imhead(imagename=cubeName,mode='get',hdkey='beammajor')
bmin = imhead(imagename=cubeName,mode='get',hdkey='beamminor')
#
# use qa tool to convert beam in radians to (1 arcsec) pixels:
bmaj = qa.convert(bmaj,'arcsec')['value']
bmin = qa.convert(bmin,'arcsec')['value']
#
toJyPerPix = 1.0 / (1.1331 * bmaj * bmin) # gaussian beam conversion = beams / pixel

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. Simdata can scale the peak surface brightness using the parameter inbright.

# correct flux density for luminosity distance
fluxScale = (dl_old/dl_new)**2 * (1.0 + z_new) / (1.0 + z_old_cmb)
# current peak flux:
peak = imstat(cubeName)['max'][0] * toJyPerPix # need to convert to Jy / pixel
# desired peak flux in Jy / pixel:
inbright = "%fJy/pixel" % (peak*fluxScale) # use python formatting convention
# inbright is a simdata parameter

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 simdata parameters indirection (change the location of M51 on the sky) and incell (change the angular scale of an input pixel). a To perform task (1), we'll just flip the sign of the declination using imhead. Notice the use of the qa tool to convert radians to sexagesimal.

# Move to southern hemisphere.
# Use qa.formxxx tool to convert rad to sexagesimal.
# For clarity, build up "indirection" string one term at a time.
# Epoch
indirection = "J2000 " # parameter for simdata
# RA
indirection += qa.formxxx(imhead(imagename=cubeName,mode='get',hdkey='crval1')['value']+'rad',format='hms') + " " 
# Dec * -1
indirection += qa.formxxx('%frad' % (-1*float(imhead(imagename=cubeName,mode='get',hdkey='crval2')['value'])),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, and mode="put" to store the new pixel scale in the model header.

# scale pixel size:  imhead returns things in rad, so convert to arcsec
oldCell = float(imhead(imagename=cubeName,mode='get',hdkey='cdelt2')['value']) * 206265 
# scale for new angular size distance
newCell = oldCell * da_old / da_new 
# Format the new pixel size for input to simdata
incell = "%farcsec" % (newCell) # parameter for simdata

Adjusting the Frequency Axis

Changing the frequency axis of the model header is just a straightforward (1+z) correction. Simdata 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 Simdata 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.

# move freq to z_new
oldFreq = float(imhead(imagename=cubeName,mode='get',hdkey='crval3')['value'])   # Hz
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 = float(imhead(imagename=cubeName,mode='get',hdkey='cdelt3')['value'])   # Hz
newDnu = abs((1.0+z_old_lsrk) /(1.0+z_new)*oldDnu) # need positive channel widths to ensure noise calc goes well
inwidth = "%fHz" % newDnu # parameter for simdata
#
# Specify the observing frequency at the center of the observing band:
incenter = "%fHz" % (newFreq + 0.5*nchan*newDnu)

Simdata

The CASA task simdata will monolithically simulate an ALMA observation, produce measurement sets with and without thermal noise, and finally produce a CLEANed image cube based on the simulated observation. Details can be found in the Simulating Observations in CASA tutorial.

The original BIMA SONG image is about 480 arcseconds across; scale this image size to the new redshift.

# estimate final image size
imSize = 480.0 * (da_old / da_new) # in arcseconds

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. Simdata needs the spacing between pointings in the mosaic; we'll require pointings spaced by half of the primary beam.

# mosaicking info
primaryBeam = 17.0 * (300e9 / newFreq) # in arcseconds; ALMA primary beam = 17 arcsec at 300 GHz
pointingspacing = "%farcsec" % (primaryBeam / 2.0) 
mapsize = "%farcsec" % imSize # how big to make the mosaic

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.

# Estimate desired beam size. BIMA SONG has 5": use 15" projected to new redshift
beamNew = 15.0 * (da_old / da_new)

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.

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.

beamNew = 4.0 * pixelSize

Now we know both the image size and pixel size in arcseconds, but simdata 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.

imSizePix = int(round(imSize / pixelSize))
if imSizePix < 256: imSizePix = 256

Let simdata decide on an appropriate ALMA configuration based on the desired beam size. Set the parameter antennalist as follows (but see Other Antenna Configurations below).

predict=True
antennalist = "alma;%farcsec" % beamNew

Now we have enough information to run simdata, and hopefully some of the python global variables that were defined above will start to make sense. The following CASA and python commands set up the remaining parameters for the simdata task.

modelimage = cubeName
integration = '10s' # watch out for memory limits vs. ability to complete mosaic here
# 10s is usually safe for large mosaics, but will require more memory
totaltime = '28800s' # 8 hr integration

If you are doing this as a tutorial, we suggest that you consider increasing the integration time so that your learning is not slowed by processing time See Etime study for details.

# make simulated images/cubes
image=True
thermalnoise="tsys-atm" # add thermal noise, produce noisy.ms measurement set
vis = '$project.noisy.ms' # clean the data with *thermal noise added* 
cell = "%farcsec" % pixelSize
imsize = [imSizePix,imSizePix]
threshold = "1.0mJy"
weighting = "natural"
stokes = 'I'
#
verbose = True
graphics="both"
overwrite = True
#
inp("simdata")
simdata()


Other Antenna Configurations

ALMA synthetic beam size as a function of array configuration number

Finally, we need to know which ALMA configuration number based on the desired angular resolution. 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.

Notice from the 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.0.1 and 3.1, here is a list of included antenna configurations.

ALMA Configuration Files
alma.csv.late.cfg
alma.csv.mid.cfg
alma.early.large.cfg
alma.early.med.cfg
alma.out01.cfg
alma.out02.cfg
alma.out03.cfg
...
alma.out27.cfg
alma.out28.cfg

There are also configuration data for the ACA, EVLA, and SMA.

Take a Break

If you have got this far, you've earned it. Simdata will be running for a while, and coffee sure sounds good right now.

Results

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.

Here is an inventory of some of the simdata products.

Filename Description
M51-ATZ-p1.ms Model measurement set sans noise
M51-ATZ-p1.noisy.ms Model measurement set with thermal noise
M51-ATZ-p1.clean.image CLEAN-deconvolved cube of M51-ATZ-p1.noisy.ms

And there are plenty of other auxiliary files.

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

imhead(imagename="M51-ATZ-p1.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.

immoments(imagename='M51-ATZ-p1.image',moments=[0,1],axis='spectral',
          excludepix=[-100,0.0006],outfile='M51-ATZ-p1.moments')

The results are shown at right.

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 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 3.1

Last checked on CASA Version 3.0.2 (r11631).

--Jack Gallimore 15:55, 30 April 2010 (UTC), updated for CASA 3.1 M. Lacy 1st December 2010