M51 at z = 0.1

This article is under construction. Watch this space!

Overview

This tutorial presents a simulation of ALMA observations of a well-known galaxy, M51, as it would be observed at redshift z = 0.1.

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 antennas, and thermal noise levels appropriate for the ALMA site. Calibration overheads are not included, nor is phase noise owing to a varying troposphere. As such, this simulation should be viewed as somewhat optimistic.

For this tutorial, we'll use the BIMA SONG 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

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

cubeName = 'm51-song'
importfits(fitsimage='NGC5194.bima12m.cm.fits', imagename=cubeName)

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 with the default cosmology; redshifts are obtained from NED.

#z's 
z_old_cmb = 0.002122 # CMB-referenced z for cosmological distances
z_old_lsrk = 0.001544 # z_obs from BIMA header
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
cubeOut = 'm51-atz' + suffix + '.im' # name for the model image (input to simdata)

Flux Density Scaling

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 -- could use imhead mode = 'get' to automate this step
bmaj = 5.82
bmin = 5.08 # note: BIMA SONG pixels = 1 arcsec, so this is the fwhm in pixels, too
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.

# correct flux density for luminosity distance
fluxScale = (dl_old/dl_new)**2 * (1.0 + z_new) / (1.0 + z_old_cmb)

Notice that there is an additional (1+z) correction because we are scaling flux densities rather than bolometric fluxes. Now we can use immath to perform the scaling. Immath needs an string expression to describe the scaling.

fluxExpression = "(IM0 * %f * %f)" % (toJyPerPix, fluxScale)

Here, fluxExpression is just a python global that will feed into immath, IM0 symbolically represents the input image (the BIMA SONG cube), and we use python formatting conventions to transfer the Jy/pix and flux scaling conversions into the expression string.

Now, generate the model cubeOut = m51-atz-p1.im.

immath(imagename=cubeName, outfile=cubeOut, mode='evalexpr',
       expr=fluxExpression) # now in Jy / pixel at new redshift

Just for bookkeeping's sake, we'll change the brightness unit in the header of the model image.

# adjust brightness units
hdval = 'Jy/pixel'
dummy = imhead(imagename=cubeOut, mode = 'put', hdkey='BUNIT', 
               hdvalue=hdval)

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. To accomplish (1), we'll just flip the sign of the declination using imhead.

# move to southern hemisphere
hdval = imhead(imagename=cubeName, mode = 'get', hdkey='CRVAL2')
hdval['value'] = -1.0 * float(hdval['value'])
dummy = imhead(imagename=cubeOut, mode = 'put', hdkey='CRVAL2', 
               hdvalue=hdval)

The global variable "dummy" is just a throwaway to store the status of the imhead operation.

Next, adjust the pixel scale for the new angular size distance.

# x-pixel
hdval = imhead(imagename=cubeName, mode = 'get', hdkey='CDELT1')
# correct for new angular size distance
hdval['value'] = (da_old / da_new)*float(hdval['value']) 
dummy = imhead(imagename=cubeOut, mode = 'put', hdkey='CDELT1', 
               hdvalue=hdval)

# y-pixel
hdval = imhead(imagename=cubeName, mode = 'get', hdkey='CDELT2')
# correct for new angular size distance
hdval['value'] = (da_old / da_new)*float(hdval['value']) 
dummy = imhead(imagename=cubeOut, mode = 'put', hdkey='CDELT2', 
               hdvalue=hdval)