View source for MRK 6: red-shifted HI absorption 6.5.2

{{Checked 6.5.2}}


Disclaimer: Due to continuous CASA software updates, GUI images may look different on more recent versions than those shown here.


[[Category: VLA]][[Category: Calibration]] [[Category: Imaging]] [[Category: Spectral Line]]

== Overview ==

This example demonstrates a trickier spectral line data set. In this case, several radio bright active galaxies were observed to look for redshifted 21 cm absorption. The data were obtained using 4IF mode, which means that the full range of velocities were split into two spectral windows (in this example, CASA identifies them as spectral windows 10 and 11). Producing the data cube will require eventually stitching these windows back together. Happily, {{tclean}} includes flexible regridding and interpolation in the spectral domain. {{tclean}} further allows you to adjust the velocity reference frame.

== Retrieve the data from the NRAO archive ==

Download your data from the [http://data.nrao.edu: VLA Archive]; in this example we'll use the publicly available survey AB658 observed in 1992. These data were published in [http://adsabs.harvard.edu/abs/1999ApJ...524..684G Gallimore et al. (1999)]. 

With the present archive defaults, you should have obtained the following VLA archive file:

<pre>
'AB658_1_48948.07214_48948.73140.exp'
</pre>

This tutorial broadly follows the techniques described in the [[Calibrating a VLA 5 GHz continuum survey|continuum survey tutorial]] and the [[NGC 5921: red-shifted HI emission|21 cm emission tutorial]], and the basic calibration steps are presented only in outline.

== Loading the Data ==

<source lang="python">
# In CASA
# Import the archive file into a CASA Measurement Set (MS).  
importvla(archivefiles='AB658_1_48948.07214_48948.73140.exp', vis='ab658.ms')

</source>

<!--
<source lang="python">

measurementSet = "ab658.ms" # Store the desired name of the measurement set
#
from glob import glob
fileList = sorted(glob("AB658*.xp?"))
#
# double check to see if the file already exists
import os
#
if not os.path.isdir("./" + measurementSet): # measurement sets are directories
    importvla(archivefiles=fileList, vis=measurementSet) # only import if ab658.ms doesn't yet exist

vis = measurementSet
mode = "summary"
listobs
</source>
-->

We will now generate a text file listing of the details for the observations using the {{listobs}}.

<source lang="python">
# In CASA
# generate a listobs file with information about the observations as well as printing some details to the terminal
listobs(vis='ab658.ms', listfile='ab658.listobs')
</source>

Next we will generate a plot of antenna positions. This will help identify a good reference antenna.

<source lang="python">
# In CASA
# plot antenna positions
plotants(vis='ab658.ms', figfile='ab658_plotants.png')

</source>

For the purposes of this tutorial, we are interested only in the source [http://en.wikipedia.org/wiki/B._E._Markarian Markarian] [http://nedwww.ipac.caltech.edu/cgi-bin/nph-objsearch?objname=mrk+6&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 6], but the reduction techniques could be applied to any of the sources in the measurement set. 

[[File:MKN6_Listobs.png|thumb|{{listobs|Listobs}} output for AB658. The data for MKN 6 and its calibrators are highlighted.]]

Using {{listobs}} and {{plotants}} we learn the source names, spectral windows, and viable reference antennas. A screenshot of the {{listobs}} output is given at right.

{| border=1
| '''Source''' || MKN6
|-
| '''Flux Density Scale Cal''' || 0134+329 = 3C48
|-
| '''Complex Gain Cal''' || 1003+830
|- 
| '''Central Antennas''' || VA27, VA09, VA23
|-
| '''Spectral Window IDs''' || 10, 11
|}

<!--
It's just as well to put that information into python global variables.

<source lang="python">
sourceName = "MKN6"
CG_Cal = "1003+830"
FDS_Cal = "0134+329"
refAnt = "VA27" # or perhaps "VA09" or "VA23"
spwStr = '10, 11'
</source>

The selection of the spectral windows, stored here in the string variable ''spwStr'', is key to getting the calibration right. Each science target in the multisource measurement set has a tuning respective of the source redshift, and the calibration observations are matched to those spectral windows. Were the spectral windows parameter (usually '''spw''') left to its default value (all available spectral windows), many calibration routines would necessarily fail.
-->

The selection of the spectral windows is key to getting the calibration right. Each science target in the multisource measurement set has a tuning respective of the source redshift, and the calibration observations are matched to those spectral windows. Were the spectral windows parameter (usually '''spw''') left to its default value (all available spectral windows), many calibration routines would necessarily fail.

== Editing the Data ==

[[File:MKN6_Amp_vs_Time.png|thumbnail|Plotms display of (uncalibrated) visibility amplitudes vs. time for the source MKN 6.]]

First inspect the data using {{plotms}}. 

<source lang="python">
# In CASA
plotms(vis='ab658.ms', field='MKN6', xaxis='time', yaxis='amp')
</source>

<!--
<source lang="python">
plotms(vis = measurementSet, field = sourceName)
</source>


Be sure to select '''Axes &rarr; X Axis &rarr; Time''' and '''Axes &rarr; Y Axis &rarr; Amp'''. 
-->

There is some obvious bad data isolated in time during the observation; see the screenshot at right. We can use the locate tool in plotms to identify specific time ranges or antennas that are causing the problems and flag the bad data with the CASA task flagdata. For more information about flagging data in general, please see our [http://casaguides.nrao.edu/index.php/VLA_CASA_Flagging Topical Guide: Flagging VLA Data]. Hopefully you find three time ranges that are bad, which will require three calls to flagdata to flag the bad baselines.


<source lang="python">
# in CASA
flagdata(vis='ab658.ms',antenna='9&12', spw='10', timerange='1992/11/22/09:52:30.00')
flagdata(vis='ab658.ms',antenna='8&18', spw='10', timerange='1992/11/22/09:58:10.00')
flagdata(vis='ab658.ms',antenna='8&9', spw='10', timerange='1992/11/22/10:28:30.00')

</source>

Replot the data to confirm flagging has done what is expected.

== Calibration ==

Calibration of HI absorption spectra follows the basic techniques outlined in the [[Calibrating a VLA 5 GHz continuum survey|continuum tutorial]] and [[NGC 5921: red-shifted HI emission|HI emission tutorial]].

First, set the flux denisty of the flux density scale calibrator. Notice that we are being careful to keep track of the spectral windows and not using the default ''spw''.

<source lang="python">
# In CASA
# set the flux denisty scale standard
setjy(vis='ab658.ms', field='0134+329', spw='10,11', model='', usescratch=True)
</source>

Bandpass calibration comes next. 

<source lang="python">
# In CASA
# bandpass calibration
bandpass(vis='ab658.ms', spw='10,11', caltable='ab658.bcal', gaintable='', gainfield='', interp='', field='0134+329', selectdata=False, bandtype='B', solint='inf', combine='scan', refant='VA27')
</source>

<!--
<source lang="python">
# set the flux density scale standard
default('setjy')
vis = measurementSet
field = FDS_Cal
spw = spwStr # this is key to getting the calibration right!
modimage = ''
setjy()
</source>

Bandpass calibration comes next. 

<source lang="python">
# set some default table names
btable = 'ab658.bcal' # bandpass
gtable = 'ab658.gcal' # gain
ftable = sourceName + '.fluxscale' # flux density scale
#
if os.path.exists(btable): rmtables(btable)
#
if os.path.exists(gtable): rmtables(gtable)
#
if os.path.exists(ftable): rmtables(ftable)
#
# bandpass calibration
default('bandpass')
vis = measurementSet
spw = spwStr # this is key to getting the calibration right!
caltable = btable
gaintable = ''
gainfield = ''
interp = ''
field = FDS_Cal
selectdata = False
opacity = 0.0
bandtype = 'B'
solint = 'inf'
combine = 'scan'
refant = refAnt
bandpass()
</source>
-->

[[File:MKN6_PlotCal_Bandpass_Cal_Curves2.png|thumb|Bandpass calibration curves. Note that the response curve is flat between channels 2 and 27.]]

In order to perform the flux density scale and complex gain calibrations, we'll need to inspect the bandpass response curve to select appropriate channels to average. We want to choose those channels where the response is relatively flat. The result of plotms is shown at right: the bandpass is flat over channels 2:27. 

<source lang="python">
# In CASA
# plot the bandpass calibration, amp vs channel
plotms(vis='ab658.bcal', spw='10~11', field='0134+329', gridrows=2, gridcols=1, rowindex=0, plotindex=0, yaxis='gainamp', xaxis='channel', coloraxis='Antenna1')

# plot the bandpass calibration, gain phase vs channel
plotms(vis='ab658.bcal', spw='10~11', field='0134+329', gridrows=2, gridcols=1, rowindex=1, plotindex=1, yaxis='gainphase', xaxis='channel', coloraxis='Antenna1', clearplots=False)
--
</source>

<!--
<source lang="python">
# plot the bandpass calibration
default('plotcal')
spw = spwStr
caltable = btable
field = FDS_Cal
subplot = 211
yaxis = 'amp'
showgui = True
plotcal()
subplot = 212
yaxis = 'phase'
plotcal()
</source>
-->

The following steps generate amplitude and complex gain solutions for the calibrators, and the fluxscale is bootstrapped from the flux density scale calibrator onto the complex gain calibrator, effectively turning the complex gain calibrator into a local flux density scale calibrator for the science target.

<source lang="python">
# In CASA
# generate amplitude and gain solutions for complex gain and flux density scale calibrators
gaincal(vis='ab658.ms', caltable='ab658.gcal', gaintable='ab658.bcal', gainfield='', interp='nearest', field='0134+329, 1003+830', 
        spw='10:2~27, 11:2~27', selectdata=False, gaintype='G', solint='inf', combine='', calmode='ap', minsnr=1.0, refant='VA27')
</source>

Note: CASA will print this message: "The the following spws have no corresponding cal spws in ab658.bcal: 0 1 2 3 4 5 6 7 8 9 12 13 14 15 16 17 18 19 20 21 22 23". This is expected because spw='10,11' are the only two spws being calibrated, spw='10,11' are not in the list presented by CASA.

<source lang="python">
# In CASA
# Bootstrap flux density scale
fluxscale(vis='ab658.ms', fluxtable='MKN6.fluxscale', caltable='ab658.gcal', reference='0134+329', transfer='1003+830')
</source>

<!--
<source lang="python">
default('gaincal')
vis = measurementSet
caltable = gtable
gaintable = btable
gainfield = ''
interp = 'nearest'
field = FDS_Cal + ',' + CG_Cal
# specify data only from the flat region of the bandpass response
spw = '10:2~27,11:2~27'
selectdata = False
opacity = 0.0
gaintype = 'G'
solint = 'inf'
combine = ''
calmode = 'ap'
minsnr = 1.0
refant = refAnt
gaincal()
#=====================================================================
#
# Bootstrap flux density scale
#
vis = measurementSet
# set the name for the output rescaled caltable
fluxtable = ftable
caltable = gtable
reference = FDS_Cal
transfer = CG_Cal
fluxscale()
</source>
-->

The task {{fluxscale}} will give the bootstrapped flux density of the complex gain calibrator. If all goes well, something like the following should appear in the log window.

<pre>
Beginning fluxscale--(MSSelection version)-------
 Found reference field(s): 0134+329
 Found transfer field(s):  1003+830
 ...
 Flux density for 1003+830 in SpW=10 (freq=1.3921e+09 Hz) is: 0.492452 +/- 0.00331148 (SNR = 148.711, N = 52)
 Flux density for 1003+830 in SpW=11 (freq=1.39425e+09 Hz) is: 0.489206 +/- 0.00287348 (SNR = 170.249, N = 52)
 ...
 Storing result in ab658.fluxscale
</pre>

We should inspect the calibration solutions to make sure they are reasonable, but we'll be [http://en.wikipedia.org/wiki/Virginia_Cavaliers cavalier] and tread forward. Apply the calibrations to the science target and split that target into its own measurement set.

<source lang="python">
# In CASA
# apply calibrations to science target
applycal(vis='ab658.ms', gaintable=['MKN6.fluxscale','ab658.bcal'], gainfield=['1003+830', '*'], interp=['linear','nearest'], field='MKN6', spw='10,11', spwmap=[], selectdata=False)
</source>

<source lang="python">
#In CASA
# split out the corrected data column (which is placed in the data column in the new MS) for our science target.
split(vis='ab658.ms', outputvis='ab658-MKN6.split.ms', field='MKN6', spw='10:2~27, 11:2~27', datacolumn='corrected')
</source>

* ''spw='10:2~27, 11:2~27'' : Note that we have selected only the inner channels, 2 through 27 for spectral windows 10, 11. 

<!--
<source lang="python">
default('applycal')
vis = measurementSet
gaintable = [ftable,btable]
gainfield = [CG_Cal,'*']
interp = ['linear','nearest']
field = sourceName
# choose the appropriate spectral windows
spw = spwStr
# by default, calibrations for a spectral window will map to the same
spwmap = []
selectdata = False
opacity = 0.0
applycal()
#
default('split')
splitms = 'ab658-MKN6.split.ms'
if os.path.isdir(splitms): rmtables(splitms)
#
vis = measurementSet
outputvis = splitms
field = sourceName
# weed out the end channels
spw = '10:2~27,11:2~27'
datacolumn = 'corrected'
split()
</source>
-->

== Continuum Subtraction ==

Ideally, we would use {{plotms}} to identify the channels where HI absorption is present (see [[NGC 5921: red-shifted HI emission#Plot the Spectrum|the HI emission tutorial.]]) The absorption line is however weak and isolated against one of a number of radio continuum sources in the data; as such, the absorption is not apparent in averaged visibilities. Instead, we'll have to produce a dirty cube using {{tclean}} and inspect the spectrum in the image domain.

=== Search for the Absorption Line ===

{{uvcontsub|Uvcontsub}} currently identifies continuum ranges by channel rather than frequency or velocity. We'll make two cubes, one for each spectral window. Note that the act of using split re-indexes the spectral windows and channels. Therefore channels 2~27 have been reindexed and are now channels 0~25 for both channels 10 & 11 (and spw='10,11' have now reindexed as spw='0,1').

<source lang="python">
# In CASA
# image cube for spectral window 0
tclean(vis='ab658-MKN6.split.ms', imagename='temp0', spw='0:0~25', specmode='cube', interpolation='nearest', niter=2000, 
      cell=['0.4arcsec', '0.4arcsec'], phasecenter='J2000 06h52m12.2 +74d25m37')
</source>

<source lang="python">
# In CASA
# make a new image cube with spectral window 1
tclean(vis='ab658-MKN6.split.ms', imagename='temp1', spw='1:0~25', specmode='cube', interpolation='nearest', niter=2000, 
      cell=['0.4arcsec', '0.4arcsec'], phasecenter='J2000 06h52m12.2 +74d25m37')
</source>

<!--
<source lang="python">
default(clean)
vis = splitms
imagename = 'temp0'
spw = '0:0~31'
mode = 'channel'
nchan = -1
start = 0
width = 1
interpolation = 'nearest'
niter = 2000
cell = ['0.4arcsec', '0.4arcsec']
phasecenter = 'J2000 06h52m12.2 +74d25m37'
clean()
#
imagename='temp1'
spw='1:0~31'
clean()
</source>
-->

It's not ordinarily necessary to give the phase center, but this is a good opportunity to use the regridding capabilities of {{tclean}}. The data were originally taken in B1950 coordinates, and the observations followed the old practice of offsetting the pointing center slightly away from the source coordinates (this is not done with JVLA observations today). Here, we specify the J2000 coordinates (looked up on [http://nedwww.ipac.caltech.edu NED]) and allow {{tclean}} to handle the precession and centering that we require.

Now we can inspect the cubes which have been named '''temp0.image''' and '''temp1.image'''. CASA will depreciate old {{viewer}} task, and no development for this task has been done for some time. So, we will instead use CARTA. 
<!--
<source lang="python">

<source lang="python">
# In CASA
viewer(infile='temp0.image', displaytype='raster')
</source>

</source>
-->

[[File:MKN6_SpectralProfile.png|thumb|Spectral profile for MKN 6 in the 2nd spectral window, '''spw = 1'''. The weak absorption line can be seen around channel 15.]]

Now, use the Spectral Profiler tool in CARTA to generate a spectrum. The source spectrum in '''spw = 1''' is shown in the figure at right.

* After opening the image, use the ellipse tool to make a region around the source. 
* Next, open the Spectral Profiler tool. Since the default will be to show the x axis in LSRK, click the gear icon to open the settings of the Spectral Profiler: 
* Change the Coordinate to Channel in the Conversion tab and Line Style to Lines in the Styling tab. 
* ''Optional:'' Save the numerical data to a file. With the mouse over the SPectral Profiler plot, there should be two icons at the lower right: one looks like a disk for saving the image, the other a table for saving the numerical data: click the table button to save a csv file. 


<!--
<source lang="python">
infile = 'temp1.image'
displaytype = 'raster'
viewer()
</source>
-->

=== Perform the Continuum Subtraction ===

Looking over the plots and text files, only channels 0~25 in either spectral window contain usable data. The absorption line is located in 1:15 and surrounding channels. Just to play it safe, we'll use the channels outside of 1:10~20 (spectral window 1, channels 10-20) to define the continuum for subtraction.

<source lang="python">
# In CASA
uvcontsub(vis='ab658-MKN6.split.ms', fitspw='0:0~25, 1:0~9; 21~25', solint='int', fitorder=1)
</source>

<!--
<source lang="python">
default(uvcontsub)
vis = 'ab658-MKN6.split.ms'
# remove old versions of the continuum subtracted data
if os.path.exists(vis + '.contsub'): rmtables(vis+'.contsub')
#
fitspw = '0:0~25,1:0~9;21~25' # define the continuum windows for each spectral window
solint = 'int'
fitorder = 1 # fit linear baselines
uvcontsub()
</source>
-->

{{uvcontsub|Uvcontsub}} produces the uv subtracted measurement set:
<!-- * ab658-MKN6.split.ms.cont (continuum based on fit to linefree channels), & -->
* ab658-MKN6.split.ms.contsub (continuum-subtracted visibilities).

== Imaging the Continuum Subtracted Cube ==

[[File:MKN6_Viewer_and_HI_Line.png|thumb|Viewer windows displaying the HI (21cm) absorption line cube of MKN 6. The viewer window shows the channel containing the strongest absorption (at 5479 km/s), and the spectral profile window shows the average spectrum around the position of the strongest absorption.]]

Inspecting the preliminary (continuum-containing) cubes, temp0.image and temp1.image, that we generated [[#Search for the Absorption Line|above]], we can determine good parameters for the velocity axis.

<source lang="python">
# In CASA
tclean(vis='ab658-MKN6.split.ms.contsub', spw='', imagename='MKN6-HIABS', specmode='cube', start='-252.349km/s', interpolation='linear',          
      outframe='BARY', veltype='optical', restfreq='1420.405752MHz', niter=2000, threshold='1.0mJy', weighting='natural', cell=['0.4arcsec', '0.4arcsec'], 
      phasecenter = 'J2000 06h52m12.2 +74d25m37')
</source>

* threshold='1.0mJy' : avoid over-cleaning since it appears that the channel rms is about 0.5 mJy/beam
* cell = ['0.4arcsec', '0.4arcsec'] : VLA-A configuration expects 1.5--2 arcsec beams. 0.4arcsec pixels ensure sampling the beam by better than 3x

<!--
<source lang="python">
default(clean)
vis = 'ab658-MKN6.split.ms.contsub'
spw=''
imagename = 'MKN6-HIABS'
mode = 'velocity'
nchan = 48
start = '5380km/s'
# This will set the channel width equal to the first channel in the first spectral window
width = ''
interpolation = 'linear'
outframe = 'bary'
veltype = 'optical'
restfreq = "1420.405752MHz"
</source>

It also appears that the channel rms is about 0.5 mJy/beam; to avoid over-cleaning, we'll set the threshold at 1 mJy/beam.

<source lang="python">
niter = 2000
threshold = '1.0mJy'
weighting = 'natural'
</source>


Set the source coordinates. We'll use defaults for the rest of the parameters; with the continuum subtracted, the data have only simple structure with no confusing sources. 

<source lang="python">
# VLA-A; expect 1.5--2 arcsec beams. 0.4arcsec pixels ensure sampling the beam by better than 3x
cell = ['0.4arcsec', '0.4arcsec']
phasecenter = 'J2000 06h52m12.2 +74d25m37'
clean()
</source>
-->

The ({{viewer}}) results are illustrated in the figure at right.

=== Imaging the Continuum Dataset ===

[[File:MKN6_Continuum.png|thumb|21cm continuum image of MKN 6. The image has been scaled, using (Data)(Adjust)(Basic Settings) '''Scaling Power Cycles''' = -1.9, to bring out artifacts apparent in the background. The 6-pointed star pattern around the source is symptomatic of residual phase-errors, which can be healed using [[Calibrating a VLA 5 GHz continuum survey#Self-calibration|self-calibration]], and the arcuate stripes across the field are sidelobes of neighboring sources that need to be cleaned using [[Imaging Flanking Fields|flanking fields]].]]

We can similarly image the continuum data with just a few changes of the {{tclean}} parameters. Here, we call tclean by selecting the spectral windows and channels pertaining to only the continuum.

<source lang="python">
# In CASA
tclean(vis='ab658-MKN6.split.ms', specmode='mfs', spw='0:0~25,1:0~9;21~25', imagename='MKN6-CONT', restfreq='1420.405752MHz', niter=2000, 
      threshold='1.0mJy', weighting='natural', cell=['0.4arcsec', '0.4arcsec'], phasecenter = 'J2000 06h52m12.2 +74d25m37')
</source>

<!--
<source lang="python">
tget clean
vis = 'ab658-MKN6.split.ms'
mode = 'mfs'
spw = '0:0~25,1:0~9;21~25'
imagename = 'MKN6-CONT'
clean()
</source>
-->

The result is shown at right. The continuum image needs some more help. Phase errors produce the 6-pointed star pattern surrounding the source and can be healed using [[Calibrating a VLA 5 GHz continuum survey#Self-calibration|self-calibration]]. The striping across the field is actually the combination of sidelobes from neighboring sources; these sidelobes can be suppressed by cleaning [[Imaging Flanking Fields|flanking fields]].

== Extract the Absorption Spectrum ==

The {{viewer|Viewer}} provides a handy first-look at the spectrum, but the extracted spectrum is an average over a spatial box. What we really want is the ''integrated'' spectrum, not the average. To do this, we'll use CASA's {{imstat}}, {{imval}}, and {{imhead}} tools to extract cube information to python variables and then let python build up the spectrum.

=== Extract the Spectrum: load the cube data into python arrays ===

First, we need to determine where in the cube to extract the spectrum; equivalently, we need to find the location of the minimum value in the cube. {{imstat|Imstat}} is the tool for this.

<source lang="python">
# In CASA
cubeStat=imstat("MKN6-HIABS.image")
</source>

The variable <tt>cubeStat</tt> is now a [http://docs.python.org/tutorial/datastructures.html#dictionaries python dictionary]. To see what it contains, just enter the variable name:

<source lang="python">
# In CASA
cubeStat
</source>

and the following summary of <tt>cubeStat</tt> appears.

<pre>
{'blc': array([0, 0, 0, 0], dtype=int32),
 'blcf': '06:52:17.165, +74.25.16.997, I, 1.39542e+09Hz',
 'max': array([ 0.00286199]),
 'maxpos': array([27, 68,  0,  6], dtype=int32),
 'maxposf': '06:52:14.485, +74.25.44.199, I, 1.39484e+09Hz',
 'mean': array([ -5.30860147e-06]),
 'medabsdevmed': array([ 0.00031676]),
 'median': array([ -3.58635123e-06]),
 'min': array([-0.01634453]),
 'minpos': array([49, 53,  0, 10], dtype=int32),
 'minposf': '06:52:12.299, +74.25.38.200, I, 1.39444e+09Hz',
 'npts': array([ 470000.]),
 'q1': array([-0.00032042]),
 'q3': array([ 0.00031312]),
 'quartile': array([ 0.00063354]),
 'rms': array([ 0.00050688]),
 'sigma': array([ 0.00050685]),
 'sum': array([-2.49504269]),
 'sumsq': array([ 0.12075419]),
 'trc': array([99, 99,  0, 46], dtype=int32),
 'trcf': '06:52:07.331, +74.25.56.597, I, 1.39093e+09Hz'}
</pre>

These entries can be accessed individually like an array variable. For example, the following commands copy the ''x'' and ''y'' positions of the absorption line to scalar variables. 

<source lang=python>
# In CASA
%cpaste

# Press Enter or Return, then copy/paste the following:
x2extract = cubeStat['minpos'][0]
y2extract = cubeStat['minpos'][1]
--
</source>

Now, extract the spectrum using {{imval}}.

<source lang=python>
# In CASA
%cpaste

# Press Enter or Return, then copy/paste the following:
extractBox = "%d,%d" % (x2extract, y2extract) # copy the position to a string
cubeSpec = imval("MKN6-HIABS.image", box=extractBox, stokes='I')
--
</source>

<!--
<source lang=python>
x2extract = cubeStat['minpos'][0]
y2extract = cubeStat['minpos'][1]
</source>

Now, extract the spectrum using {{imval}}.

<source lang=python>
extractBox = "%d,%d" % (x2extract, y2extract) # copy the position to a string
cubeSpec = imval("MKN6-HIABS.image", box=extractBox, stokes='I')
</source>
-->

<tt>cubeSpec</tt> is another [http://docs.python.org/tutorial/datastructures.html#dictionaries python dictionary], and the spectrum is stored in <tt>cubeSpec['data']</tt>.

<div style="background-color: #dddddd;">
'''Comment:''' The absorption line source is unresolved by this observation, or, to put it another way, is entirely contained within a synthesized beam. The extraction at a single spatial pixel gives the integrated spectrum and therefore changes units from surface brightness (mJy/beam) to flux density (mJy). However, the data have been forced into a rectangular grid, and there's no guarantee that the absorption minimum centers within a pixel. Better to use {{imfit}} to determine more accurately the position of the minimum, and perhaps re-{{tclean}} with a new phasecenter. To avoid complicating the present tutorial, we'll just stick with the minimum in the current {{tclean}} grid.
</div>

=== Generate the frequency axis ===

We can derive the frequency axis from the header.

<source lang="python">
# In CASA
%cpaste

# Press Enter or Return, then copy/paste the following:
import numpy # make sure vector and array arithmetic options are loaded
cubeHead = imhead("MKN6-HIABS.image", mode="list")
nSpec = cubeStat['trc'][3] + 1 # get the number of frequency channels. 
f0 = float(cubeHead['crval4']) # reference freq in Hz
df = float(cubeHead['cdelt4']) # channel width in Hz
i0 = cubeHead['crpix4'] # reference pixel
freqSpec = (numpy.arange(nSpec) - i0)*df + f0
--
</source>

<!--
<source lang="python">
import numpy # make sure vector and array arithmetic options are loaded
cubeHead = imhead("MKN6-HIABS.image", mode="list")
nSpec = cubeStat['trc'][3] + 1 # get the number of frequency channels. 
f0 = float(cubeHead['crval4']) # reference freq in Hz
df = float(cubeHead['cdelt4']) # channel width in Hz
i0 = cubeHead['crpix4'] # reference pixel
freqSpec = (numpy.arange(nSpec) - i0)*df + f0
</source>
-->

=== Plot the spectrum using matplotlib ===

[[File:MKN6_Absorption_Spectrum.png|thumb|Absorption spectrum of MKN 6 in frequency space.]]

Now we have the spectrum stored in <tt>freqSpec</tt>, <tt>cubeSpec['data']</tt>. Use [http://matplotlib.sourceforge.net/api/pyplot_api.html matplotlib.pyplot] to plot the spectrum (see figure at right).

<source lang="python">
# In CASA
%cpaste


# Press Enter or Return, then copy/paste the following:
import matplotlib.pyplot as plt
plt.clf() # clear the plot (figure)
plt.plot(freqSpec/1.e9, cubeSpec['data'], 'k-')
plt.xlabel("Freq (GHz)")
plt.ylabel("Flux Density (mJy)")
plt.show()
--
</source>

<!--
<source lang="python">
import matplotlib.pyplot as plt
plt.clf() # clear the plot (figure)
plt.plot(freqSpec/1.e9, cubeSpec['data'], 'k-')
plt.xlabel("Freq (GHz)")
plt.ylabel("Flux Density (mJy)")
plt.show()
</source>
-->

=== A publication-quality figure: plot the spectrum against velocity ===

Nice, now let's make a velocity axis instead. We'll use the optical convention to convert frequency to velocity; the radio convention is offered as an alternative but is commented out. 

<source lang="python">
# In CASA
%cpaste

# Press Enter or Return, then copy/paste the following:
import matplotlib.pyplot as plt
fRest = cubeHead['restfreq'][0]
velocity = 299792.458 * (fRest/freqSpec - 1.0) # optical convention, km/s
# velocity = 299792.458 * (1.0 - freqSpec/fRest) # radio convention, km/s
--
</source>

<!--
<source lang="python">
fRest = cubeHead['restfreq'][0]
velocity = 299792.458 * (fRest/freqSpec - 1.0) # optical convention, km/s
# velocity = 299792.458 * (1.0 - freqSpec/fRest) # radio convention, km/s
</source>
-->

[[File:MKN6_Absorption_Spectrum_Publication_Ready.png|thumb|Publication quality spectrum of MKN 6.]]

Plot the spectrum again, but this time try for publication-ready quality (see figure at right). This example is somewhat more clever to highlight the use of TeX typography in <tt>matplotlib</tt>.  Note that you will have to have a compatible version of LaTeX installed for this to be successful.

<source lang="python">
# In CASA
# set up fonts
%cpaste


# Press Enter or Return, then copy/paste the following:
from matplotlib import rc
rc('font',**{'family':'serif','serif':['Times']})
rc('text', usetex=True)
plt.clf() # clear the plot (figure)
plt.plot(velocity, cubeSpec['data'], 'k-')
# the r prefix to the following format strings prevents \ from triggering an escape sequence, 
# like \n for line feed
plt.xlabel(r"$cz$ (km s$^{-1}$)",fontsize=16) 
plt.ylabel(r"$S_{\nu}$ (mJy)",fontsize=16)
plt.ylim(-0.023,0.005)
plt.savefig("mkn6-abs.eps") # save the figure as encapsulated postscript, suitable for journals
plt.savefig("mkn6-abs.png") # save the figure as a bitmap, better suited for astro-ph
--
</source>

<!--
<source lang="python">
# set up fonts
from matplotlib import rc
rc('font',**{'family':'serif','serif':['Times']})
rc('text', usetex=True)
plt.clf() # clear the plot (figure)
plt.plot(velocity, cubeSpec['data'], 'k-')
# the r prefix to the following format strings prevents \ from triggering an escape sequence, 
#     like \n for line feed
plt.xlabel(r"$cz$ (km s$^{-1}$)",fontsize=16) 
plt.ylabel(r"$S_{\nu}$ (mJy)",fontsize=16)
plt.ylim(-0.023,0.005)
plt.savefig("mkn6-abs.eps") # save the figure as encapsulated postscript, suitable for journals
plt.savefig("mkn6-abs.png") # save the figure as a bitmap, better suited for astro-ph
</source>
-->

=== Save the spectrum to a text file ===

It's handy to have a simple, columnar text file to inspect and analyze later. Or perhaps you might want to analyze the spectrum outside of CASA. Numpy makes it easy to save and restore text files.

<source lang="python">
# In CASA
# First, combine the spectrum into a single array
%cpaste

# Press Enter or Return, then copy/paste the following:
spectrum = numpy.vstack((velocity, cubeSpec['data']))
numpy.savetxt("mkn6.txt", numpy.transpose(spectrum))
--
</source>

<!--
<source lang="python">
# First, combine the spectrum into a single array
spectrum = numpy.vstack((velocity, cubeSpec['data']))
numpy.savetxt("mkn6.txt", numpy.transpose(spectrum))
</source>
-->

Here's a listing of the first few lines of ''mkn6.txt''. 

<pre>
5.367264765747175261e+03 -2.500247210264205933e-04
5.388621623358961187e+03 -2.288779069203883410e-04
5.409981470534922664e+03  4.886248498223721981e-04
5.431344307902921173e+03  3.546753214322961867e-05
5.452710136090820924e+03 -1.000610413029789925e-03
5.474078955726684399e+03 -2.162145450711250305e-03
5.495450767438905132e+03 -1.413836842402815819e-03
5.516825571855946691e+03 -9.065763442777097225e-04
5.538203369606402703e+03 -3.753502853214740753e-03
5.559584161319202394e+03 -1.479859836399555206e-02
</pre>

Here's how to restore that text file back into CASA (or python, for that matter).

<source lang="python">
# In CASA
%cpaste

# Press Enter or Return, then copy/paste the following:
# import numpy # if you haven't already
spectrum = numpy.loadtxt("mkn6.txt",dtype="float")
x = spectrum[:,0]
y = spectrum[:,1]
--
</source>

{{Checked 5.5.0}}

<!--
<source lang="python">
# import numpy # if you haven't already
spectrum = numpy.loadtxt("mrk6.txt",dtype="float")
x = spectrum[:,0]
y = spectrum[:,1]
</source>
-->