View source for VLA high frequency Spectral Line tutorial - IRC+10216-CASA6.2.0

[[Category:EVLA]][[Category:Calibration]][[Category:Spectral Line]]


'''Please use CASA 6.2.0 for this tutorial'''  

== Overview ==
[[Image:irc10216_dust.jpg|200px|thumb|Figure 1: VLT V-band image of IRC+10216 showing dust rings out to a radius of 90" by Le&atilde;o et al. (2006, A&A, 455, 187).]]

This tutorial describes the data reduction for two spectral lines observed toward the Asymptotic Giant Branch (AGB) star [https://en.wikipedia.org/wiki/CW_Leonis IRC+10216]. 
In this 36 GHz VLA observation, an 8 MHz wide subband (spectral window, abbreviated as ''spw'' in CASA) using 64 channels was observed in each of two 1 GHz wide basebands, with one of the subbands centered on the HC3N line and the other on the SiS line.  The raw data were loaded into CASA, where zero visibility values and shadowed data were flagged/removed (note that this already may be done by the [https://go.nrao.edu/arch-info NRAO Data Archive]).  The data were then '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.manipulation.split.html#casatasks.manipulation.split split]''' to only keep the antennas delivering good quality Ka-band data (i.e., antennas without a Ka-band receiver and X-band pointing data was discarded).  Furthermore, the '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.manipulation.split.html#casatasks.manipulation.split split]''' operation was used to average the original 1-second visibility integrations to 10 seconds.  The latter produces a significantly smaller data set for processing of this tutorial, though should not be done in general on any other data set.

IRC+10216 (CW Leo, see Figure 1) is a 16th magnitude star in the optical, but the brightest star in the sky at a wavelength of 5 micrometer (60 THz). It was discovered during the first survey of the infrared sky, carried out by Bob Leighton and Gerry Neugebauer in 1965.  This Asymptotic Giant Branch star is a Mira-type variable going through prodigious episodic mass loss.  The dust condensed from the atmosphere during the mass loss is responsible for the millimeter and infrared emission; the radio continuum emission emerges from the stellar photosphere.  Molecules form along with the dust, and a steady state chemistry occurs in the dense inner regions (Tsuji 1973 A&A 23, 411).  As the density of the circumstellar material drops during the expansion of the outflow, the chemistry freezes while the molecules continue their long coast outward into the Galaxy.  As the shell thins, ultraviolet light from the ambient Galactic radiation field penetrates the dust and initiates new chemistry in the gas.

SiS, a simple molecule created in the dense inner envelope chemistry, is photodissociated as it coasts out into the shell.  It is observed as centrally concentrated emission in this tutorial. Interestingly, in the transition at 18 GHz, the SiS line shape is much different from what is seen here.  At the extreme velocities in the profile, very bright narrow emission is seen which has been interpreted as maser emission.  More VLA observations of this line can be found in the NRAO Data Archive.

HC3N, a more complex molecule, is created by photochemistry that is triggered when atoms and pieces of molecules, destroyed by ultraviolet radiation, undergo the next phase of chemistry in the shell. In addition to its rotational modes, HC3N has many vibrational modes which may be excited.  Owing to this, it can re-radiate energy absorbed from ultraviolet radiation more effectively than some molecules with a single bond.  Eventually, it too is destroyed. However, during its brief existence, its rise to abundance in the envelope results in a ring of emission, which is what is observed in this VLA tutorial. Cordiner & Millar (2009, ApJ, 697, 68) describe a new chemical model for the shell, which also takes into account the variation of mass loss by the star.  These authors show that in addition to purely chemical effects, local gas and dust density peaks play a role in shaping the observed emission.

==Obtaining the Data==

The post-split averaged data can be downloaded from [http://casa.nrao.edu/Data/EVLA/IRC10216/day2_TDEM0003_10s_norx.tar.gz http://casa.nrao.edu/Data/EVLA/IRC10216/day2_TDEM0003_10s_norx.tar.gz] (download data size is 1.1GB, extracted it will be about 2.0 GB)

We are providing this starting data set, rather than the original data set for two reasons.  First, many of these initial processing steps can be rather time consuming (> 1 hr).  Second, while necessary, many of these steps are not fundamental to the calibration and imaging process, which is the focus of this tutorial.  For completeness, however, here are the steps that were taken from the original data set to produce the data set above.

If you want to start the data reduction from the beginning, from the NRAO Archive select the TDEM0003_sb1345754_1.55312.131578217595 data set.

* The initial Science Data Model (SDM) file was converted into a measurement set (MS).
* Basic data flagging was applied to account for shadowing of the antennas.  In the D array configuration antennas may be blocked (shadowed) by other antennas in the array, depending upon the direction and elevation of the source.
* The data were averaged from the original 1-second correlator visibility integration time to 10-second averages.  In the D configuration, the fringe rate is relatively slow and time-average smearing is less of a concern.

All of these steps can be done directly from the NRAO Archive's Download page, by selecting '''CASA MS''' as the download format (it's a good idea to also check the '''Create MS or SDM tar file''' box), checking the '''Apply flags generated during observing''' box, and setting '''Time Averaging''' to 10s. Then, as noted above in the Overview, the Ka-band data for this source was split from the larger TDEM0003 MS to create the smaller data set that this tutorial covers. The last step, splitting out the Ka-band data, must be done by the user within CASA; it cannot be done using the NRAO archive interface.

Once the download is complete, unzip and unpack the file (within its own working directory, in which then CASA is run):

<pre style="background-color:lightgrey;">
# in a terminal, outside of CASA:
tar -xzf day2_TDEM0003_10s_norx.tar.gz
</pre>

== How to Use This CASA Guide==

[[Image:listob-inp-CASA5.5.0.png|200px|thumb|Figure 2: Inputs for the task '''listobs''' as run in this tutorial]]

Please use CASA 6.2.0 for this tutorial. <!-- For internal purposes, typing ''casa -ls'' in a linux window shows the available versions and the current version and to explicitly change the current version type, ''casa -r 6.2.0-124''). -->

There are at least three different ways to interact with CASA, described in more detail in [[Getting Started in CASA]]. In this guide we provide the pseudo-interactive method for every step and the interactive method for some of the steps. Since it is possible to use both methods, take care not to run the same task with identical parameters twice using both methods.

* Interactively examining task inputs. In this mode, type ''' ''taskname'' ''' to load the task, '''inp''' to examine the inputs (see Figure 2), and '''go''' once those inputs have been set to your satisfaction. Inputs that deviate from the default value are shown in blue (instead of black for default values) and bad inputs are colored red; those red values need to be fixed before a '''go''' is meaningful. The input parameters themselves are changed one by one, e.g., ''selectdata=True''. Summaries of the inputs to various tasks used in the data reduction below are provided, to illustrate which parameters need to be set. More detailed help can be obtained on any task by typing '''help ''taskname'''''. Once a task is run, the set of inputs are stored and can be retrieved via '''tget ''taskname'''''; subsequent runs will overwrite the previous '''tget''' file. To reset a task to its default settings type, '''default taskname'''. 
<pre style="background-color: #d8bfd8;">
# Interactive CASA
<default|tget> taskname
inp 
parameter1 = value1
parameter2 = value2
(etc) 
inp   (Always double check the input parameters before running the task; fix the parameter values in red font.)  
go
</pre>

* Pseudo-interactively via task function calls. In this case, all of the desired inputs to a task are provided at once on the CASA command line. This tutorial is made up of such calls, which were developed by looking at the inputs for each task and deciding what is needed to be changed from the default values. For task function calls, only parameters that you want to be different from their defaults need to be set.
<source lang="python">
# Pseudo-interactive CASA
taskname('input parameters')
</source>

* Non-interactively via a script. A series of task function calls can be combined together into a script, and run from within CASA via '''execfile('scriptname.py')'''. This and other CASA Tutorial Guides have been designed to be extracted into a script via the script extractor by using the method described within the [[Extracting scripts from these tutorials]] page. Should you use the script generated by the script extractor for this CASA Guide, be aware that it will require some small amount of interaction related to the plotting, occasionally suggesting that you close the graphics window and hitting return in the terminal to proceed. It is in fact unnecessary to close the graphics windows (it is suggested that you do so purely to keep your desktop uncluttered).

==Initial Inspection and Flagging==

=== Observing Logs ===

For all VLA observations, the operators keep an observing log.  To access the operator logs, go to the [http://www.vla.nrao.edu/cgi-bin/oplogs.cgi observing log website] and enter a date range covering only 26 Apr 2010, then select TDEM0003. Pertinent information from this observation is repeated below:

<pre style="background-color: #E0FFFF;">
INFORMATION FROM OBSERVING LOG:
Top left:
  Date of the observation: 26-April-2010
  Antenna configuration is "D"
  Antenna 10 (ea10) is not part of the 27 element array (it is the 28th antenna undergoing routine maintenance)
Operator comments:
  There are no Ka-band receivers on ea11, ea13, ea14, ea16, ea17, ea18, ea26 (those still had to be commissioned)
  Antenna ea06 is also taken out of the array (for major maintenance)
  Antenna 15 has a focus timeout in C band (used for pointing - be careful when these solutions are applied to Ka observations for this antenna)
  Antennas ea10, ea12, ea22 do not have good baseline positions at the time of the observations (can be corrected later)
</pre>

All of that information should be kept in mind during the calibration, e.g., if unavailable antenna receivers from the list above do appear in the data, the data is just noise and should be flagged. For the data set that you have downloaded, these antennas have already been removed.  Also, keep an eye on other messages such as loss of data, subreflector problems, or any other such issues. If they were not captured in the online flagging, they should be carefully inspected and treated appropriately during calibration or flagged. 

To start, we will look at the content of the raw data. Task '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.information.listobs.html#casatasks.information.listobs listobs]''' provides almost all relevant observational parameters such as correlator setup (frequencies, bandwidths, channel number and widths, polarization products), sources, scans, scan intents, and antenna locations:

<pre style="background-color: #d8bfd8;">
# Interactive CASA
default listobs
inp
vis='day2_TDEM0003_10s_norx'
inp
go
</pre>

<source lang="python">
# Pseudo-interactive CASA
listobs(vis='day2_TDEM0003_10s_norx')
</source>

Below we have cut and pasted the output from the logger that is listed after the observed scan listing. Note that ''Field'' refers to an observed phase center, i.e., a position in the sky. A ''Source'' refers to a Field using a specific frequency setting. Note, '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.information.listobs.html#casatasks.information.listobs listobs]''' and several other tasks will now return dictionaries in the terminal as a convenient way of obtaining machine-readable output (as opposed to parsing the log, etc.).

<pre style="background-color: #fffacd;">
Fields: 4
ID   Code Name                RA               Decl           Epoch   SrcId      nRows
2    D    J0954+1743          09:54:56.823626 +17.43.31.22243 J2000   2          65326
3    NONE IRC+10216           09:47:57.382000 +13.16.40.65999 J2000   3         208242
5    F    J1229+0203          12:29:06.699729 +02.03.08.59820 J2000   5          10836
7    E    J1331+3030          13:31:08.287984 +30.30.32.95886 J2000   7           5814
Spectral Windows:  (2 unique spectral windows and 1 unique polarization setups)
SpwID  Name      #Chans   Frame   Ch0(MHz)  ChanWid(kHz)  TotBW(kHz) CtrFreq(MHz)  Corrs          
0      Subband:0     64   TOPO   36387.229       125.000      8000.0  36391.1670   RR  RL  LR  LL
1      Subband:0     64   TOPO   36304.542       125.000      8000.0  36308.4795   RR  RL  LR  LL
Sources: 10
ID   Name                SpwId RestFreq(MHz)  SysVel(km/s) 
0    J1008+0730          0     0.03639232     -0.026       
0    J1008+0730          1     0.03639232     -0.026       
2    J0954+1743          0     0.03639232     -0.026       
2    J0954+1743          1     0.03639232     -0.026       
3    IRC+10216           0     0.03639232     -0.026       
3    IRC+10216           1     0.03639232     -0.026       
5    J1229+0203          0     0.03639232     -0.026       
5    J1229+0203          1     0.03639232     -0.026       
7    J1331+3030          0     0.03639232     -0.026       
7    J1331+3030          1     0.03639232     -0.026       
Antennas: 19:
ID   Name  Station   Diam.    Long.         Lat.                Offset from array center (m)                ITRF Geocentric coordinates (m)        
	                                                           East         North     Elevation               x             y               z
0    ea01  W09       25.0 m   -107.37.25.2  +33.53.51.0       -521.9407     -332.7782       -1.1977 -1601710.017000 -5042006.928200  3554602.355600
1    ea02  E02       25.0 m   -107.37.04.4  +33.54.01.1          9.8247      -20.4292       -2.7808 -1601150.059500 -5042000.619800  3554860.729400
2    ea03  E09       25.0 m   -107.36.45.1  +33.53.53.6        506.0591     -251.8666       -3.5832 -1600715.948000 -5042273.187000  3554668.184500
3    ea04  W01       25.0 m   -107.37.05.9  +33.54.00.5        -27.3562      -41.3030       -2.7418 -1601189.030140 -5042000.493300  3554843.425700
4    ea05  W08       25.0 m   -107.37.21.6  +33.53.53.0       -432.1158     -272.1493       -1.5032 -1601614.091000 -5042001.655700  3554652.509300
5    ea07  N06       25.0 m   -107.37.06.9  +33.54.10.3        -54.0667      263.8720       -4.2292 -1601162.593200 -5041829.000000  3555095.890500
6    ea08  N01       25.0 m   -107.37.06.0  +33.54.01.8        -30.8810       -1.4664       -2.8597 -1601185.634945 -5041978.156586  3554876.424700
7    ea09  E06       25.0 m   -107.36.55.6  +33.53.57.7        236.9058     -126.3369       -2.4443 -1600951.588000 -5042125.911000  3554773.012300
8    ea12  E08       25.0 m   -107.36.48.9  +33.53.55.1        407.8394     -206.0057       -3.2252 -1600801.916000 -5042219.371000  3554706.449900
9    ea15  W06       25.0 m   -107.37.15.6  +33.53.56.4       -275.8288     -166.7451       -2.0590 -1601447.198000 -5041992.502500  3554739.687600
10   ea19  W04       25.0 m   -107.37.10.8  +33.53.59.1       -152.8599      -83.8054       -2.4614 -1601315.893000 -5041985.320170  3554808.304600
11   ea20  N05       25.0 m   -107.37.06.7  +33.54.08.0        -47.8454      192.6015       -3.8723 -1601168.786100 -5041869.054000  3555036.936000
12   ea21  E01       25.0 m   -107.37.05.7  +33.53.59.2        -23.8638      -81.1510       -2.5851 -1601192.467800 -5042022.856800  3554810.438800
13   ea22  N04       25.0 m   -107.37.06.5  +33.54.06.1        -42.5986      132.8623       -3.5431 -1601173.953700 -5041902.660400  3554987.536500
14   ea23  E07       25.0 m   -107.36.52.4  +33.53.56.5        318.0523     -164.1848       -2.6960 -1600880.570000 -5042170.388000  3554741.457400
15   ea24  W05       25.0 m   -107.37.13.0  +33.53.57.8       -210.0944     -122.3885       -2.2581 -1601377.008000 -5041988.665500  3554776.393400
16   ea25  N02       25.0 m   -107.37.06.2  +33.54.03.5        -35.6245       53.1806       -3.1345 -1601180.861480 -5041947.453400  3554921.628700
17   ea27  E03       25.0 m   -107.37.02.8  +33.54.00.5         50.6647      -39.4832       -2.7249 -1601114.365500 -5042023.153700  3554844.945600
18   ea28  N08       25.0 m   -107.37.07.5  +33.54.15.8        -68.9057      433.1889       -5.0602 -1601147.940400 -5041733.837000  3555235.956000
</pre>

Source J1008+0730, with ID 0, does not show up in the Field listing as this source did not generate any unflagged visibilities in the data set for these Ka-band observations. It was used for pointing in C-band only.

In addition to source names, antenna names, etc., CASA indexes each of these parameters starting from 0. To specify, e.g., a source field, we can use the index or the name (both are accepted). <!-- Note: &#151; is HTML code for em dash LOS:2017 but does not work in my version. ALSO m-dashes should be avoided in favor of using proper running sentences so I am removing them all -->

'''Note:''' The Rest Frequency and Systemic Velocity are wrong in the listobs log by a factor 10^6 and 1000, respectively, given the quoted units (MHz) and (km/s). This was due to a temporary error in the VLA Observation Preparation Tool (OPT) that has subsequently been fixed.  Because the sky frequencies are correct, and we set the rest frequency explicitly later in the deconvolution stage, this does not present a problem for the data reduction. 

<pre style="background-color: #E0FFFF;">
Summary of Observing Strategy
Gain Calibrator: J0954+1743, field id=2
Bandpass Calibrator: J1229+0203, field id=5
Flux Calibrator: J1331+3030 (3C286), field id=7
Target: IRC+10216,  field id=3
Ka-band spws = 0,1
</pre> 

[[Image:CASA5.7.0-antpos.png|200px|thumb|Figure 3: Antenna locations from running plotants ]]
Next, we will look at a graphical plot of the antenna locations using '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.visualization.plotants.html#casatasks.visualization.plotants plotants]''' (see Figure 3) and save an image copy to disk for later reference. This will be useful for selecting a reference antenna. Typically a good choice is an antenna close to the center of the array, but not too close to the center that it suffers from shadowing (as this could be the case in D-array configuration). Unless it shows problems after inspection of the data, we provisionally choose ea02 as our reference antenna.

[[Image:casa5.4.0-f04_elevationvstime.png|200px|thumb|Figure 4: Elevation as a function of time (after selecting colorize by field).]]

<pre style="background-color: #d8bfd8;">
# Interactive CASA
default plotants
inp
vis='day2_TDEM0003_10s_norx'
figfile='ant_locations.png'
inp
go
</pre>

<source lang="python">
# Pseudo-interactive CASA
plotants(vis='day2_TDEM0003_10s_norx',figfile='ant_locations.png')
</source>

Next we will look at the elevation as a function of time for all sources (see Figure 4). While it's not the case for these data, you may want to flag if the elevation is very low (usually at the start or end of an observation). Also, how near in elevation your flux calibrator is to your target will impact your ultimate absolute flux density calibration accuracy. 

Unfortunately, the target and flux density calibrator are not particularly well-matched in elevation for this observation, as you can see by plotting the elevation for each source (indicated by a different color). We will be using data in spw 0, channels 4 to 60 for this plot as given by the parameter ''spw = '0:4~60' ''; the general CASA selection syntax is described in the [https://casadocs.readthedocs.io/en/v6.2.0/notebooks/visibility_data_selection.html?highlight=visibilty%20data%20selection Measurement Selection Syntax]. Note that because we have set the channel averaging parameter to 64 (as a convenience to always average all channels in a spw), but only select 57 channels (4 to 60), CASA will issue a warning about this in the terminal. However, as the extra channels are not selected to contain additional data, this is nothing to be concerned about. RR and LL are the parallel hand correlations:

<pre style="background-color: #d8bfd8;">
# Interactive CASA
default plotms
inp
vis='day2_TDEM0003_10s_norx'
xaxis='time'
yaxis='elevation'
spw='0:4~60'
correlation='RR,LL'
avgchannel='64'
coloraxis='field'
inp
go
</pre>

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='day2_TDEM0003_10s_norx',
       xaxis='time',yaxis='elevation',correlation='RR,LL',
       avgchannel='64',spw='0:4~60', coloraxis='field')
</source>

Thus, we are strongly dependent on the opacity and gaincurve corrections to get the flux scale right for these data. (This is something to keep in mind when planning observations.) We will deal with this after the initial inspection.

[[Image:casa5.4.0-f05_raw_amptime.png|200px|thumb|Figure 5: Result of plotms]]
[[Image:casa5.4.0-f06_flag_zero.png|200px|thumb|Figure 6: Zooming in and marking region (hatched box)]]

Next, we will look at all the source amplitudes as a function of time (see Figure 5):

<pre style="background-color: #d8bfd8;">
# Interactive CASA
yaxis='amp'
inp
go
</pre>

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='day2_TDEM0003_10s_norx',
       xaxis='time',yaxis='amp',correlation='RR,LL',
       avgchannel='64',spw='0:4~60', coloraxis='field')
</source>

Now zoom in on the region very near zero amplitude for sources J0954+1743 and IRC+10216 between 3h30 and 4h30 (see Figure 6). To zoom, select the 
Zoom tool [[File:casaplotms-zoom.png]] in the lower left corner of  the '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.visualization.plotms.html#casatasks.visualization.plotms plotms]''' GUI, then you can left click to draw a box (right click to unzoom). 
Look for the low values (you may want to zoom in a few times to really see the suspect points clearly). Now use the Mark Region [[Image:MarkRegionsButton.png]] and Locate [[File:casaplotms-locate-tool.png]] buttons (positioned along the bottom of the GUI) to see which antenna is causing problems. The output will be shown in the logger. Since all the located baselines include ea12, this is likely the responsible antenna.

<pre style="background-color: #98FB98;">
IMPORTANT NOTES ON PLOTMS:

* When using the locate button, it is important to have only selected a modest number 
of points with the mark region tool (see example of marked region in the thumbnail), 
otherwise the response will be very slow and possibly hang the tool 
(all of the information will be output to your terminal window, not the logger). 

* Throughout the tutorial, when you are done marking/locate, use the Clear Regions 
tool to get rid of the marked box before plotting other things. 

* After flagdata command flagging, you have to force a complete reload of the cache 
to look at the same plot again with the new flags applied. To do this, either check 
the "reload" box in the lower left before clicking Plot, or do Shift+Plot.

* Fields in the plotms GUI tabs can be altered and the data replotted
without having to close and restart plotms. After making changes, check the
Reload box and then click the Plot button adjacent to the Reload box.


</pre> 

Now click the clear region [[File:casaplotms-clear-region.png]] button, then go back to the zoom button to zoom in further to note exactly what the time range is: 03:41:00~04:10:00.

Check the other sideband by changing spw to 1:4~60. You will have to
re-zoom. If you have trouble, click on the Mark icon and then back to
zoom.  In spw 1, ea07 is bad from the beginning until after the next
pointing run: within the Data tab of plotms, set the timerange parameter to 03:21:40~04:10:00.  To see this, compare the amplitudes when antenna is set to 'ea07' and when it is set to one of the other antennas, such as 'ea08'.

If you set antenna to 'ea12' and zoom in on this initial time range, you can also see that ea12 is bad during the same time range as for spw 0.  You can also see this by entering '!ea07' in the antenna parameter, which removes ea07 from the plot (in CASA, ! deselects). We will also flag some bad data from antenna ea08.

We will now use '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.flagging.flagdata.html#casatasks.flagging.flagdata flagdata]''' to flag the three bad antennas (ea12, ea07, and ea08) for the appropriate time range and spw:

<pre style="background-color: #d8bfd8;">
# Interactive CASA
default flagdata
inp
vis='day2_TDEM0003_10s_norx'
mode='list'
inpfile=["field='2,3' antenna='ea12' timerange='03:41:00~04:10:00'","field='2,3' antenna='ea07,ea08' timerange='03:21:40~04:10:00' spw='1'"]
inp
go
</pre>

<source lang="python">
# Pseudo-interactive CASA
flagdata(vis='day2_TDEM0003_10s_norx',
         mode='list', 
         inpfile=["field='2,3' antenna='ea12' timerange='03:41:00~04:10:00'",
                  "field='2,3' antenna='ea07,ea08' timerange='03:21:40~04:10:00' spw='1'"])
</source>


Task '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.flagging.flagdata.html#casatasks.flagging.flagdata flagdata]''' in ''mode='list' ''can take several commands separated by commas. The data selection will be connected by AND for each string. 

Note that because the chosen time range is limited to fields 2 and 3, the field parameter is not really
needed but the above is a nice example of the task syntax.

Now remove the !ea07 from the antenna parameter and replot both spw, zooming in to
be sure that all obviously low points are gone. Also zoom in and
check 3C286 (J1229+0203 is already obvious because it is so bright). 

[[Image:casa5.4.0-f07_IRC10216_uvdist_raw.png|200px|thumb|Figure 7: Amplitude vs. uv-distance for IRC+10216, both spw]]

Next, look more closely at IRC+10216:

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='day2_TDEM0003_10s_norx',field='3',
       xaxis='time',yaxis='amp',correlation='RR,LL',
       avgchannel='64',spw='0~1:4~60', coloraxis='spw')
</source>

We can see there are some noisy high points. So now try:

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='day2_TDEM0003_10s_norx',field='3',
       xaxis='uvdist',yaxis='amp',correlation='RR,LL',
       avgchannel='64',spw='0~1:4~60', coloraxis='spw')
</source>

Now we can see that most of the high points on IRC+10216 are due to large scale emission on short baselines (see Figure 7), but there is still some noisy stuff. For a target like this with extended emission it's best to wait until later to decide what to do about it.  We will not be able to get adequate calibration for antennas that are truly bad (even if they don't stand out here), so these will be obvious later.

== The Observing Log, Antenna Position Corrections, Opacities, Gaincurves, Antenna Efficiencies and other Calibration "Priors"==

=== Antenna Position Corrections ===

When antennas are moved in the array, their position may not be known very well until a baseline solution has been determined at the telescope. In that case, offset positions will be published a few days after the observations and stored in a [http://www.vla.nrao.edu/astro/archive/baselines/ machine readable table] that CASA can access. This step is to determine if there are offsets and apply them if present. As mentioned in the observing log above, antennas ea10, ea12, and ea22 do not have good baseline positions. 
Antenna ea10 was not in the array, but for the other two antennas we need to check for
any improved baseline positions that were derived after the observations were taken.
In CASA, these corrections can be derived by '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.gencal.html#casatasks.calibration.gencal gencal]''' for the VLA. The parameter ''caltype='antpos' ''will use antenna position offset values and create a calibration table for the delay changes that this correction involves. If the 'antenna' parameter is not specified, '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.gencal.html#casatasks.calibration.gencal gencal]''' will query the VLA webpages for the offsets ([http://www.vla.nrao.edu/astro/archive/baselines/ VLA Baseline Corrections page]).

<pre style="background-color: #d8bfd8;">
# Interactive CASA
default gencal
inp
vis='day2_TDEM0003_10s_norx'
caltable='antpos.cal'
caltype='antpos'
inp
go
</pre>

<source lang="python" style="color: red">
# Pseudo-interactive CASA
gencal(vis='day2_TDEM0003_10s_norx',caltable='antpos.cal',
       caltype='antpos',
       antenna='')
</source>

The above '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.gencal.html#casatasks.calibration.gencal gencal]''' returns the antenna offsets in the logger:

<pre style="background-color: #fffacd;">
Determine antenna position offests from the baseline correction database
offsets for antenna ea01 :  0.00000   0.00300   0.00000
offsets for antenna ea02 : -0.00080   0.00000   0.00000
offsets for antenna ea03 : -0.00280   0.00000   0.00000
offsets for antenna ea05 :  0.00000   0.00280   0.00000
offsets for antenna ea12 : -0.01000   0.00450  -0.00170
offsets for antenna ea22 : -0.02570   0.00270  -0.01900
offsets for antenna ea23 : -0.00140   0.00000   0.00000
offsets for antenna ea24 : -0.00150   0.00000   0.00000
offsets for antenna ea27 :  0.00000   0.00190  -0.00160
</pre>

Apparently, more antennas than the ones noted by the operator in the log needed, and now benefit from updated positions.  However, ea12 and ea22 clearly have the largest corrections.

'''Please note: '''If you are reducing VLA data taken before 2010 March 01, you will need to set '' caltype='antposvla' ''. Before this date the automated lookup will not work and you will have to specify the antenna positions explicitly.

=== Gaincurves and Antenna Efficiencies ===


The gaincurve describes how each antenna behaves as a function of elevation, for each receiver band. Currently, only gaincurves for the VLA are available<!-- BAD LINK - HIDDEN FOR TIME BEING (see [http://www.vla.nrao.edu/astro/guides/vlbivla/current/node18.html] for the incorporated models)-->. This option should not be used with any other telescopes. 

Gaincurves are generated into a separate calibration table using the task '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.gencal.html#casatasks.calibration.gencal gencal]''' with the '' caltype='gceff' ''. This option will also include antenna efficiencies in the resulting calibration table:

<source lang="python">
# Pseudo-interactive CASA
gencal(vis='day2_TDEM0003_10s_norx',caltable='gaincurve.cal',
       caltype='gceff')
</source>

<!-- ''' Note that this method replaces the ''gaincurve=True'' option in the calibrations tasks in CASA 3.4 and earlier ''' -->

=== Opacities ===

Opacity corrections are typically only applied for high frequency (e.g., Ku, K, Ka, & Q band) observations, where the atmospheric opacity affects the flux density scale if the sources (including calibrators) are observed at different elevations. For lower frequencies the effect is marginal, and opacity correction is usually skipped.
The opacity of the observation can be computed from a seasonal model and/or weather station information. After the zenith opacities are derived, they will be recomputed for the correct elevation of the data automatically using <math>e^{(-\csc[el]\tau_z)}</math> in '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.gaincal.html#casatasks.calibration.gaincal gaincal]''', '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.applycal.html#casatasks.calibration.applycal applycal]''', '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.bandpass.html#casatasks.calibration.bandpass bandpass]''', etc. 

To start, we want to plot the opacity of the atmosphere at 36 GHz at the time this observation was taken so it can be corrected in subsequent calibrations. We will use the task '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.visualization.plotweather.html#casatasks.visualization.plotweather plotweather]''' to plot the weather conditions at the time of the observation and calculate the atmospheric opacities based on these data in combination with a seasonal model that contains long-term statistics at the VLA site (see Figure 8). The parameter 'seasonal_weight=0.5' gives both equal weights: 

[[Image:CASA6.2.0-plotweather.png|200px|thumb|right|Figure 8: plotweather output]]

We will be running '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.visualization.plotweather.html#casatasks.visualization.plotweather plotweather]''' in a way that will assign the opacity list (one entry for each spw in ascending order) to the variable myTau:

<source lang="python">
# Pseudo-interactive CASA
myTau = plotweather(vis='day2_TDEM0003_10s_norx', doPlot=True)
</source>

With this input, the task will:

* return the opacity values to the logger;

<pre style="background-color: #fffacd;">
SPW : Frequency (GHz) : Zenith opacity (nepers)
 0  :   36.387  :  0.038
 1  :   36.305  :  0.038
wrote weather figure: day2_TDEM0003_10s_norx.plotweather.png
</pre> 

* create a file '''day2_TDEM0003_10s_norx.plotweather.png''' (in the current working directory) with the elevation of the sun, wind speed and direction, temperature, and precipitable water vapor (PWV) as functions of time over the observation. View this file with your preferred image viewer like gthumb, xv, or Preview;

* and assign the myTau variable to the list of opacities per spectral window. 

<source lang="python">
# Pseudo-interactive CASA 
# Type the variable by itself
myTau
</source>
It echoes the following output, the values set per spectral window:
<pre style="background-color: #fffacd;">
  Out[..]: [0.038307558882357123, 0.038092235934944034]
</pre> 

Note that you may get a SEVERE error warning when running '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.visualization.plotweather.html#casatasks.visualization.plotweather plotweather]''' concerning the leap second being out of date:

<pre style="background-color: #fffacd;">
SEVERE MeasTable	Leap second table TAI_UTC seems out-of-date.
SEVERE MeasTable	Until the table is updated (see the CASA documentation or your system admin),
SEVERE MeasTable	times and coordinates derived from UTC could be wrong by 1s or more.
</pre> 

This happens when the CASA repository is not up to date. This repository can be updated manually by following the instructions found at:

[https://casadocs.readthedocs.io/en/v6.2.0/notebooks/external-data.html?highlight=casa%20data%20repository#Updating-via-repository-copy CASA Data Respository]. For the purposes of this tutorial you do not need to update the CASA repository, but keep in mind that you may need to do this on your own data.



We can now create a calibration table for the opacities via '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.gencal.html#casatasks.calibration.gencal gencal]''' with the'' calmode='opac' ''parameter; we can either input the opacities directly or use the myTau variable:

<source lang="python">
# Pseudo-interactive CASA
gencal(vis='day2_TDEM0003_10s_norx',caltable='opacity.cal',
       caltype='opac',
       spw='0,1',
       parameter=myTau)
</source>
<!-- ''' Note that this method replaces the ''opacity'' option in the calibrations tasks in CASA 3.4 and earlier ''' -->

==Set Up the Model for the Flux Calibrator==

Next, we set the model for the flux calibrator. Depending on your observing frequency and angular resolution you can do this several ways. In the past, we typically used a point source (constant flux) model for the flux calibrator, possibly with a uvrange cut-off if necessary. More recently, model images for the most common flux calibrators have been made available for use in cases where the sources are somewhat resolved. This is most likely to be true at higher frequencies and at higher resolutions (more extended arrays). 

The VLA flux calibrator models that are available in CASA can be checked through '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.imaging.setjy.html#casatasks.imaging.setjy setjy]''':

<pre style="background-color: #d8bfd8;">
# Interactive CASA
default setjy
vis='day2_TDEM0003_10s_norx'
listmodels=True
inp
go
</pre>

<source lang="python">
# Pseudo-interactive CASA
setjy(vis='day2_TDEM0003_10s_norx',listmodels=True)
</source>

The terminal (not the logger) will now show the models, e.g., 3C286_C.im, 3C48_K.im, etc. Task '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.imaging.setjy.html#casatasks.imaging.setjy setjy]''' will search in the current working directory for images that may contain models, as well as in a CASA directory where known calibrator models are stored.

We will pick the Ka ('A') band model of 3C286: 3C286_A.im.
Task '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.imaging.setjy.html#casatasks.imaging.setjy setjy]''' scales the total flux in the model image to that appropriate for your individual spectral window frequencies according to the calibrator's flux and reports this number to the logger. It is a good idea to save this information for your records.

<source lang="python">
# Pseudo-interactive CASA
setjy(vis='day2_TDEM0003_10s_norx',field='7',spw='0~1',scalebychan=True,model='3C286_A.im')
</source>

The logger output for each spw is:
<pre style="background-color: #fffacd;">
J1331+3030 (fld ind 7) spw 0  [I=1.7469, Q=0, U=0, V=0] Jy @ 3.6387e+10Hz, (Perley-Butler 2017)
J1331+3030 (fld ind 7) spw 1  [I=1.7499, Q=0, U=0, V=0] Jy @ 3.6305e+10Hz, (Perley-Butler 2017)
Using model image /home/casa/data/distro/nrao/VLA/CalModels/3C286_A.im
Scaling spw(s) [0, 1]'s model image by channel to  I = 1.74988, 1.74824, 1.7466 Jy @(3.63045e+10, 3.63499e+10, 3.63953e+10)Hz (LSRK) for visibility prediction (a few representative values are shown).
</pre>

The absolute fluxes for the frequencies of all spectral window channels have now been determined, and we can proceed to the bandpass and complex gain calibrations.

==Bandpass and Delay==

Before determining the bandpass solution, we need to inspect phase and amplitude
variations with time and frequency on the bandpass calibrator to
decide how best to proceed. We limit the number of antennas to make
the plot easier to see. We chose ea02 for the reference antenna as it seems like a good
candidate. Let's try a single baseline to antenna ea23:

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='day2_TDEM0003_10s_norx',field='5',
       xaxis='channel',yaxis='phase',correlation='RR',
       avgtime='1e8',spw='0:4~60',antenna='ea02&ea23')
</source>
[[Image:casa5.4.0-f09_nobandpass_phase.png|200px|thumb|Figure 9: Phase as a function of channel for ea02.]]

The phase variation is modest ~10 degrees.  Now expand to all baselines that include ea02, and add an extra dimension of color (see Figure 9):

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='day2_TDEM0003_10s_norx',field='5',
       xaxis='channel',yaxis='phase',correlation='RR',
       avgtime='1e8',spw='0:4~60',antenna='ea02', coloraxis='antenna2')
</source>
[[Image:casa5.4.0-f10_nobandpass_phasetime.png|200px|thumb|Figure 10: Phase as a function of time for all baselines with antenna ea02.]]
From this we can see that the phase variation across the bandpass for each baseline to antenna ea02 is modest. Next check LL, and then spw=1 for both correlations. Also check
other antennas, if you like.

Now look at the phase as a function of ''time'' (note that there is a 360 degree wrap ambiguity on this baseline):

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='day2_TDEM0003_10s_norx',field='5',
       xaxis='time',yaxis='phase',correlation='RR',
       avgchannel='64',spw='0:4~60',antenna='ea02&ea23')
</source>

Expand to all antennas with ea02:

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='day2_TDEM0003_10s_norx',field='5',
       xaxis='time',yaxis='phase',correlation='RR',
       avgchannel='64',spw='0:4~60',antenna='ea02', coloraxis='antenna2')
</source>

We can see the phase variations are smooth (see Figure 10), but do vary significantly over the 5 minutes of observation and in most cases by a few 10s of degrees (zoom in to see this better).

The conclusion from this investigation is that, especially for high frequency observations, we need to correct the phase variations with time before solving for the bandpass to prevent decorrelation of the vector averaged bandpass solution. Since the phase variation as a function of channel is modest (as noted in Figure 9), we can average over several (or all) channels to increase the signal to noise of the phase vs. time solution. If the phase variation, as a function of channel is larger, we may need to restrict averaging to only a few channels to prevent introducing delay-based closure errors as can happen from averaging over non-bandpass corrected channels with large phase variations.

First, let's work on the delay calibration. The delay is the linear slope of phase across frequency. From the plot we saw that it was rather modest (the phases over frequency per baseline were aligned almost horizontal) and the bandpass calibration will certainly take care of it. Nevertheless, it is best to derive a delay calibration first and then calculate the bandpass. The delay calibration is an antenna based calibration solution and it can be derived in '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.gaincal.html#casatasks.calibration.gaincal gaincal]''' with parameter'' gaintype='K' ''. 

<pre style="background-color: #d8bfd8;">
# Interactive CASA
default gaincal
inp
vis='day2_TDEM0003_10s_norx'
caltable='delays.cal'
field='5'
refant='ea02'
gaintype='K'
gaintable=['antpos.cal','gaincurve.cal','opacity.cal']
inp
go
</pre>

<source lang="python">
# Pseudo-interactive CASA
gaincal(vis='day2_TDEM0003_10s_norx', caltable='delays.cal', field='5', 
        refant='ea02', gaintype='K', gaintable=['antpos.cal','gaincurve.cal','opacity.cal'])
</source>

(Note, the warning "No VLATrDelCorr keyword in the antpos caltable; turning trop delay correction OFF" can be ignored for this dataset. However, if your data was taken between 2016 August 09 to 2016 November 14, it may be affected by the [https://science.nrao.edu/facilities/vla/data-processing/vla-atmospheric-delay-problem VLA Atmospheric Delay Problem], and the logger will report that the correction was turned ON.) 


We use the strong bandpass calibrator for the solution. It will be extrapolated in time to all observations.

This call also exemplifies how the CASA calibration table system works. Calibration tables are incremental. So we supply the previous 'antpos.cal', 'gaincurve.cal', and 'opacity.cal' tables as input to '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.gaincal.html#casatasks.calibration.gaincal gaincal]'''. The following calibration step will require the parameter'' gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delay.cal'] ''and so on. 


Now we proceed to the actual bandpass calibration. 
The first step is to correct for decorrelation. Since the bandpass calibrator is quite strong we do the phase-only
solution on the integration time of 10 seconds (parameter ''solint='int' '').

<source lang="python">
# Pseudo-interactive CASA
gaincal(vis='day2_TDEM0003_10s_norx',caltable='bpphase.gcal',
        field='5',spw='0~1:20~40',
        refant='ea02',calmode='p',solint='int',minsnr=2.0,
        gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delays.cal'])
</source>

Next we will plot the solutions with '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.visualization.plotms.html#casatasks.visualization.plotms plotms]''' (see Figure 11).
[[Image:bpphase-PhasevsTime-CASA5.5.0.png|200px|thumb|Figure 11: Phase only calibration before bandpass. The four lines of dots are both polarizations in both spw.]]

<pre style="background-color: #d8bfd8;">
# Interactive CASA
default plotms
inp
vis='bpphase.gcal'
gridrows=3
gridcols=3
xaxis='time'
yaxis='phase'
iteraxis='antenna'
plotrange=[0,0,-180,180]
coloraxis='corr'
inp
go
</pre>

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='bpphase.gcal',gridrows=3,gridcols=3,xaxis='time',yaxis='phase',
       iteraxis='antenna',coloraxis='corr',plotrange=[0,0,-180,180])
</source> 

These solutions will appear in the CASA plotter GUI. If you closed it after plotting the antennas positions above, it should reopen. If it is still open from before, the new plots should just appear. After you are done looking at the first set of plots, push the Next button on the GUI to see the next set of antennas. When using plotcal, the plot will not clear unused panels, so the last page may show some of the older panels that are not overwritten with new data; however plotms does its own thing and will rescale the remaining panels of the gridrows-by-gridcols plot to look differently from the first few (neater) panels.

Now we can apply this phase solution on-the-fly while determining the bandpass solutions on the timescale of the bandpass calibrator scan (note parameter '' solint='inf' '', which is the default). 

<pre style="background-color: #d8bfd8;">
# Interactive CASA
default bandpass
inp
vis='day2_TDEM0003_10s_norx'
caltable='bandpass.bcal'
field='5'
refant='ea02'
solnorm=False
gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delays.cal','bpphase.gcal']
inp
go
</pre>

<source lang="python">
# Pseudo-interactive CASA
bandpass(vis='day2_TDEM0003_10s_norx',caltable='bandpass.bcal',field='5',
        refant='ea02',solint='inf',solnorm=False,
        gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delays.cal','bpphase.gcal'])
</source>

'''A few words about solint and combine:'''

The use of parameter'' solint='inf' ''in '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.bandpass.html#casatasks.calibration.bandpass bandpass]''' will derive one bandpass solution for the whole J1229+0203 scan. Note that if there had been two observations of the bandpass calibrator, this command would have combined the data from both scans to form one bandpass solution, because the default of the combine parameter for '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.bandpass.html#casatasks.calibration.bandpass bandpass]''' is ''combine='scan' ''. To solve for one bandpass for each bandpass calibrator scan we will use ''combine=' ' ''in the bandpass call. In all calibration tasks, regardless of solint, scan boundaries are only crossed when ''combine='scan' ''. Likewise, field (spw) boundaries are only crossed if parameter ''combine='field' '' or parameter ''combine='spw' '', the latter two are not generally good ideas for bandpass solutions.  

Plot the solutions, amplitude (see Figure 12) and phase (see Figure 13):
[[Image:CASA5.7.0-Fig12-bandpass-ampsol.jpg|200px|thumb|Figure 12: Amplitude Bandpass solutions (the different colors are the RR and LL polarization products)]]
[[Image:CASA5.7.0-Fig13-bandpass-PhaseVsChan.jpg|200px|thumb|Figure 13: Phase Bandpass solutions]]

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='bandpass.bcal',gridrows=3,gridcols=3,xaxis='chan',
       yaxis='amp',iteraxis='antenna',coloraxis='corr')

plotms(vis='bandpass.bcal',gridrows=3,gridcols=3,xaxis='chan',
       yaxis='phase',iteraxis='antenna',coloraxis='corr')
</source>

Note the scale for the phases, they are very small since the delays were already taken out in a previous calibration step (stored in delays.cal and used when it is included in gaintable).

This step isn't necessary from a calibration perspective, but we can go ahead with checking the bandpass calibration on the bandpass calibrator by running '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.applycal.html#casatasks.calibration.applycal applycal]'''.  In the future, we hope to plot corrected data on-the-fly without this '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.applycal.html#casatasks.calibration.applycal applycal]''' step. Subsequent applycals will overwrite this one, so no need to worry.

[[Image:CASA5.7.0-F14-applybandpass-phase.jpg|200px|thumb|Figure 14: Phase as a function of channel, plotting the corrected data (after Custom and adjusting "Style" to 3 in the Display tab on the far left.)]]

<pre style="background-color: #d8bfd8;">
# Interactive CASA
default applycal
inp
vis='day2_TDEM0003_10s_norx'
field='5'
gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delays.cal','bandpass.bcal']
gainfield=['','','','5','5']
calwt=False
inp
go
</pre>

<source lang="python">
# Pseudo-interactive CASA
applycal(vis='day2_TDEM0003_10s_norx',field='5',
        gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delays.cal','bandpass.bcal'],
        gainfield=['','','','5','5'],calwt=False)
</source>

Task '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.applycal.html#casatasks.calibration.applycal applycal]''' works like a matrix. The first entries of the ''gaintable'' and ''gainfield'' lists are to be used together, so are the second entries etc. (except for the opacity list, which is referring to spws). All will be applied to the 'field' selection (field='5' in the first line). In the above example, 'antpos.cal' (gaintable='antpos'), which is a field independent solution (gainfield=' '), is applied to field '5'. The 'bandpass.bcal' (gaintable='bandpass.bcal') that was obtained for field '5' (gainfield='5', the bandpass observation) is also applied to field '5', etc. Again, lists within the lists are fine.  

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='day2_TDEM0003_10s_norx',field='5',xaxis='channel',yaxis='phase',
       ydatacolumn='corrected',correlation='RR',
       avgtime='1e8',spw='0:4~60',antenna='ea02', coloraxis='antenna2')

plotms(vis='day2_TDEM0003_10s_norx',field='5',xaxis='channel',yaxis='amp',
       ydatacolumn='corrected',correlation='RR',
       avgtime='1e8',spw='0:4~60',antenna='ea02', coloraxis='antenna2')
</source>

Note, the phase and amplitude as a function of channel are now very flat (see Figure 14).

==Gain Calibration==

Now that we have a bandpass solution to apply, we can solve for the antenna-based phase and amplitude gain calibration. Since the phase changes on a much shorter timescale than the amplitude, at least for the higher frequencies and/or longer baselines, we will solve for them separately. If the phase changes significantly over time and if the un-corrected phase were averaged over this timescale, the amplitude would decorrelate. Note that we now also re-solve for the gain solutions of the bandpass calibrator to derive new solutions that are corrected for the bandpass shape. Since the bandpass calibrator will not be used again, while not strictly necessary, it is useful to check its calibrated flux density. We use a minimum signal-to-noise of 2 as it seems to be a good compromise for using good data without rejecting too many solutions (parameter ''minsnr = 2 ''). 

<source lang="python">
# Pseudo-interactive CASA
gaincal(vis='day2_TDEM0003_10s_norx',caltable='intphase.gcal',field='2,5,7',
        spw='0~1:4~60',solint='int',refant='ea02',minsnr=2.0,calmode='p',
        gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delays.cal','bandpass.bcal'])
</source>
[[Image:CASA5.7.0-F15-intphase-PhasevsTime.jpg|200px|thumb|Figure 15: Plot of phase solutions on an integration time.]]

Here parameter ''solint='int' ''coupled with parameter ''calmode='p' ''will derive a single phase solution for each integration (recall, we averaged the data for every 10 seconds). Note that the bandpass table is applied on-the-fly before solving for the phase solutions, however the bandpass is NOT applied to the data permanently until applycal is run later on.

Notice that several solutions are rejected due to SNR < 2 (printed to terminal). For the most part, it is only one or two solutions out of >30 so this isn't too worrying. Take note if you see large numbers of rejected solutions per integration; this is likely an indication that solint is too short for the S/N of the data.

Now we will look at the phase solution, and note the scatter within a scan time (see Figure 15):

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='intphase.gcal',gridrows=3,gridcols=3,xaxis='time',yaxis='phase',
       iteraxis='antenna',coloraxis='corr',plotrange=[0,0,-180,180])
</source>

Although ''solint='int' ''(i.e., the integration time of 10 seconds) is the best choice to apply before for solving for the amplitude solutions, it is not a good idea to use this to apply to the target. This is because the phase-scatter within a scan can dominate the interpolation between calibrator scans. Instead, we also solve for the phase on the scan time, ''solint='inf' ''(but ''combine=' ' '', since we want one solution per scan) for application to the target later on. Unlike the bandpass task, the default of the combine parameter in the '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.gaincal.html#casatasks.calibration.gaincal gaincal]''' task is ''combine=' ' ''.
[[Image:CASA5.7.0-F16-scanphase-PhasevsTime.jpg|200px|thumb|Figure 16: Plot of phase solutions on a scan time.]]

<source lang="python">
# Pseudo-interactive CASA
gaincal(vis='day2_TDEM0003_10s_norx',caltable='scanphase.gcal',field='2,5,7',
        spw='0~1:4~60',solint='inf',refant='ea02',minsnr=2.0,calmode='p',
        gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delays.cal','bandpass.bcal'])
</source>

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='scanphase.gcal',gridrows=3,gridcols=3,xaxis='time',yaxis='phase',
       iteraxis='antenna',coloraxis='corr',plotrange=[0,0,-180,180])
</source>

Note that there are no failed solutions here because of the added S/N acquired by the longer solint (see Figure 16).
Alternatively, instead of making a separate phase solution for application to the target, we can also run '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.smoothcal.html#casatasks.calibration.smoothcal smoothcal]''' to smooth the solutions derived on the integration time.

Next we apply the bandpass and parameter ''solint='int' ''phase-only calibration solutions on-the-fly to derive amplitude solutions. 
Here the use of ''solint='inf' '', not ''combine=' ' '', will result in one solution per scan interval.

<source lang="python">
# Pseudo-interactive CASA
gaincal(vis='day2_TDEM0003_10s_norx',caltable='amp.gcal',field='2,5,7',
        spw='0~1:4~60',solint='inf',refant='ea02',minsnr=2.0,calmode='ap',
        gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delays.cal','bandpass.bcal','intphase.gcal'])
</source>
[[Image:CASA5.7.0-F17-amp-gcal-PhasevsTime.jpg|200px|thumb|Figure 17: Plot of residual phase solutions on a scan time]]

Now let's look at the resulting phase solutions (see Figure 17).  Since we have taken out the phase as best we can by applying the ''solint='int' ''phase-only solution, this plot will give a good idea of the residual phase error. If we see scatter of more than a few degrees here, we should consider going back and looking for more data to flag, in particular for the times where these errors are large, or for other effects that may explain this behavior (and flag as necessary).

<source lang="python">
# Pseudo-interactive CASA       
plotms(vis='amp.gcal',gridrows=3,gridcols=3,xaxis='time',yaxis='phase',
       iteraxis='antenna',coloraxis='corr',plotrange=[-1,-1,-20,20])
</source>

Note that we have restricted the plot range so that all antennas are plotted with the same scale (instead of auto-scaled; for the X-axis the ''-1,-1'' values still allow autoscaling) (see Figure 18).  We can see that both antenna ea12 (all times) and ea23 (first 1/3 of observation) show particularly large residual phase noise.
[[Image:CASA5.7.0-F18-ampgcal-AmpvsTime.jpg|200px|thumb|Figure 18: Plot of amplitude solutions on a scan time]]

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='amp.gcal',gridrows=3,gridcols=3,xaxis='time',yaxis='amp',
       iteraxis='antenna',coloraxis='corr',plotrange=[-1,-1,0,1])
</source>

Note that the amplitude solutions for ea12 are very low, in particular when we set ''plotrange=[-1,-1,0,1]''; this is another indication that this antenna is dubious.

Next we use the flux calibrator (whose flux density was set in '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.imaging.setjy.html#casatasks.imaging.setjy setjy]''' above) to derive the flux of the other calibrators. The flux table can be written as an incremental table, just like all other calibration tables. In that case it won't replace the amp.gcal but both tables, amp.gcal as well as flux.cal need to be carried onward.  

<pre style="background-color: #d8bfd8;">
# Interactive CASA
default fluxscale
inp
vis='day2_TDEM0003_10s_norx'
caltable='amp.gcal'
fluxtable='flux.cal'
reference='7'
incremental=True
inp
go
</pre>

<source lang="python">
# Pseudo-interactive CASA
fluxscale(vis='day2_TDEM0003_10s_norx',caltable='amp.gcal',
          fluxtable='flux.cal',reference='7',incremental=True)
</source>

It is a good idea to note the derived flux densities:

<pre style="background-color: #fffacd;">
Found reference field(s): J1331+3030
Found transfer field(s):  J0954+1743 J1229+0203

Flux density for J0954+1743 in SpW=0 (freq=3.63912e+10 Hz) is: 0.220848 +/- 0.0307827 (SNR = 7.1744, N = 38)
Flux density for J0954+1743 in SpW=1 (freq=3.63085e+10 Hz) is: 0.227301 +/- 0.0298975 (SNR = 7.60268, N = 38)
Flux density for J1229+0203 in SpW=0 (freq=3.63912e+10 Hz) is: 27.2414 +/- 0.244997 (SNR = 111.191, N = 38)
Flux density for J1229+0203 in SpW=1 (freq=3.63085e+10 Hz) is: 26.5297 +/- 0.256679 (SNR = 103.358, N = 38)
</pre>

==Applycal and Inspect==

Now we apply the calibration to each source with the derived tables and the appropriate source for that particular calibration. For the calibrators, all bandpass solutions are the same as for the bandpass calibrator (id=5), and the phase and amplitude calibration as derived from their own visibility data. 

'''Note:''' In all '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.applycal.html#casatasks.calibration.applycal applycal]''' steps we set ''calwt=False''. It is very important to turn off this parameter which determines if the weights are calibrated along with the data. Data from antennas with better receiver performance and/or longer integration times should have higher weights, and it can be advantageous to factor this information into the calibration. Before the VLA upgrade, meaningful weights were available for each visibility. However, at the time of this observation during commissioning, the VLA was not yet recording the information necessary to calculate meaningful weights. Since these data weights are used at the imaging stage you can get strange results from having ''calwt=True'' when the input weights are themselves not meaningful, especially for self-calibration on resolved sources (e.g., flux density calibrator and target). 

For more recent data, the switched power information is recorded, but we currently do not recommend using this information to calculate data weights without exercising considerable caution and in particular not for this older data set.

Next we apply the antenna position, antenna gaincurve, and opacity calibration, which are independent on source. We further apply the delay and bandpass as derived for source "5" (the other sources have no solutions for this in the table), and the short-term phase, amplitude and flux scaling that was determined from each calibrator source itself.

<pre style="background-color: #d8bfd8;">
# Interactive CASA
# For the gain/phase calibrator (field '2'):
default applycal
inp
vis='day2_TDEM0003_10s_norx'
field='2'
gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delays.cal','bandpass.bcal','intphase.gcal','amp.gcal','flux.cal']
gainfield=['','','','5','5','2','2','2']
calwt=False
inp
go

# For the bandpass calibrator (field '5'):
field='5'
gainfield=['','','','5','5','5','5','5']
inp
go

# For the flux calibrator (field '7'):
field='7'
gainfield=['','','','5','5','7','7','7']
inp
go
</pre>

<source lang="python">
# Pseudo-interactive CASA
# For the gain/phase calibrator (field '2'):
applycal(vis='day2_TDEM0003_10s_norx',field='2',
        gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delays.cal','bandpass.bcal','intphase.gcal','amp.gcal','flux.cal'],
        gainfield=['','','','5','5','2','2','2'],
        calwt=False)

# For the bandpass calibrator (field '5'):
applycal(vis='day2_TDEM0003_10s_norx',field='5',
        gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delays.cal','bandpass.bcal','intphase.gcal','amp.gcal','flux.cal'],
        gainfield=['','','','5','5','5','5','5'],
        calwt=False)

# For the flux calibrator (field '7'):
applycal(vis='day2_TDEM0003_10s_norx',field='7',
        gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delays.cal','bandpass.bcal','intphase.gcal','amp.gcal','flux.cal'],
        gainfield=['','','','5','5','7','7','7'],
        calwt=False)
</source>

For the target (field id=3) we apply the bandpass (id=5), and the calibration from the gain calibrator (id=2):

<source lang="python">
# Pseudo-interactive CASA
# For the target source IRC10216 (field '3'):
applycal(vis='day2_TDEM0003_10s_norx',field='3',
        gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delays.cal','bandpass.bcal','scanphase.gcal','amp.gcal','flux.cal'],
        gainfield=['','','','5','5','2','2','2'],
        calwt=False)
</source>

Now inspect the corrected data:
[[Image:CASA5.7.0-F19-applycal_inspect.jpg|200px|thumb|Figure 19: Plot of calibrated amplitudes over time.]]

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='day2_TDEM0003_10s_norx',
       xaxis='time',yaxis='amp',ydatacolumn='corrected',
       field='5',spw='0:4~60',correlation='RR,LL',
       avgchannel='64',coloraxis='antenna1')
</source>

This plot shows some data deviating from the average amplitudes (see Figure 19). Use methods previously described to mark a region for a small number of deviant data points and click Locate. We will find that ea12 is responsible.

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='day2_TDEM0003_10s_norx',
       xaxis='time',yaxis='amp',ydatacolumn='corrected',
       field='2',spw='0:4~60',correlation='RR,LL',
       avgchannel='64',coloraxis='antenna2')
</source>

Here we see some problems with high points. Mark some regions and locate in '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.visualization.plotms.html#casatasks.visualization.plotms plotms]''' to find out which antennas and in which spws. Pay special attention to antennas that have been called out already as showing some dubious behavior.

What we find is that ea07, for which we flagged spw=1 above, is also bad for the same timerange in spw=0. This was not obvious in the raw data because spw=0 was adjusted in the on-line system by a gain attenuator, while spw=1 was not. So a lack of power on this antenna can look like very low (and obvious) amplitudes in spw=1 but not for spw=0. Looking carefully we will see that ea07 is actually pretty noisy throughout.
[[Image:CASA5.7.0-F20-ea12.jpg|200px|thumb|Figure 20: Plot of antenna ea12 by itself]]
[[Image:CASA5.7.0-F21-ea23.jpg|200px|thumb|Figure 21: Plot of antenna ea23 by itself]]

From the Locate we also find that ea12 and ea23 show some high points; to see this, re-plot baselines using each of the antennas individually (see Figures 20 and 21):

<source lang="python">
# Pseudo-interactive CASA
# Plot of ea12:
plotms(vis='day2_TDEM0003_10s_norx',
       xaxis='time',yaxis='amp',ydatacolumn='corrected',
       field='2',spw='0:4~60',correlation='RR,LL',
       avgchannel='64',antenna='ea12', coloraxis='antenna2')

# Plot of ea23:
plotms(vis='day2_TDEM0003_10s_norx',
       xaxis='time',yaxis='amp',ydatacolumn='corrected',
       field='2',spw='0:4~60',correlation='RR,LL',
       avgchannel='64',antenna='ea23', coloraxis='antenna2')
</source>

It may be a good idea to flag ea12 completely. It's just a bit noisy all around and ea23 is pretty noisy during the first scans between the initial and second pointing. Recall that these antennas had more noisy calibration solutions.

[[Image:CASA5.7.0-F22-target-uvdist.jpg|200px|thumb|Figure 22: IRC+10216 as a function of uv-distance.]]
Now let's see how the target looks. Because the target has resolved structure, it's best to look at it as a function of uvdistance. We will exclude the three antennas we already know have problems.

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='day2_TDEM0003_10s_norx',
       xaxis='uvdist',yaxis='amp',ydatacolumn='corrected',
       field='3',spw='0:4~60',correlation='RR,LL',
       avgchannel='64',antenna='!ea07;!ea12;!ea23', coloraxis='antenna2')
</source>

The color indicates that the spikes are caused by a single antenna (see Figure 22). Use Zoom, Mark, and Locate to see which one. Also look at spw=1.

The bad antenna turns out to be ea28; to confirm, re-plot with ''antenna=!ea28'':

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='day2_TDEM0003_10s_norx',field='3',ydatacolumn='corrected',
       xaxis='uvdist',yaxis='amp',correlation='RR,LL',
       avgchannel='64',spw='0:4~60',antenna='!ea07;!ea12;!ea23;!ea28', coloraxis='antenna2')
</source>

This looks much better.  To see if it's restricted to a certain time, plot the following:

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='day2_TDEM0003_10s_norx',field='3',ydatacolumn='corrected',
       xaxis='time',yaxis='amp',correlation='RR,LL',
       avgchannel='64',spw='0:4~60',antenna='ea28', coloraxis='antenna1')
</source>

Baselines with ea28 clearly show issues until about two-thirds of the way through the observation. Plot another distant antenna to compare. 

We will flag it all, since it's hanging far out on the north arm by itself.

The additional data we've identified as bad need to be flagged, and then all the calibration steps will need to be run again.

<pre style="background-color: #d8bfd8;">
# Interactive CASA
# Completely flag antennas ea07, ea12, and ea28:
default flagdata
inp
vis='day2_TDEM0003_10s_norx'
mode='list'
inpfile=["antenna='ea07,ea12,ea28'","antenna='ea23' timerange='03:21:40~04:10:00'"]
inp
go
</pre>

<source lang="python">
# Pseudo-interactive CASA
flagdata(vis='day2_TDEM0003_10s_norx',
         mode='list',
         inpfile=["antenna='ea07,ea12,ea28'",
                   "antenna='ea23' timerange='03:21:40~04:10:00'"])
</source>

==Redo Calibration after more Flagging==

After flagging, we will now repeat the same calibration as steps above.  Here, we append  '''_redo''' to the table names to distinguish them from the first round, in case we want to compare with previous versions.   

<source lang="python">
# Pseudo-interactive CASA 
gaincal(vis='day2_TDEM0003_10s_norx',caltable='bpphase_redo.gcal',
        field='5',spw='0~1:20~40',solint='int',
        refant='ea02',minsnr=2.0,calmode='p',
        gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delays.cal'])

bandpass(vis='day2_TDEM0003_10s_norx',caltable='bandpass_redo.bcal',
        field='5',
        refant='ea02',solint='inf',solnorm=True,
        gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delays.cal','bpphase_redo.gcal'])

gaincal(vis='day2_TDEM0003_10s_norx',caltable='intphase_redo.gcal',
        field='2,5,7',spw='0~1:4~60',solint='int',
        refant='ea02',minsnr=2.0,calmode='p',
        gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delays.cal','bandpass_redo.bcal'])

gaincal(vis='day2_TDEM0003_10s_norx',caltable='scanphase_redo.gcal',
        field='2,5,7',spw='0~1:4~60',solint='inf',
        refant='ea02',minsnr=2.0,calmode='p',
        gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delays.cal','bandpass_redo.bcal'])

gaincal(vis='day2_TDEM0003_10s_norx',caltable='amp_redo.gcal',
        field='2,5,7',spw='0~1:4~60',solint='inf',
        refant='ea02',minsnr=2.0,calmode='ap',
        gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delays.cal','bandpass_redo.bcal','intphase_redo.gcal'])

fluxscale(vis='day2_TDEM0003_10s_norx',caltable='amp_redo.gcal',
          fluxtable='flux_redo.cal',reference='7',incremental=True)
</source>

<pre style="background-color: #fffacd;">
Flux density for J0954+1743 in SpW=0 (freq=3.63912e+10 Hz) is: 0.238452 +/- 0.0285143 (SNR = 8.36253, N = 32)
Flux density for J0954+1743 in SpW=1 (freq=3.63085e+10 Hz) is: 0.229846 +/- 0.0290786 (SNR = 7.9043, N = 32)
Flux density for J1229+0203 in SpW=0 (freq=3.63912e+10 Hz) is: 26.3016 +/- 0.285113 (SNR = 92.2499, N = 32)
Flux density for J1229+0203 in SpW=1 (freq=3.63085e+10 Hz) is: 25.7648 +/- 0.292648 (SNR = 88.0403, N = 32)
</pre>

Feel free to pause here and remake the calibration solution plots from above, just be sure to put in the revised table names.

==Redo Applycal and Inspect==

Now apply all of the new calibrations, which will overwrite the old ones.  These commands are identical to those above, with the exception of the '''_redo''' part of each calibration filename.

<source lang="python">
# Pseudo-interactive CASA
applycal(vis='day2_TDEM0003_10s_norx',field='2',
        gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delays.cal','bandpass_redo.bcal','intphase_redo.gcal','amp_redo.gcal','flux_redo.cal'],
        gainfield=['','','','5','5','2','2','2'],
        calwt=False)

applycal(vis='day2_TDEM0003_10s_norx',field='5',
        gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delays.cal','bandpass_redo.bcal','intphase_redo.gcal','amp_redo.gcal','flux_redo.cal'],
        gainfield=['','','','5','5','5','5','5'],
        calwt=False)

applycal(vis='day2_TDEM0003_10s_norx',field='7',
        gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delays.cal','bandpass_redo.bcal','intphase_redo.gcal','amp_redo.gcal','flux_redo.cal'],
        gainfield=['','','','5','5','7','7','7'],
        calwt=False)

applycal(vis='day2_TDEM0003_10s_norx',field='3',
        gaintable=['antpos.cal','gaincurve.cal','opacity.cal','delays.cal','bandpass_redo.bcal','scanphase_redo.gcal','amp_redo.gcal','flux_redo.cal'],
        gainfield=['','','','5','5','2','2','2'],
        calwt=False)
</source>

[[Image:CASA5.7.0-F23-gaincal-corrflag.jpg|200px|thumb|Figure 23: Gain calibrator after further flagging and recalibration]]
[[Image:CASA5.7.0-F24-target-corrflag.jpg|200px|thumb|Figure 24: IRC+10216 after further flagging and recalibration (after selecting colorize by spw).]]

Now you can inspect the calibrated data again (see Figure 23). Except for random scatter, things look pretty good.

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='day2_TDEM0003_10s_norx',
       xaxis='time',yaxis='amp',ydatacolumn='corrected',
       field='2',spw='0:4~60',correlation='RR,LL',
       avgchannel='64',coloraxis='antenna2')
</source>

You can use the Mark and Locate buttons to assess that the remaining scatter seems random, i.e., no particular antenna or time range appears to be responsible (see Figure 24).

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='day2_TDEM0003_10s_norx',
       xaxis='uvdist',yaxis='amp',ydatacolumn='corrected',
       field='3',spw='0~1:4~60',correlation='RR,LL',
       avgchannel='64',coloraxis='spw')
</source>

==Split==

Now we split the data into individual files. This is not strictly necessary, as we can select the appropriate fields in later clean stages, but it is safer in case for example we get confused with later processing and want to fall back to this point (this is especially a good idea if we plan to do continuum subtraction or self-calibration later on). It also makes smaller individual files in case we want to copy to another machine or want to share the data with a colleague.

Here, we split off the target (field 3). Note that to avoid a bug introduced in CASA 5.4.0 causing excessive flagging of the data (due to the combination of uvcontsub and statwt further on), we will not include the cross-hand polarization products (RL and LR) in '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.manipulation.split.html#casatasks.manipulation.split split]'''. As we will only make Stokes I maps of the sources, RL and LR are not used anyway in the below.  

<pre style="background-color: #d8bfd8;">
# Interactive CASA
default split
inp
vis='day2_TDEM0003_10s_norx'
outputvis='IRC10216'
field='3'
correlation='RR,LL'
inp
go
</pre>

<source lang="python">
# Pseudo-interactive CASA
split(vis='day2_TDEM0003_10s_norx',outputvis='IRC10216',field='3',correlation='RR,LL')
</source>

<!--

<span style="color:red">

</span>

To reinitialize the scratch columns for use by later tasks, we need to run ''{{clearcal}}'' for both new datasets.

<source lang="python">
# In CASA
clearcal(vis='J0954')
#
clearcal(vis='IRC10216')
</source>
-->
This concludes the calibration phase of the data reductions.

== UV Continuum Subtraction and Setting Up for Self-Calibration==

[[Image:CASA5.7.0-F25-irc10216-uvspec.jpg|200px|thumb|Figure 25: UV-plot of the spectral line signal in both spw for IRC+10216.]]

Now we can make a vector averaged uv-plot of the calibrated target spectral line data. It is important to note that we will only see signal in such a plot if (1) the data are well calibrated, and (2) there is significant signal near the phase center of the observations or if the line emission (or absorption) is weak but extended. If this isn't true for your data, you won't be able to see the line signal in such a plot and will need to make an initial (dirty or lightly cleaned) line+continuum cube to determine the line-free channels. Generally, this is the recommended course for finding the line-free channels more precisely than is being done here, as weak line signal would not be obvious in this plot. 

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='IRC10216',field='',ydatacolumn='data',
       xaxis='channel',yaxis='amp',correlation='RR',
       avgtime='1e8',avgscan=True,spw='0~1:4~60',antenna='',coloraxis='spw')
</source>

In the Display tab, change the Unflagged Points Symbol to Custom and Style to circle or diamond and 3 pixels. 

The horned profile typical of an expanding shell should be seen, as in Figure 25. From this plot, we can guess that strong line emission is restricted to channels 18 to 47 (zoom in if necessary to see exactly what the channel numbers are). 

In the Data tab, under Averaging, we can also click on All Baselines to average all baselines, but this is a little harder to see.

Now we want to use the line-free channels to create a model of the continuum emission that can be subtracted to create a line-only data set. We want to refrain from going too close to the edges of the band. These channels are typically noisy, and we don't want to get too close to the line channels because we could only see strong line emission in the vector averaged uv-plot and there may be weaker line emission in the channels next to the channel interval with strong lines.

We can now either specify the line free channels via parameters ''fitspw='0~1:4~14;51~59' ''and ''excludechans=False ''in '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.manipulation.uvcontsub.html#casatasks.manipulation.uvcontsub uvcontsub]''', or do the inverse selection via:

<pre style="background-color: #d8bfd8;">
# Interactive CASA
default uvcontsub
inp
vis='IRC10216'
fitspw='0~1:15~50'
excludechans=True
want_cont=True
inp
go
</pre>

<source lang="python">
# Pseudo-interactive CASA
uvcontsub(vis='IRC10216',fitspw='0~1:15~50',excludechans=True,want_cont=True)
</source>

which excludes the line channels from the continuum fit. (Note, it is ok to ignore the warning about "Forcing use of OLD VisibilityIterator".)

The parameter ''want_cont=True'' will produce two new data sets; IRC10216.contsub is the continuum subtracted line data, and IRC10216.cont is the continuum estimate.  Note, however, that it is still a multi-channel data set.

==Velocity Systems and Doppler corrections==

'''Note that for spectral line observations where the Doppler shift of the line is greater than the width of a channel, we recommend using the cvel2 task.  This step will correct for the Doppler shift during the observation.  For this particular tutorial, the Doppler shift of the line is not wider than a channel, therefore we may skip this step.'''

The VLA WIDAR correlator does not support Doppler ''tracking'', rather, Doppler ''setting'' is possible which will use the calculated sky frequency based on the velocity of a source at the start of the observation. The sky frequency is then fixed throughout that observing run. Typically, a fixed frequency is better for the calibration of interferometric data. The downside, however, is that a spectral line may shift over one or more channels during an observation. Task '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.imaging.tclean.html#casatasks.imaging.tclean tclean]''', takes care of such a shift when re-gridding the visibilities in velocity space (default is LSRK) to form an image. Sometimes, in particular when adding together different observing runs, it may be advisable to re-grid all data sets to the same velocity grid, combine all data to a single file, then Fourier transform and deconvolve. The tasks '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.manipulation.cvel2.html#casatasks.manipulation.cvel2 cvel2]''', '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.manipulation.concat.html#casatasks.manipulation.concat concat]''', and '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.imaging.tclean.html#casatasks.imaging.tclean tclean]''' serve this purpose respectively. The following run of '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.manipulation.cvel2.html#casatasks.manipulation.cvel2 cvel2]''' shows an example on how the parameters of '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.manipulation.cvel2.html#casatasks.manipulation.cvel2 cvel2]''' may be set. 

The '''IRC10216.contsub''' visibility spans the following channel range (see also the listobs output in the first part of the tutorial): 

<pre style="background-color: #d8bfd8;">
# Interactive CASA
default vishead
inp
vis='IRC10216.contsub'
mode='summary'
inp
go
</pre>

<source lang="python">
# Pseudo-interactive CASA
vishead(vis='IRC10216.contsub', mode='summary')
</source>

<pre style="background-color: #fffacd;">
SpwID  Name      #Chans   Frame   Ch0(MHz)  ChanWid(kHz)  TotBW(kHz) CtrFreq(MHz)  Corrs          
0      Subband:0     64   TOPO   36387.229       125.000      8000.0  36391.1670   RR  LL
1      Subband:0     64   TOPO   36304.542       125.000      8000.0  36308.4795   RR  LL
</pre>

For spw 0, this corresponds to about 1 km/s channel width. If we want to image the HC3N spectral line with a rest frequency of 36.39232 GHz over a velocity range of &minus;50km/s to 0km/s and a channel width of 5 km/s, we may decide to re-grid the visibilities in '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.manipulation.cvel2.html#casatasks.manipulation.cvel2 cvel2]''' as:

<pre style="background-color: #d8bfd8;">
# Interactive CASA
default cvel2
inp
vis='IRC10216.contsub'
outputvis='IRC10216.contsub-cveled'
mode='velocity'
nchan=10
start='-50km/s'
width='5km/s'
restfreq='36.39232GHz'
outframe='LSRK'
veltype='optical'
inp
go
</pre>

<source lang="python">
# Pseudo-interactive CASA
cvel2(vis='IRC10216.contsub', outputvis='IRC10216.contsub-cveled', 
     mode='velocity',nchan=10, start='-50km/s', width='5km/s',
     restfreq='36.39232GHz',outframe='LSRK', veltype='optical') 
</source>

This will create a new data set where the data is binned into the new grid. Since all data in measurement sets are stored in frequency space, an inspection with '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.information.vishead.html#casatasks.information.vishead vishead]''' now gives:

<source lang="python">
# Pseudo-interactive CASA
vishead(vis='IRC10216.contsub-cveled', mode='summary')
</source>

<pre style="background-color: #fffacd;">
SpwID  Name      #Chans   Frame   Ch0(MHz)  ChanWid(kHz)  TotBW(kHz) CtrFreq(MHz)  Corrs          
0      Subband:0     10   LSRK   36392.927       606.974      6070.6  36395.6588   RR  LL
</pre>

After the '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.manipulation.cvel2.html#casatasks.manipulation.cvel2 cvel2]''' step, the data can then be combined with other observations via '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.manipulation.concat.html#casatasks.manipulation.concat concat]''' and imaged in '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.imaging.tclean.html#casatasks.imaging.tclean tclean]''' with parameter ''specmode='cubedata' '' to conserve that velocity system and grid.

Note that '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.manipulation.cvel2.html#casatasks.manipulation.cvel2 cvel2]''' can also Hanning smooth the data, if needed (as an alternative to the stand-alone '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.manipulation.hanningsmooth.html#casatasks.manipulation.hanningsmooth hanningsmooth]''' task).

==Image the Spectral Line Data==

The step above was for information purposes, for the tutorial we now go back to the split off data of the target source IRC+10216 with 64 channels of which channel index 18 (the 19th channel) to channel index 47 (the 48th channel) clearly have line emission (seen in Figure 25). To be conservative, we assume a few channels extra on either side and we will use spw='15~50' for the channel index range with line emission.

Before any imaging, it is good practice to remove the effects of relative noise scatter introduced by flagging uneven bits in the visibility data between the channels and times. Of course line emission affects this so the range of channels with line emission should be excluded in the noise determination. The task to weigh the data using noise statistics is '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.manipulation.statwt.html#casatasks.manipulation.statwt statwt]'''. Note that this split-off data set does not have a corrected data column anymore so we are using '' datacolumn = 'data' ''.

<pre style="background-color: #d8bfd8;">
# Interactive CASA
default statwt
inp
vis='IRC10216.contsub'
fitspw='0~1:15~50'
excludechans=True
datacolumn='data'
inp
go
</pre>

<source lang="python">
# Pseudo-interactive CASA
statwt(vis='IRC10216.contsub',fitspw='0~1:15~50',excludechans=True,datacolumn='data')
</source>

First we make images from the continuum-subtracted, calibrated spectral line data. Because the spectral line emission from IRC+10216 has significant extended emission, it is very important to run '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.imaging.tclean.html#casatasks.imaging.tclean tclean]''' interactively, and make a clean mask. To make the cube a bit smaller and stay away from noisy edge channels we restrict the channel range using the spw parameter. Also, because we may want to self-calibrate the data, we have to save the model column (which is not the default).

'''Note that interrupting [https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.imaging.tclean.html#casatasks.imaging.tclean tclean] by Ctrl+C or by closing the viewer window may corrupt your visibilities or your CASA session.  We are currently implementing a command that will nicely exit to prevent this from happening, but for the moment try to avoid Ctrl+C.'''

<pre style="background-color: #98FB98;">Unfortunately, at this time there is a bug in CASA tclean with on-the-fly velocity conversion that renders the first two channels in the output image cube useless. Please ignore the first few channels of the output image.
</pre>

If '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.manipulation.cvel2.html#casatasks.manipulation.cvel2 cvel2]''' task was used in a previous step, then set the parameter specmode='cubedata', otherwise leave specmode='cube' within the '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.imaging.tclean.html#casatasks.imaging.tclean tclean]''' task.

<pre style="background-color: #d8bfd8;">
# Interactive CASA
default tclean
inp
vis='IRC10216.contsub'
spw='0:7~58'
imagename='IRC10216_HC3N.cube'
imsize=300
cell=['0.4arcsec']
specmode='cube'
outframe='LSRK'
restfreq='36.39232GHz'
perchanweightdensity=True
pblimit=-0.0001
weighting='briggs'
robust=0
niter=100000
threshold='3.0mJy'
interactive=True
savemodel='modelcolumn'
inp
go
</pre>

<source lang="python">
# Pseudo-interactive CASA
tclean(vis='IRC10216.contsub',spw='0:7~58',imagename='IRC10216_HC3N.cube',
       imsize=300,cell=['0.4arcsec'],specmode='cube',outframe='LSRK',restfreq='36.39232GHz',
       perchanweightdensity=True,pblimit=-0.0001,weighting='briggs',robust=0,
       niter=100000,threshold='3.0mJy',interactive=True,savemodel='modelcolumn')
</source>
<!-- old version used clean
clean(vis='IRC10216.contsub',imagename='IRC10216_HC3N.cube_r0.5',
      imagermode='csclean',
      imsize=300,cell=['0.4arcsec'],spw='0:7~58',
      mode='velocity',interpolation='linear',
      restfreq='36.39232GHz',outframe='LSRK',
      weighting='briggs',robust=0.5,
      interactive=True,
      threshold='3.0mJy',niter=100000)
-->

<!-- Lorant for reference, was in the old clean version; tclean uses gridder='standard'
* ''imagermode = csclean ''will invoke the Cotton-Schwab cleaning algorithm and the data will be regridded into a new output velocity frame, correcting for Doppler shifts of the line during the run (VLA data for each observation is always at a fixed sky frequency). The iterations are chosen as a high value to allow many clean cycles when needed. Typically, however, the threshold will kick in earlier and stop the cleaning process.  
-->
(The [https://casaguides.nrao.edu/index.php?title=VLA_CASA_Imaging VLA Imaging topical guide] has more information on the settings of the individual parameters of tclean. 

It will take a little while to grid the data, but the '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.visualization.viewer.html#casatasks.visualization.viewer viewer]''' will open when it's ready to start an interactive clean. Use the Channels tape deck buttons at the bottom (sometimes it will load on the side) of the Viewer display GUI to step through to the channel with the most extended (in angular size) emission (like channel index 22, 23 or 24 out of the 52 in the cube). 
[[Image:CASA5.7.0-F26-HC3N-cube.jpg|200px|thumb|Figure 26: Channel 30 shown (index 23 out of the 52 in the cube; the first few channels are noise only or unusable) for the HC3N cube in the interactive viewer, after adjusting the colors by entering a data range of -0.03 to 0.06 using the wrench icon, with the (green) polygon showing the mask drawn with the polygon tool. Double click inside it to finalize the clean area (then it will turn into a white line).]]

'''Important:''' In the green menu on top, '''select all channels and all polarizations for the clean mask''' and select the polygon tool [[File:Polygon_btn.png]]. The lines drawn by the polygon tool will be green, making it difficult to see on a green image.  The color map can be changed under the Data Display Options (click on the wrench icon or in the top menu, select Data, then select Adjust Data Display).  In the image, make a single mask that applies to all channels by clicking on the image to drop anchor points (see Figure 26) and drawing lines between the points; double click for the last point after which the anchor points can be adjusted (be careful not to click outside as that removes the polygon). Once the polygon region is satisfactory, '''double click inside it''' to save the mask region. The polygon will then turn white. Note, that if we had the time and patience, we could make a clean mask for each channel, and this would create a slightly better result. 

After making the mask, check that the emission in all the other channels fits within the mask you made using the tape deck controls to move back and forth. If you are unhappy with your mask, or perhaps you forgot to select all channels and all polarizations, you may remove the mask by selecting the Erase radio button located in the left-hand side of bright green box, then select the Rectangle drawing tool, draw a box over the mask, and then double click inside the box. To create a new mask, select the Add radio button (located just above the Erase radio button), then select the Polygon drawing tool to re-draw the mask. If you only need to remove part of the mask, you can add to the existing mask by drawing new polygons. Feel free to experiment with these tools.

'''Note''': If an interactive clean is started, and a mask is not created, then tclean will stop when you tell it to go on because it has nothing to clean. There is no default mask for interactive clean. If a non-interactive clean is started, and a mask is not created, then tclean will use the entire image as the mask.

To continue with '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.imaging.tclean.html#casatasks.imaging.tclean tclean]''' use one of the Next action buttons in the green area on the Viewer Display GUI: The red X [[File:clean-stop.png]] will stop '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.imaging.tclean.html#casatasks.imaging.tclean tclean]''' where it is; the blue arrow [[File:clean-continue.png]] will stop the interactive part of '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.imaging.tclean.html#casatasks.imaging.tclean tclean]''', but continue to clean non-interactively until reaching the stopping niter ('''iterations x cycles''') or threshold, whichever comes first. The green arrow [[File:clean-redo.png]] will clean until it reaches the iterations parameter on the left side of the green area. When the interactive viewer comes back use the tape deck to recheck that the mask encompasses what you think is real emission. If your mouse has a middle button, then by default it controls the image stretch.

For this example, the threshold has been set to ''threshold = '3mJy' ''to protect you from cleaning too deeply. With a careful clean mask you can clean too close to the thermal noise limit (the actual observed rms noise limit and not the theoretical one calculated for the proposal, as flagging, weather, etc., can affect what you actually get). It is ALWAYS best to clean each channel in a cube to a specific threshold than to stop by simply using the ''niter'' parameter, which can leave each channel cleaned to different levels. There are many ways to determine a suitable threshold. One way is to make a dirty image (''niter = 0''), open the cube using the viewer, go to a line free channel, select the box region tool, make a box near the field center about the size of your source, and double click inside. The rms noise of that channel will appear in the Statistics tab. Try a few different boxes and average the results; this is a good estimate of the rms per channel assuming your data are not dynamic range limited (i.e., noise can be higher in channels with strong signal). This is the absolute minimum for threshold and with no mask you probably shouldn't clean deeper than 3x this rms. 

Another method of getting these statistics would be to launch the version of the viewer that is run from the OS command line rather than the CASA command line. You can open this, and run other shell commands, from inside CASA using '!':

<source lang="python">
# In CASA
!casaviewer '<Name of dirty cube>' &
</source>

Like before, go to a line free channel, select the box region tool, make a box near the field center about the size of your source, and double click inside. Using this command line version of the viewer, the rms noise of that channel will appear in the terminal window from which the viewer was launched.

Keep cleaning by using the green Next Action arrow [[File:clean-redo.png]] until the residual displayed in the viewer looks noise-like. To speed things up, you might change the max_cycleniter parameter in the viewer to something like 300 or larger. This parameter can also be set in the task command. You will notice that in this particular case there are residuals that cannot be cleaned. These residuals are due to the extended, resolved out structure on size scales larger than the array is sensitive to (Largest Angular Scale (LAS) that the array is sensitive to can be calculated from the shortest baseline length), and potential residual phase and amplitude calibration errors. We will explore this in a few sections with self-calibration. 

Note that if you are unhappy with the cube and want to restart the cleaning, remove the original cube and all references to it in the memory by using '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.manipulation.rmtables.html#casatasks.manipulation.rmtables rmtables]'''.

<pre style="background-color: #98FB98;">
**Only use rmtables if you are unhappy with the data cube results.**
(Note, rm -rf will not properly clear the memory of the tables.)

# Pseudo-interactive CASA
rmtables(tablenames='IRC10216_HC3N.cube.*')
</pre>

Repeat the process for the SiS line using the call below. Note that the emission for this line is in spectral window index 1 and is less extended than the HC3N. The latter has to do with the different excitation situations of the two different lines. The SiS is excited closer to the central star than the HC3N. You may also notice that cleaning this cube takes much more iterations than the HC3N cube so you may want to increase the iterations in the viewer (at the left in the green area at the top) or just stop before it ends by itself. '''Remember to select All Channels and All Polarizations.'''

<source lang="python">
# Pseudo-interactive CASA
tclean(vis='IRC10216.contsub',spw='1:7~58',imagename='IRC10216_SiS.cube',
       imsize=300,cell=['0.4arcsec'],specmode='cube',outframe='LSRK',restfreq='36.30963GHz',
       perchanweightdensity=True,pblimit=-0.0001,weighting='briggs',robust=0,
       niter=100000,threshold='3.0mJy',interactive=True,savemodel='modelcolumn')
</source>
<!-- old version used clean
clean(vis='IRC10216.contsub',imagename='IRC10216_SiS.cube_r0.5',
      imagermode='csclean',
      imsize=300,cell=['0.4arcsec'],spw='1:7~58',
      mode='velocity',interpolation='linear',
      restfreq='36.30963GHz',outframe='LSRK',
      weighting='briggs',robust=0.5,
      interactive=True, 
      threshold='3.0mJy',niter=100000)
-->

You can look at both cubes using the viewer, and use the tape deck controls to play the cube as a movie.
<source lang="python">
# Pseudo-interactive or Interactive CASA
viewer
</source>

==Image the Continuum data==

Below, the use of parameter ''mode='mfs' ''will make a single multi-frequency synthesis image out of the specified spw/channels. Again, make an interactive clean mask, for which (using interactive mode) the number of iterations has to be larger than zero. Since no threshold is set, you will need to stop cleaning when the residual looks noise-like using the red X Next Action button [[File:clean-stop.png]] (it will be done when the viewer comes back the second time). The continuum for IRC10216 is very weak but interesting, tracing the photosphere of the AGB star. 

The continuum data set produced with parameter ''wantcont=True'' in uvcontsub is the model fit. To image the continuum itself, use the line-free channels of the original data, not the model fit.

<source lang="python">
# Pseudo-interactive CASA
tclean(vis='IRC10216',spw='0~1:4~14,0~1:51~60',imagename='IRC10216.36GHzcont',
       imsize=300,cell=['0.4arcsec'],specmode='mfs',weighting='briggs',
       pblimit=-0.0001,robust=0,niter=100000,interactive=True)
</source>
<!-- old clean version
clean(vis='IRC10216',imagename='IRC10216.36GHzcont',
      mode='mfs',imagermode='csclean',
      imsize=300,cell=['0.4arcsec'],spw='0~1:5~14,0~1:48~59',
      weighting='briggs',robust=0.5,
      interactive=True)
-->
Now look at the result in the viewer:

<source lang="python">
# Pseudo-interactive or Interactive CASA
viewer
</source>

==Image Analysis and Viewing==

Next, we will make integrated intensity maps (moment 0) and intensity-weighted velocity maps (moment 1). For HC3N, we will also produce a velocity dispersion, peak flux, and median map. All are derived with '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.analysis.immoments.html#casatasks.analysis.immoments immoments]'''. In order to do this, we will want to know what channels the line emission starts and ends on and the rms noise in a single channel.

First, open with the viewer file IRC10216_HC3N.cube.image:

<source lang="python">
# Pseudo-interactive CASA 
viewer('IRC10216_HC3N.cube.image')
</source>

Use the tape deck controls to see which channels have significant line emission. For HC3N, the line channel range in the cube is 11 to 40, which is 18 to 47 in the original 64-channel uv-data. However, when imaging, the seven beginning and five ending channels were discarded by selecting 0:7~58, which resulted in a 52 channel image cube; this is also the same for SiS. 

[[Image:CASA5.7.0-F27-HC3N-mom0.jpg|200px|thumb|Figure 27: HC3N moment 0 map.]]
[[Image:CASA5.7.0-F28-SiS-mom0.jpg|200px|thumb|Figure 28: SiS moment 0 map.]]

Then use the tape deck to go to a line free channel, select the box region tool and make a box. When you double click in the box, the image statistics for the channel you are using will be calculated. If not already showing a Regions panel in the viewer, you can see the results by going to the menu and select View and then Regions. This opens up the Regions panel, which has a "Statistics" tab; click it to read the rms. Move the box around a bit to see what the variation in rms noise is. You should get something around 2 mJy/beam or below (typically 1.4 mJy/beam is seen). Note that the rms is much worse in channels with strong emission because of the low dynamic range of these data. To make the box tool to go away, i.e., if you want to make a new one, press the Escape key. 

Now let's make the moment 0 and moment 1 maps. For moment 0, it's best to limit the calculation to image channels with significant signal in them, but not to apply a flux cutoff, as this will bias the derived integrated intensities upward. For large cubes, when the beam changes as a function of frequency, [https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.analysis.immoments.html#casatasks.analysis.immoments immoments] will convolve all planes to a common beam. '''The warnings about the kernel sampling may be ignored.'''

<pre style="background-color: #d8bfd8;">
# Interactive CASA
default immoments
inp
imagename='IRC10216_HC3N.cube.image'
chans='11~40'
outfile='IRC10216_HC3N.cube.image.mom0'
inp
go

imagename='IRC10216_SiS.cube.image'
outfile='IRC10216_SiS.cube.image.mom0'
inp
go
</pre>

<source lang="python">
# Pseudo-interactive CASA
immoments(imagename='IRC10216_HC3N.cube.image',chans='11~40',
          outfile='IRC10216_HC3N.cube.image.mom0')

immoments(imagename='IRC10216_SiS.cube.image',chans='11~40',
          outfile='IRC10216_SiS.cube.image.mom0')
</source>

Now use the viewer (see Figures 27 & 28):

<source lang="python">
# Pseudo-interactive CASA
viewer('IRC10216_HC3N.cube.image.mom0')

viewer('IRC10216_SiS.cube.image.mom0')
</source>

For moment 1, it is essential to apply a conservative flux cutoff to limit the calculation to high signal-to-noise areas. Here we use about 5 sigma.

[[Image:VLT-HC3Nmom0-SiSmom0-36GHzcont-CASA5.5.0.png|200px|thumb|Figure 29: The VLT V-band image from Le&atilde;o et al. (2006, A&A, 455,187) showing dust rings, overlaid with white HC3N moment 0 contours, black SiS moment 0 contours, and red 36 GHz continuum contours.]]

<source lang="python">
# Pseudo-interactive CASA
immoments(imagename='IRC10216_HC3N.cube.image',moments=[1],
          chans='11~40',excludepix=[-100,0.01],
          outfile='IRC10216_HC3N.cube.image.mom1')

immoments(imagename='IRC10216_SiS.cube.image',moments=[1],
          chans='11~40',excludepix=[-100,0.01],
          outfile='IRC10216_SiS.cube.image.mom1')
</source>

Finally, we will do velocity dispersion, peak flux, and median map in a single step for HC3N. Task '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.analysis.immoments.html#casatasks.analysis.immoments immoments]''' can perform these steps even though the latter two are not moments in a mathematical sense. Check the help file to find out the options. Peak flux and median are produced with the moment parameter set to 8 and 3.

<source lang="python">
# Pseudo-interactive CASA
immoments(imagename='IRC10216_HC3N.cube.image',moments=[2,8,3],
          chans='11~40',excludepix=[-100,0.01],
          outfile='IRC10216_HC3N.cube.image.extramoms')
</source>

The above task will create the following files:
* IRC10216_HC3N.cube.image.extramoms.weighted_dispersion_coord
* IRC10216_HC3N.cube.image.extramoms.median
* IRC10216_HC3N.cube.image.extramoms.maximum

Now use the '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.visualization.viewer.html#casatasks.visualization.viewer viewer]''' to further explore the images.


For fun, you can download the VLT V-band image at http://casa.nrao.edu/Data/EVLA/IRC10216/irc_fors1_dec_header.fits kindly provided by Izan Le&atilde;o and overlay the moment images and 36 GHz continuum (see Figure 29). More information about the dust properties can be found in the Le&atilde;o et al. (2006) paper http://adsabs.harvard.edu/abs/2006A%26A...455..187L.

Note that it should be possible to create [https://casa.nrao.edu/docs/cookbook/casa_cookbook008.html#section%3Adisplay.image.pvdiagram position velocity cuts] using the viewer. Masking the data cube to extract the emission is described in the [[Masking_images_for_analysis|Masking Images for Analysis]] page.

== Spectrum Fitting ==

Frequently, we may want to fit Gaussians or polynomials to the spectral line in the data cube. This can be done with CASA's '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.analysis.specfit.html#casatasks.analysis.specfit specfit]''' task; '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.analysis.specfit.html#casatasks.analysis.specfit specfit]''' can fit those functions to an average spectrum defined by some area, or, alternatively, for each pixel. We will fit an average spectrum.

=== Fitting an Average Spectrum ===

[[Image:CASA5.7.0-F30-HC3N-spectral-region.jpg|200px|thumb|Figure 30: The viewer showing the cube at plane 23 (#29 of 64). The magenta box marks the region that is averaged for the spectral profile.]]

[[Image:CASA6.2.0-F31-HC3N-spectral-profile.jpg|200px|thumb|Figure 31: Spectral profile of the HC3N line within the magenta box.]]

First, we want to inspect the spectrum. Load the image into the '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.visualization.viewer.html#casatasks.visualization.viewer viewer]''' (select the HC3N image cube) and mark a region with the mouse button that is assigned to the rectangular '''R''' region in the tool bar. It is probably best to do this at a channel that shows the entire extent of the source (see Figure 30). When done, double click in the center of the box to see information about the image cube in the terminal and then select '''spectral profile''' from the '''Tools''' menu. The average spectrum will be displayed in a separate panel (see Figure 31).

To fit this profile in '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.analysis.specfit.html#casatasks.analysis.specfit specfit]''', we need a region file outlining the 2-D region that is averaged (the magenta box in the viewer screenshot).
In the following, we use the new CASA region format (CASA 3.3 and higher) that is described at the [[CASA Region Format | CASA Region Format]] page. Following the guidelines on that page, we create a file named  '''specfit.crtf''' that describes a box with its [[x1,y1],[x2,y2]] corners in J2000 RA DEC coordinates.  Use any text editor (e.g., emacs, vi, etc.) to create the file '''specfit.crtf''' with the CRTF identifier line as well as the following box parameters:

<pre style="background-color: #fffacd;">
#CRTFv0
box[[09:47:59.2, 13.16.24], [09:47:55.8, 13.17.09]]
</pre>

We will fit two Gaussians to the two peaks of the spectrum (the scientific merit is debatable).  A file with initial values for the fit can be provided via the parameter ''estimates'' (see ''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.analysis.specfit.html#casatasks.analysis.specfit specfit]''' for more information). Here we will let CASA determine the start values by itself:  

<pre style="background-color: #d8bfd8;">
# Interactive CASA
default specfit
inp
imagename='IRC10216_HC3N.cube.image'
region='specfit.crtf'
estimates=''
ngauss=2
inp
go
</pre>

<source lang="python">
# Pseudo-interactive CASA
myfit = specfit(imagename='IRC10216_HC3N.cube.image',region='specfit.crtf',estimates='',ngauss=2) 
</source>

Note that for the latter (pseudo-interactive CASA), the output is stored in a Python dictionary called "myfit", as well as printed to the CASA logger.  You should get something similar to the following output, the results will vary based on your cleaning regions, depths, etc.:

<pre style="background-color: #fffacd;">
Results for component 0:
    Type     : GAUSSIAN
    Peak     : 4.59 +/- 0.38 mJy/beam
    Center   : -16.60 +/- 0.41 km/s
               15.85 +/- 0.40 pixel
    FWHM     : 9.04 +/- 0.99 km/s
               8.77 +/- 0.96 pixel
    Integral : 44.2 +/- 6.1 mJy/beam.km/s
Results for component 1:
    Type     : GAUSSIAN
    Peak     : 3.82 +/- 0.32 mJy/beam
    Center   : -33.98 +/- 0.59 km/s
               32.72 +/- 0.57 pixel
    FWHM     : 13.4 +/- 1.5 km/s
               13.1 +/- 1.5 pixel
    Integral : 54.7 +/- 7.7 mJy/beam.km/s
</pre>
We can see this seems to have caught the two peaks reasonably well.

<!--
<pre style="background-color: #fffacd;">
Results for component 0:
    Type     : GAUSSIAN
    Peak     : 5.85 +/- 0.49 mJy/beam
    Center   : -16.22 +/- 0.38 km/s
    15.47 +/- 0.37 pixel
    FWHM     : 8.13 +/- 0.90 km/s
    7.89 +/- 0.87 pixel
    Integral : 50.7 +/- 7.1 mJy/beam.km/s
Results for component 1:
    Type     : GAUSSIAN
    Peak     : 4.64 +/- 0.38 mJy/beam
    Center   : -33.33 +/- 0.64 km/s
    32.09 +/- 0.62 pixel
    FWHM     : 14.9 +/- 1.7 km/s
    14.5 +/- 1.6 pixel
    Integral : 74 +/- 10 mJy/beam.km/s
</pre>

=== Spectral Fitting pixel by pixel ===

[[Image:CASA5.7.0-F32-irc-amp0.jpg|200px|thumb|Figure 32: Output image of '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.analysis.specfit.html#casatasks.analysis.specfit specfit]''' of the amplitude of the first Gaussian component.]]

Task '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.analysis.specfit.html#casatasks.analysis.specfit specfit]''' can also fit Gaussians to every spectrum in each single spatial pixel. The following command will do this within the '''spectfit.crtf''' region defined above:

<pre style="background-color: #d8bfd8;">
# Interactive CASA
default specfit
inp
imagename='IRC10216_HC3N.cube.image'
region='specfit.crtf'
ngauss=2
multifit=True
amp='fit-amp.image'
center='fit-center.image'
fwhm='fit-fwhm.image'
logresults=True
inp
go
</pre>

<source lang="python">
# Pseudo-interactive CASA
pixfit = specfit(imagename='IRC10216_HC3N.cube.image', region='specfit.crtf', ngauss=2,
         multifit=True,  amp='fit-amp.image', center='fit-center.image', fwhm='fit-fwhm.image')
</source>

In this example, '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.analysis.specfit.html#casatasks.analysis.specfit specfit]''' will produce three images per Gaussian, images that map the best fit values of the Gauss peaks (amplitudes), velocity centers, and full widths at half maximum. For the first Gaussian, the image to the right displays the amplitude image '''fit.amp.image_0'''.

-->

== Self-Calibration ==

The many different aspects of self-calibration could fill several CASA guides. Here we describe a simple process for this particular relatively low S/N data (low S/N per channel, at least).

While running '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.imaging.tclean.html#casatasks.imaging.tclean tclean]''' above, the model column for each channel will have been filled with the clean model (if you made a Fourier transform of this model, you would see an image of the clean components). 

We choose to do the self cal on the spw=1 SiS line data because it has the strongest emission in a single channel and is a bit more compact than the HC3N data. We will run '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.gaincal.html#casatasks.calibration.gaincal gaincal]''' specifying the channel in the uv-data that has the brightest peak in the image (use the '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.visualization.viewer.html#casatasks.visualization.viewer viewer]''' to figure out which channel this is for spw=1), and note the peak flux. Since we started the image with a channel range, we need to account for the fact that the image channel numbers do not map exactly to the uv-data channel numbers (they are off by seven, so that channel index 12 (or e.g., 39) in the image is channel index 19 (or 46) in the uv-data). 

Next we need to understand the S/N of the data. In particular, to self-cal, we need enough signal on a single baseline within the chosen solution interval, ''solint'',  to get a S/N of at least 3. We calculated above an average rms noise of about 2 mJy/beam/channel for the entire time range (about 95 minutes of on-source time) and all antennas (16). We can use our knowledge of the radiometer equation (see [http://go.nrao.edu/vla-rms VLA Sensitivity]) where rms scales as 1/sqrt(time * #baselines), and the number of baselines= N(N-1)/2 and N=# of antennas. So the rms noise on one baseline, for one 10 second integration in this observation is given by:

<math> {\rm RMS(baseline)} = {\rm 2\ mJy\ beam^{-1}\ channel^{-1}} \sqrt{ \frac{95\times 60\ {\rm sec}}{10\ {\rm sec}}\times\frac{16\times 15}{2\times 1}}\ {\sim}\ 500\ {\rm mJy\ beam^{-1}\ channel^{-1}}</math>

[[Image:CASA5.7.0-F32-timeonsource.jpg|200px|thumb|Figure 32: Plot to estimate the time on source.]]

The 95 minutes of on-source time can be estimated from a plot where we can sum up the amount of time on a source (see Figure 32):

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='day2_TDEM0003_10s_norx',field='3',ydatacolumn='corrected',
       xaxis='time',yaxis='amp',correlation='RR,LL',
       avgchannel='64',spw='1:4~60',antenna='')
</source>

This analysis suggests that the rms noise on one baseline, for one 10 second integration, is about 500 mJy. In contrast, the peak flux density in the strongest SiS channel is only about 200 mJy (this can be checked using the '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.visualization.viewer.html#casatasks.visualization.viewer viewer]'''). Since the emission is fairly compact, most baselines will see approximately this peak flux. This is why we choose the more compact of the two possible lines.  Thus, a 10 second solution interval is not enough to get a SNR of at least 3 on a 200 mJy peak.   We need to use a solint large enough so that the rms noise is not worse than about 1/3 of 200 mJy. Thus, a solint of 10 minutes is about the shortest we can use and be reasonably confident of the solutions. 

Now we run '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.gaincal.html#casatasks.calibration.gaincal gaincal]''' with the solint we have determined. Note that because our desired solint is more than the scan time, we need to include parameter ''combine='scan' ''.

<source lang="python">
# Pseudo-interactive CASA
gaincal(vis='IRC10216.contsub',caltable='pcal_ch19one_10min',
        spw='1:19',calmode='p',solint='10min',combine='scan',
        refant='ea02',minsnr=3.0)
</source>

[[Image:CASA5.7.0-F33-selfcal-PhasevsTime.jpg|200px|thumb|Figure 33: Phase-only self-calibration solutions with 10 minute solint (second page of the display). ]]
Next, look at the solutions:

<source lang="python">
# Pseudo-interactive CASA
plotms(vis='pcal_ch19one_10min',gridrows=3,gridcols=3,
         xaxis='time',yaxis='phase',iteraxis='antenna',
         coloraxis='corr',plotrange=[0,0,-50,50])
</source>

For some antennas we can see some trends away from zero: ea08 and ea27 are examples, and we can also see some smaller variations with time (see Figure 33).

Now we explore whether applying this solution actually improves matters. To do this, run '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.applycal.html#casatasks.calibration.applycal applycal]''' to apply the solutions to the line data set for both spw.  We will set the interpolation parameter to use linearPD, this scales the phases by the frequency offset for the different spw.  We only set linearPD for spw=0, since we derived the solutions for spw=1 in the above '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.gaincal.html#casatasks.calibration.gaincal gaincal]''' task. Use parameter ''spwmap'' to tell it that the solutions derived for spw=1 should be applied to both spw=0 and spw=1.  Again, it's important to set parameter ''calwt=False''.

<source lang="python">
# Pseudo-interactive CASA
applycal(vis='IRC10216.contsub',field='',spw='0,1',interp=['linearPD',''],
        gaintable=['pcal_ch19one_10min'],spwmap=[[1,1]],calwt=False)
</source>

'''Note:''' In this example, we ran the self-cal steps on the full uv continuum subtracted spectral line data set. For a more complex iterative self-calibration procedure, you may find it easier to split off the channel/spw you want to experiment on with '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.manipulation.split.html#casatasks.manipulation.split split]''', and then do all the imaging ('''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.imaging.tclean.html#casatasks.imaging.tclean tclean]''') and '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.gaincal.html#casatasks.calibration.gaincal gaincal]''' steps with it. The '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.gaincal.html#casatasks.calibration.gaincal gaincal]''' tables created on the single channel can still be applied with '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.applycal.html#casatasks.calibration.applycal applycal]''' to the multi-channel/spw data set. If you do this, keep in mind that once split, the single-channel data will have its spw ID reset to 0 (you can check with '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.information.listobs.html#casatasks.information.listobs listobs]'''), no matter what spw it came from. Thus, in order to '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.calibration.applycal.html#casatasks.calibration.applycal applycal]''' with it, you would need parameter ''spwmap=<nowiki>[[0,0]]</nowiki>''.

To save time, we will use the clean mask we made before and run in a non-interactive mode. We can use a mask over again as long as the number of channels in the '''[https://casadocs.readthedocs.io/en/v6.2.0/api/tt/casatasks.imaging.tclean.html#casatasks.imaging.tclean tclean]''' call haven't changed. We can change cell or imsize and it will still do the right thing.

<source lang="python">
# Pseudo-interactive CASA
tclean(vis='IRC10216.contsub',spw='0:7~58',imagename='IRC10216_HC3N.cube.pselfcal',
       imsize=300,cell=['0.4arcsec'],specmode='cube',outframe='LSRK',restfreq='36.39232GHz',
       weighting='briggs',pblimit=-0.0001,perchanweightdensity=True,robust=0,
       niter=100000,threshold='3.0mJy',mask='IRC10216_HC3N.cube.mask')
</source>

<source lang="python">
# Pseudo-interactive CASA
tclean(vis='IRC10216.contsub',spw='1:7~58',imagename='IRC10216_SiS.cube.pselfcal',
     imsize=300,cell=['0.4arcsec'],specmode='cube',outframe='LSRK',restfreq='36.30963GHz',
     weighting='briggs',pblimit=-0.0001,perchanweightdensity=True,robust=0,
     niter=100000,threshold='3.0mJy',mask='IRC10216_SiS.cube.mask')
</source>

Now investigate the original and self-calibrated images in the viewer.  We will find that even this single self-cal step significantly improves the images. Try opening both versions of the SiS image cubes. Select a bright channel via the tape deck controls, then use the wrench and pwrench tools to make a plot setting the same image range for both cubes, and two panels in x, then to see both images of that channel side-by-side click the blink toggle.

[[Image:CASA5.7.0-F34-SiS-comparison.jpg|600px|]]

Figure 34: Original and self-cal SiS images for channel index 13/52 (zoomed in once), notice the increase in peak flux.

Repeat for HC3N (see Figure 35):

[[Image:CASA5.7.0-F35-HC3N-comparison.jpg|600px|]]

Figure 35: Original and self-cal HC3N images for channel index 23/52 (zoomed in once).

Now we can redo the moment images with the improved cubes.  Be sure to change the output file names. 

If you want to use the self-calibration solutions of the line data on the continuum, use the following:

<source lang="python">
# Pseudo-interactive CASA
applycal(vis='IRC10216',field='',spw='0,1',gaintable=['pcal_ch19one_10min'],spwmap=[[1,1]],calwt=False)

tclean(vis='IRC10216',spw='0~1:4~14,0~1:51~60',imagename='IRC10216.36GHzcont.pselfcal',
       imsize=300,cell=['0.4arcsec'],specmode='mfs',weighting='briggs',niter=100000,
       pblimit=-0.0001,robust=0,threshold='0.5mJy',mask='IRC10216.36GHzcont.mask')
</source>
Note that using the self-calibration of the extended line does not necessarily make the continuum better; it may seem like the imperfection of using the extended line channel has imprinted on the continuum.


Questions about this tutorial? Please contact the [https://science.nrao.edu/observing/helpdesk NRAO Science Helpdesk].

[[Main Page | &#8629; '''CASAguides''']]

<!-- --[[User:Cbrogan|Crystal Brogan]]    
--additions:[[User:Jott| Juergen Ott]], [[User:Mkrauss| Miriam Krauss]], [[User:Dshepher| Deb Shepherd]]
--corrections: [[User:Lsjouwerman| Lorant Sjouwerman]] (4.4.0, 2015/06/24)
-- 4.5 update: [[User:Jott| Juergen Ott]] (2015/11/3)
--content update: [[User:Hmedlin| Heidi Medlin]] (4.5.2, 2016/02/26) 
--content update: [[User: Tperreau| Tony Perreault]] (4.5.2, 2016/03/01)
-- 4.6 update and corrections: [[User:Lsjouwerman| Lorant Sjouwerman]] (2016/05/05)
-- 5.0 update and corrections: [[User:Lsjouwerman| Lorant Sjouwerman]] (2017/09/20)
-- 5.0 added Interactive CASA method: [[User:Hmedlin| Heidi Medlin]] (2018/06/27)
-- 5.4 all new figs with new names, added plotms commands for plotcal (but not replaced plotcal), mentioned issues with uvcontsub/statwt, viewer (regions) and tclean first channel out (and mom0 image, not mentioned here), converted clean to tclean, used new flux scale, added rmtables and contimuum calibrated with line: [[User:Lsjouwerman| Lorant Sjouwerman]] (2018/11/30)
-- 5.5 update and corrections (removed plotcal task and images and updated many tclean parms): [[User:Hmedlin| Heidi Medlin]] (2019/06/25), minor text mods [[User:Lsjouwerman| Lorant Sjouwerman]] (2019/06/14)
-- 5.7 update new figures from plotms & tclean runs (Fig 12 - 35 inclusive) from bandpass calibration to end of tutorial; set parameter solnorm=False in bandpass task: [[User:Tperreau| Tony Perreault]] (2020/10/15)
-- 6.2 updated all links to CASA6.2.0 and two figures (plotweather & spectral profile) [[User:Hmedlin| Heidi Medlin]] (2021/06/23) 
-->

{{Checked 6.2.0}}