EVLA high frequency Spectral Line tutorial - IRC+10216 part1: Difference between revisions

From CASA Guides
Jump to navigationJump to search
No edit summary
(Undo revision 11261 by Jott (talk))
 
(22 intermediate revisions by 2 users not shown)
Line 1: Line 1:
[[Category:EVLA]][[Category:Calibration]][[Category:Spectral Line]]
[[Category:EVLA]][[Category:Calibration]][[Category:Spectral Line]]


This is the first part of the EVLA IRC+10216 data reduction tutorial. [http://casaguides.nrao.edu/index.php?title=EVLA_high_frequency_Spectral_Line_tutorial_-_IRC%2B10216_CASA3.4  CASA 3.4 version]


'''Note: this tutorial will use the current casapy-stable version of CASA. Most of the online documentation, however, will is for the 3.2.1 release (including the direct links of tasks to the CASA webpages). Please refer to the inline help of the tasks for the additional parameters. Or try to get a hand of the latest cookbook prepared for this event.'''


<!--
This is the first part of the EVLA IRC+10216 data reduction tutorial. [http://casaguides.nrao.edu/index.php?title=EVLA_high_frequency_Spectral_Line_tutorial_-_IRC%2B10216_part2  IRC+10216 Part 2: continuum subtraction, imaging, selfcal] 
'''
* added: opacity section
* added: cvel section (but needs final inputs)
* added; hanning smoothing?
* added: link to VLA gain curves
* continuum subtraction in the image plane?
* added: link to pV diagrams
* added: link to masking
* iterations in plotms
'''
-->


== Overview ==
== Overview ==
Line 21: Line 10:


This tutorial describes the data reduction for two spectral lines observed toward the AGB star IRC+10216.  
This tutorial describes the data reduction for two spectral lines observed toward the AGB star IRC+10216.  
In this EVLA OSRO1 mode observation one subband was observed in each of two basebands, the subbands were centered on the HC3N and SiS lines near 36 GHz. The raw data were loaded into CASA with {{importevla}}, where zero and shadowed data were flagged. Then the data were "{{split}}", so we could average from the native 1 second integrations to 10 seconds, select only antennas with Ka-band receivers, and select only spectral windows (called spw in CASA) with Ka-band data. This produces a significantly smaller dataset for processing.
In this EVLA "OSRO1" mode observation one subband was observed in each of two basebands, with the subbands centered on the HC3N and SiS lines near 36 GHz. The raw data were loaded into CASA with {{importevla}}, where zero and shadowed data were flagged (note that this is done by the [https://archive.nrao.edu/archive/advquery.jsp NRAO Data Archive]). Then the data were {{split}}, so we could average from the native 1-second integrations to 10 seconds, select only antennas with Ka-band receivers, and select only spectral windows (called spw in CASA) with Ka-band data. This produces a significantly smaller dataset for processing.


IRC+10216 (CW Leo) is the brightest star in the sky at 5 microns but only 16th magnitude visually.  It was discovered during the first survey
IRC+10216 (CW Leo) is the brightest star in the sky at 5 microns but only 16th magnitude visually.  It was discovered during the first survey of the infrared sky, carried out by Bob Leighton and Gerry Neugebauer in 1965.  An Asymptotic Giant Branch star, it 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 emission; the continuum emission seen at radio wavelengths probes the actual 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 material drops, the chemistry freezes.  But the molecules continue their long coast outward into the Galaxy, and as the shell thins ultraviolet light from the ambient galactic radiation field penetrates and initiates a new chemistry in the gas.
of the infrared sky, carried out by Bob Leighton and Gerry Neugebauer in 1965.  An Asymptotic Giant Branch star, it 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 emission; the continuum emission seen at radio wavelengths probes the actual 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 material drops, the chemistry freezes.  But the molecules continue their long coast outward into the Galaxy, and as the shell thins ultraviolet light from the ambient galactic radiation field penetrates and initiates a 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.  The result is the
SiS, a simple molecule created in the dense inner envelope chemistry is photodissociated as it coasts out into the shell.  The result is the centrally condensed emission we see in the present observations. Interestingly, in the next lower transition at 18 GHz, the line shape is much different from what we see here.  At the extreme velocities in the profile, very bright narrow emission is seen which has been interpreted as maser emission.  The interested student can find EVLA observations of this line in the archive.
centrally condensed emission we see in the present observations. Interestingly, in the next lower transition at 18 GHz, the line shape is
much different from what we see here.  At the extreme velocities in the profile, very bright narrow emission is seen which has been interpreted
as maser emission.  The interested student can find EVLA observations of this line in the archive.


HC3N, a much more complex species, is created by the photochemistry which becomes active as atoms and pieces of molecules destroyed by
HC3N, a much more complex species, is created by the photochemistry which becomes active as atoms and pieces of molecules destroyed by ultraviolet radiation undergo the next phase of chemistry in the shell. HC3N has many vibrational modes which may be excited in addition to its rotational modes.  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, but during its brief existence its rise to abundance in the envelope results in a ring of emission, which is what is observed in this image made with the EVLA. A recent model by Cordiner & Millar (2009, ApJ, 697, 68) describes a new chemical model for the shell, which also takes into account the variation of mass loss by the star.  They show that in addition to purely chemical effects, local gas and dust density peaks play a role in shaping the observed emission.
ultraviolet radiation undergo the next phase of chemistry in the shell. HC3N has many vibrational modes which may be excited in addition to
its rotational modes.  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, but during its brief existence its rise to abundance in the envelope results in a ring of emission, which is what is observed in this image made with the EVLA. A recent model by Cordiner & Millar (2009, ApJ, 697, 68) describes a
new chemical model for the shell, which also takes into account the variation of mass loss by the star.  They show that in addition to
purely chemical effects, local gas and dust density peaks play a role in shaping the observed emission.


==Getting the data==
==Getting the data==


The post-split averaged data can be downloaded from http://casa.nrao.edu/Data/EVLA/IRC10216/day2_TDEM0003_10s_norx.tar
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] (data size: 1.1GB)


Once the download is complete, unpack the file:
Once the download is complete, unzip and unpack the file:


<source lang='bash'>
<source lang='bash'>
# in a terminal, outside of CASA:
# in a terminal, outside of CASA:
tar -xvf day2_TDEM0003_10s_norx.tar
tar -xzvf day2_TDEM0003_10s_norx.tar.gz
</source>
</source>


Line 61: Line 41:


For all EVLA observations, the operators keep an observing log. You can look at
For all EVLA observations, the operators keep an observing log. You can look at
the observing logs at the observing log [[http://www.vla.nrao.edu/cgi-bin/oplogs.cgi website]]. Pertinent information from this observation is repeated below.
the observing logs at the observing log [http://www.vla.nrao.edu/cgi-bin/oplogs.cgi website] by entering a date range covering only 26 Apr 2010, then selecting "TDEM0003". Pertinent information from this observation is repeated below:
 
<pre style="background-color: #E0FFFF;">
<pre style="background-color: #E0FFFF;">
INFORMATION FROM OBSERVING LOG:
INFORMATION FROM OBSERVING LOG:
Date of the observation: 26-April-2010
Date of the observation: 26-April-2010
There are no Ka-band receivers on ea11, ea13, ea14, ea16, ea17, ea18, ea26   
There are no Ka-band receivers on ea11, ea13, ea14, ea16, ea17, ea18, ea26   
Antennas ea10, ea06 are out of the array
Antenna ea06 is out of the array
Antenna ea12 is newly back
Antenna 15 has a focus timeout
The pointing is often bad on ea15
Antennas ea10, ea12, ea22 do not have good baseline positions
Antennas ea10, ea12, ea22 do not have good baseline positions
</pre>
</pre>


As mentioned in the log, antennas ea10, ea12, and ea22 do not have good baseline positions.
As mentioned in the log, 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
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.
any improved baseline positions that were derived after the observations were taken.


In CASA, we need to insert these corrections by hand using '''{{gencal}}'''. The resulting table will need to be supplied as a "prior" calibration to all subsequent calibration steps. The corrections can be ascertained from the [http://www.vla.nrao.edu/astro/archive/baselines/ EVLA/VLA Baseline Corrections page]. Be sure to carefully read to the bottom of that  
In CASA, we need to insert these corrections by hand using {{gencal}}. The resulting table will need to be supplied as a "prior" calibration to all subsequent calibration steps. The corrections can be ascertained from the [http://www.vla.nrao.edu/astro/archive/baselines/ EVLA/VLA Baseline Corrections page]. Be sure to carefully read to the bottom of that  
webpage to see how to calculate the additive corrections. The current case is simple as there is only a single correction for each antenna. In the future we will implement an automated lookup of the corrections.  
webpage to see how to calculate the additive corrections. The current case is simple as there is only a single correction for each antenna. In the future, CASA will include an automated lookup of the corrections.  


<source lang="python">
<source lang="python">
Line 88: Line 68:


'''Please note: '''if you are reducing VLA data taken before March 1, 2010, you need to set caltype='antposvla'.
'''Please note: '''if you are reducing VLA data taken before March 1, 2010, you need to set caltype='antposvla'.


There are two other types of "priors" that we will need to apply at each calibration stage described below:
There are two other types of "priors" that we will need to apply at each calibration stage described below:


# Gaincurve -- the gaincurve describes how each antenna behaves as a function of elevation, for each receiver band. Currently only gaincurves for the VLA/EVLA are available (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.
# Opacity correction -- the opacity of the observation can be computed from a seasonal model and/or weather station information. We are planning to have a task available for this information. At the moment, the [[CASA_EVLA_Scripts]] page hosts a script and a contributed task to display the weather information and to calculate the zenith opacities for each spectral window. 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 {{gaincal}}, {{applycal}}, {{bandpass}} etc.


(1) Gaincurve -- the gaincurve describes how each antenna behaves as a function of elevation, for each receiver band. Currently only gaincurves for the VLA/EVLA are available (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.
Now we want to plot the opacity of the atmosphere at 36 GHz at the time this observation was taken so it can be corrected for in
subsequent calibrations. There is not currently a task in CASA to do this however, an external task, plotWX.py, has been written that will do this. In the following we will download the contributed task, and bring it into CASA. We will then use the results to determine the opacity of the atmosphere at the time these observations were taken.


Download these two files:


(2) Opacity correction -- the opacity of the observation can be computed from a seasonal model and/or weather station information. We are planning to have a task available for this information. At the moment, the [[CASA_EVLA_Scripts]] page hosts a script to show the weather information and to calculate the zenith opacities for each spectral window. 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 {{gaincal}}, {{applycal}}, {{bandpass}} etc.  
[[File:Task_plotWX.py]] and [[File:PlotWX.xml]]. The first file contains the actual python code and the second the interface and inline help for CASA. Download the files and note that for obscure reasons the wiki will always capitalize the uploaded files abut we need to change this. In a regular terminal:


Download the script: [[File:Script_plotWX.py]]  and
<source lang='bash'>
# in a terminal, outside of CASA:
wget http://casaguides.nrao.edu/images/e/e4/PlotWX.xml
mv PlotWX.xml plotWX.xml
wget http://casaguides.nrao.edu/images/2/24/Task_plotWX.py
mv Task_plotWX.py task_plotWX.py
</source>


<pre style="background-color: #98FB98;">
Now, we go back to CASA and build the task itself with '''buildmytask''' (it is a shell command and not a CASA command, but we will run it from within CASA; it will create a few *.pyc and *cli files as well as mytasks.py). And then add it into your CASA session by running '''mytasks.py''':
replace this line at the top of the script


myMS='MSname.ms'
<source lang="python">
# In CASA
!buildmytasks
execfile('mytasks.py')
</source>


with
Now your brand-new '''plotWX''' task is available:


myMS='day2_TDEM0003_10s_norx'
<pre style="background-color: #fffacd;">
CASA <7>: inp plotWX
--------> inp(plotWX)
#  plotWX :: Plot elements of the weather table for a given MS
vis                =        ''       #  MS name
seasonal_weight    =        0.5        #  weight of the seasonal model
doPlot              =      True        #  set this to True to create a plot
async              =      False        #  If true the taskname must be started using plotWX(...)
</pre>


</pre>
'''seasonal_weight''' is a parameter that combines global seasonal weather conditions (measured over many years) with the weather information from the weather station obtained during the observations. 0.5 will combine both with equal weight.


and run the script.  
[[Image:day2_TDEM0003_10s_norx.plotWX.png|200px|thumb|right|plotWX weather table figure]]


Back in CASA, run the script via:
We will be running plotWX in a way that will assign the opacity list (one entry for each spw in ascending order) to the variable myTau:


<source lang="python">
<source lang="python">
# In CASA
# In CASA
execfile('Script_plotWX.py')
myTau = plotWX(vis='day2_TDEM0003_10s_norx', doPlot=T)
</source>
</source>


This will return two outputs:
With this input, the task
 
A file '''"day2_TDEM0003_10s_norx.plotWX.png"'''  with the elevation of the sun, the wind speed and direction, the temperature, precipitable water vapor (PWV) as functions of time over the observation, as well as the zenith opacities as a function of frequency (based on the average weather during the observations as well as long-term seasonal data).


 
* returns the opacity values to the logger:
[[Image:day2_TDEM0003_10s_norx.plotWX.png|200px|thumb|right|plotWX weather table figure]]
 
In addition, the task returns the opacity values to the logger:


<pre style="background-color: #fffacd;">
<pre style="background-color: #fffacd;">
Line 137: Line 131:
</pre>  
</pre>  


and it also returns a python variable with the opacities listed in the sequence of ascending spectral window:
* creates a file '''"day2_TDEM0003_10s_norx.plotWX.png"'''  with the elevation of the sun, the wind speed and direction, the temperature, 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 assigns the myTau variable to the list of opacities per spectral window:  


<source lang="python">
<source lang="python">
# In CASA
# In CASA
print myTau
myTau
[0.041019209411983566, 0.040779609355637569]
[0.041019209411983663, 0.040779609355637694]
</source>
</source>


The calibration task opacity parameters can now take either the variable directly '''opacity=myTau''' or the list of opacities like '''opacity=[0.0410, 0.0408]''' or, if a single value is good enough to correct for all spectral windows, one can simple use '''opacity=0.041'''.
The calibration task opacity parameters can now take either the variable directly '''opacity = myTau''' or the list of opacities like '''opacity = [0.0410, 0.0408]''' or, if a single value is good enough to correct for all spectral windows, one can simple use '''opacity = 0.041'''.
 


==Initial Inspection and Flagging==
==Initial Inspection and Flagging==
Line 206: Line 201:
</pre>
</pre>


In addition to source names, antenna names, etc. CASA indexes each of these parameters staring from 0. When, say a field is specified, one can use the index or the name, both is accepted.  
In addition to source names, antenna names, etc. CASA indexes each of these parameters starting from 0. To specify, e.g., a source field, one can use the index or the name -- both are accepted.  


 
Note that 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 EVLA Observing Tool 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.  
Note that 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 EVLA Observing Tool 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 reductions.  


<pre style="background-color: #E0FFFF;">
<pre style="background-color: #E0FFFF;">
Summary of Observing Strategy
Summary of Observing Strategy
Gain Calibrator: J0954+1743 field id=2
Gain Calibrator: J0954+1743, field id=2
Bandpass Calibrator: J1229+0203   field id=5
Bandpass Calibrator: J1229+0203, field id=5
Flux Calibrator: J1331+3030 (3C286) field id=7
Flux Calibrator: J1331+3030 (3C286), field id=7
Target: IRC+10216  field id=3
Target: IRC+10216, field id=3
Ka-band spws = 0,1
Ka-band spws = 0,1
</pre>  
</pre>  
Line 224: Line 218:
in case you want it later. This will be useful for picking a reference antenna --
in case you want it later. This will be useful for picking a reference antenna --
typically a good choice is an antenna close to the center of the array. Unless it  
typically a good choice is an antenna close to the center of the array. Unless it  
shows problems after inspection of the data, we provisionally chose ea02.
shows problems after inspection of the data, we provisionally choose ea02.


[[Image:elevationvstime.png|thumb|Elevation as a function of time (after selecting colorize by field).]]
[[Image:elevationvstime.png|thumb|Elevation as a function of time (after selecting colorize by field).]]
Line 232: Line 226:
</source>
</source>


Next, let's look at the elevation as a function of time for all sources. It's not the case for these data, but if the elevation is very low (usually at start or end of track) you may want to flag. Also, how near in elevation your flux calibrator is to your target will impact your ultimate absolute flux calibration accuracy. Unfortunately, the target and flux calibrator are not particularly well-matched for this observation, as you can show by plotting the elevation for each source (each sources has a different colors). We will be using data in spectral window 0, channels 4 to 60 for this plot as given by the '''spw='0:4~60'''' parameter, the general CASA selection syntax is described in the [http://www.aoc.nrao.edu/~sbhatnag/misc/msselection/msselection.html Measurement Selection Syntax Document]:
Next, let's look at the elevation as a function of time for all sources. It's not the case for these data, but if the elevation is very low (usually at start or end of track) you may want to flag. Also, how near in elevation your flux calibrator is to your target will impact your ultimate absolute flux calibration accuracy.  
 
Unfortunately, the target and flux calibrator are not particularly well-matched for this observation, as you can show by plotting the elevation for each source (each sources has a different colors). We will be using data in spectral window 0, channels 4 to 60 for this plot as given by the '''spw = '0:4~60'''' parameter; the general CASA selection syntax is described in the [http://www.aoc.nrao.edu/~sbhatnag/misc/msselection/msselection.html Measurement Selection Syntax Document]:


<source lang="python">
<source lang="python">
Line 246: Line 242:
[[Image:Zoom1_mark.png|thumb|Zooming in and marking region (hatched box)]]
[[Image:Zoom1_mark.png|thumb|Zooming in and marking region (hatched box)]]


Next, let's look at all the source amplitudes as a function of time.
Next, let's look at all the source amplitudes as a function of time:
 
<source lang="python">
<source lang="python">
# In CASA
# In CASA
Line 253: Line 250:
       avgchannel='64',spw='0:4~60', coloraxis='field')
       avgchannel='64',spw='0:4~60', coloraxis='field')
</source>
</source>


Now zoom in on the region very near zero amplitude for sources J0954+1743 and IRC+10216. To zoom, select the  
Now zoom in on the region very near zero amplitude for sources J0954+1743 and IRC+10216. To zoom, select the  
Zoom tool in lower left corner of  the {{plotms}} GUI, then you can left click to draw a box.  
Zoom tool [[File:casaplotms-zoom.png]] in lower left corner of  the {{plotms}} GUI, then you can left click to draw a box.  
Look for the low values (you may want to zoom a few times to really see the suspect points clearly). Now use the Mark Region and Locate buttons (located along the bottom of the GUI) to see which antenna is causing problems. The output is be shown in the logger. Since all the "located" baselines include ea12, this is the responsible antenna.
Look for the low values (you may want to zoom 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 (located along the bottom of the GUI) to see which antenna is causing problems. The output is be shown in the logger. Since all the "located" baselines include ea12, this is the responsible antenna.


<pre style="background-color: #98FB98;">
<pre style="background-color: #98FB98;">
Line 263: Line 259:


* When using the locate button it is important to have only selected a modest number  
* 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 thumnail),  
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  
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).  
(all of the information will be output to your terminal window, not the logger).  
Line 271: Line 267:


* After flagdata command flagging, you have to force a complete reload of the cache  
* 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 change anything
to look at the same plot again with the new flags applied. To do this, either check the "force reload" box in the lower left, or do Shift+Plot.
in the plotms GUI (the colorize parameter, antenna, spw, anything) and hit the
plot button.


* If the plotms tool does get hung during a plot try clicking the cancel button on the  
* If the plotms tool hangs during a plot, try clicking the cancel button on the  
load progress GUI, and/or if you see a "table locked" message try typing  
load progress GUI, and/or if you see a "table locked" message try typing  
clearstat on the CASA command line.
"clearstat" on the CASA command line.


* Occasionally plotms will get into a strange state that you cannot clear from inside.  
* Occasionally, plotms will get into a strange state that you cannot clear from inside.  
We are working on these issues, but for now, when all else fails, exit from casapy and  
We are working on these issues, but for now, when all else fails, exit from casapy and  
restart.   
restart.   
Line 285: Line 279:
</pre>  
</pre>  


Now click the clear region button, and 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.
Now click the clear region [[File:casaplotms-clear-region.png]] button, and 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
Check the other sideband by changing spw to 1:4~60. You will have to
rezoom. If you have trouble, click on the Mark icon and then back 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 next
zoom.  In spw 1, ea07 is bad from the beginning until after next
pointing run: 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'.
pointing run: 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 intial timerange, 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' for antenna, which removes ea07 from the plot (in CASA selection, ! deselects).  
If you set antenna to 'ea12' and zoom in on this initial timerange, 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' for antenna, which removes ea07 from the plot (in CASA selection, "!" deselects).  


We can set up a flagging command to get both bad antennas for the
We can set up a flagging command to get both bad antennas for the
Line 309: Line 303:


Note that because  the chosen timerange is limited to fields 2 and 3, the field parameter is not really
Note that because  the chosen timerange is limited to fields 2 and 3, the field parameter is not really
needed; however, flagdata will run fastest if you put as many constraints as possible.
needed; however, flagdata will run fastest if you put as many constraints as possible (since it will traverse less data).


Now remove the !ea07 from antenna and replot both spw, zooming in to
Now remove the !ea07 from antenna and replot both spw, zooming in to
Line 317: Line 311:
[[Image:IRC10216_uvdist1.png|thumb|Amplitude vs. uv-distance for IRC+10216, both spw]]
[[Image:IRC10216_uvdist1.png|thumb|Amplitude vs. uv-distance for IRC+10216, both spw]]


Lets look more closely at IRC+10216:
Let's look more closely at IRC+10216:


<source lang="python">
<source lang="python">
Line 326: Line 320:
</source>
</source>


You can see a
You can see a that there are some noisy high points. But now try:
that there are some noisy high points. But now try


<source lang="python">
<source lang="python">
Line 336: Line 329:
</source>
</source>


Most of the high points on IRC+10216 are due to large scale emission on short baselines, 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.
Now, we can see that most of the high points on IRC+10216 are due to large scale emission on short baselines, 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.


==Set Up the Model for the Flux Calibrator==
==Set Up the Model for the Flux Calibrator==
Line 343: Line 336:
the flux calibrator, possibly with a uvrange cutoff if necessary. More recently for the VLA/EVLA,  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 flux calibrator, possibly with a uvrange cutoff if necessary. More recently for the VLA/EVLA,  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).  


Currently, CASA contains models for the most common calibrators at the most common frequencies but not yet the new EVLA frequency bands S, Ku, and Ka. These will be added as son as they become available. One may check the availability of calibration models in {{setjy}}:
Currently, CASA contains models for the most common calibrators at the most common frequencies but not yet the new EVLA frequency bands S and Ka. These will be added as soon as they become available. One may check the availability of calibration models in {{setjy}}:
 


<source lang="python">
<source lang="python">
Line 351: Line 343:
</source>
</source>


The logger will now show the available models, e.g. 3C286_C.im, 3C48_K.im etc.
The terminal will now show the available models, e.g. 3C286_C.im, 3C48_K.im etc. ({{setjy}} will search in the working directory for images that may contain models, as well as in a CASA directory where known calibrator models are stored.)


As mentioned above, the Ka band does not yet have a model incorporated at this time as the full configuration cycle has not ended yet (the A-configuration data has not yet obtained as of Aug 2011). However, models of Ka band that are good for D, C, and B configurations are available from the [https://science.nrao.edu/facilities/evla/data-processing/flux-calibrator-models-for-new-evla-bands  EVLA calibrator model webpage for the new frequency bands].  
As mentioned above, the Ka band does not yet have a model incorporated at this time as the full configuration cycle has not ended yet. However, models for Ka-band that are good for D, C, and B configurations are available from the [https://science.nrao.edu/facilities/evla/data-processing/flux-calibrator-models-for-new-evla-bands  EVLA calibrator model webpage for the new frequency bands].  


From this page, please download the file 36.68 GHz CASA image of 3C286, this is the [http://www.aoc.nrao.edu/evlacal/3C286_36.68G.model.tgz direct link.] To do so, you may run  
From this page, please download the file 36.68 GHz CASA image of 3C286, this is the [http://www.aoc.nrao.edu/evlacal/3C286_36.68G.model.tgz direct link.] To do so, you may run  


<source lang='bash'>
<source lang='bash'>
Line 367: Line 358:


{{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.
{{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">
<source lang="python">
Line 377: Line 367:
<pre style="background-color: #fffacd;">
<pre style="background-color: #fffacd;">
The logger output for each spw is:
The logger output for each spw is:
J1331+3030  rawspwid=  0  [I=1.776, Q=0, U=0, V=0] Jy, (Perley-Butler 2010)
J1331+3030  rawspwid=  0  [I=1.7762, Q=0, U=0, V=0] Jy, (Perley-Butler 2010)
J1331+3030  rawspwid=  1  [I=1.779, Q=0, U=0, V=0] Jy, (Perley-Butler 2010)
J1331+3030  rawspwid=  1  [I=1.7794, Q=0, U=0, V=0] Jy, (Perley-Butler 2010)
</pre>
</pre>


The absolute fluxes for the frequencies have now been determined and one can proceed to the bandpass and complex gain calibrations.
The absolute fluxes for the frequencies have now been determined, and one can proceed to the bandpass and complex gain calibrations.


==Bandpass==
==Bandpass==
Line 398: Line 388:
</source>
</source>
[[Image:Nobandpass_phase.png|thumb|Phase as a function of channel for ea02 (after Custom and upping "Style" to 3.)]]
[[Image:Nobandpass_phase.png|thumb|Phase as a function of channel for ea02 (after Custom and upping "Style" to 3.)]]
The phase variation is modest ~10 degrees.  Now expand to all baselines that include ea02, then hit the "Plot" button.
 
If you find the points difficult to see, go to the "Display" tab, set "Unflagged Points Symbol" to "Custom", and increase the "Style" number, then replot.  The phase variation is modest ~10 degrees.  Now expand to all baselines that include ea02, and add an extra dimension of color:


<source lang="python">
<source lang="python">
Line 412: Line 403:
other antennas if you like.
other antennas if you like.


Now look at the phase as a function of time.
Now look at the phase as a function of ''time'':


<source lang="python">
<source lang="python">
Line 421: Line 412:
</source>
</source>


 
Expand to all antennas with ea02:
Expand to all antennas with ea02


<source lang="python">
<source lang="python">
Line 430: Line 420:
       avgchannel='64',spw='0:4~60',antenna='ea02', coloraxis='antenna2')
       avgchannel='64',spw='0:4~60',antenna='ea02', coloraxis='antenna2')
</source>
</source>
You may want to select "Custom" under "unflagged points symbol" and then change Style from "2" to "3" under the "Display" tab, "Unflagged Points Symbol" .


You can see that the phase variations are smooth, but do vary
You can see that the phase variations are smooth, but do vary
Line 446: Line 434:
channels to prevent introducing delay-based closure errors as can happen from averaging over
channels to prevent introducing delay-based closure errors as can happen from averaging over
non-bandpass corrected channels with large phase variations.
non-bandpass corrected channels with large phase variations.


Since the bandpass calibrator is quite strong we do the phase-only
Since the bandpass calibrator is quite strong we do the phase-only
Line 453: Line 440:
Remember that previously we determined the '''myTau''' variable as a list of opacities. It's repeated here but may still be in your CASA session.  
Remember that previously we determined the '''myTau''' variable as a list of opacities. It's repeated here but may still be in your CASA session.  


[[Image:Prebp_phasecal2.png|thumb|Phase only calibration before bandpass. The 4 lines are both polarizations in both spw, unfortunately two of them get the same color green at the moment.]]
[[Image:Prebp_phasecal2.png|thumb|Phase only calibration before bandpass. The 4 lines are both polarizations in both spw, unfortunately two of them get the same color green.]]
<source lang="python">
<source lang="python">
# In CASA
# In CASA
Line 464: Line 451:
</source>
</source>


Plot the solutions
Plot the solutions:


<source lang="python">
<source lang="python">
Line 474: Line 461:
These solutions will appear in the CASA plotter gui. If you closed it after plotting the antennas 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.
These solutions will appear in the CASA plotter gui. If you closed it after plotting the antennas 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.


Next we can apply this phase solution on the fly while determining
Now we can apply this phase solution on the fly while determining
the bandpass solutions on the timescale of the bandpass calibrator scan (solint='inf').  
the bandpass solutions on the timescale of the bandpass calibrator scan (solint='inf').  


We also use the opacity list now instead of myTau - but both options will work.
We also use the opacity list now instead of myTau -- but both options will work.
 
[[Image:Bandpass_amp.png|thumb|Amplitude Bandpass solutions]]
[[Image:Bandpass_phase1.png|thumb|Phase Bandpass solutions]]


<source lang="python">
<source lang="python">
Line 493: Line 477:


The use of  solint='inf' in {{bandpass}} will derive one bandpass
The use of  solint='inf' in {{bandpass}} will derive one bandpass
solution for the whole J1229+0203 scan. Note that if there had been two observations of the bandpass calibrator (for example), this command would have combined the data from both scans to form one bandpass solution, because the default of the combine parameter '''for {{bandpass}}''' is combine='scan'. To solve for one bandpass for each bandpass calibrator scan you would also need to include 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 combine='field' (combine='spw'), the latter two are not generally good ideas for bandpass solutions.   
solution for the whole J1229+0203 scan. Note that if there had been two observations of the bandpass calibrator (for example), this command would have combined the data from both scans to form one bandpass solution, because the default of the combine parameter for '''{{bandpass}}''' is combine='scan'. To solve for one bandpass for each bandpass calibrator scan you would also need to include 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 combine='field' (combine='spw'), the latter two are not generally good ideas for bandpass solutions.   


Plot the solutions, amplitude and phase:
Plot the solutions, amplitude and phase:
[[Image:Bandpass_amp.png|thumb|Amplitude Bandpass solutions]]
[[Image:Bandpass_phase1.png|thumb|Phase Bandpass solutions]]


<source lang="python">
<source lang="python">
Line 509: Line 495:
</source>
</source>


Note that phases on ea12 look noiser than on other antennas. This
Note the scale for ea02 (our reference antenna), and that the phases on ea12 look noiser than on other antennas. This jitter could indicate bad pointing since ea12 had just come back in the array.
jitter could indicate bad pointing; note that ea12 had just come back in
the array.


This step isn't necessary from a calibration perspective, but if you
This step isn't necessary from a calibration perspective, but if you
want to go ahead and check the bandpass calibration on the bandpass
want to go ahead and check the bandpass calibration on the bandpass
calibrator you can run {{applycal}} here. In future we hope to plot
calibrator you can run {{applycal}} here. In the future we hope to plot
corrected data on-the-fly without this {{applycal}} step. Later applycals
corrected data on-the-fly without this {{applycal}} step. Subsequent applycals
will overwrite this one, so no need to worry.
will overwrite this one, so no need to worry.


Line 528: Line 512:
</source>
</source>


Similar to {{flagdata}}, {{applycal}} works like a matrix. The first entries in the 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. In the above example, 'antpos.cal' from any field is applied to source '5', and the 'bandpass.cal' that was obtained for field '5' (the bandpass observation) is also applied to field '2'. Again, lists within the lists are fine.   
Similar to {{flagdata}}, {{applycal}} works like a matrix. The first entries in the 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. In the above example, 'antpos.cal' from any field is applied to source '5', and the 'bandpass.cal' that was obtained for field '5' (the bandpass observation) is also applied to field '5'. Again, lists within the lists are fine.   


<source lang="python">
<source lang="python">
Line 548: Line 532:
==Gain Calibration==
==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, we will solve for them separately. In particular, if the phase changes significantly over a scan time, the amplitude would be decorrelated, if the un-corrected phase were averaged over this timescale. Note that we re-solve for the gain solutions of the bandpass calibrator, so we can derive new solutions that are corrected for the bandpass shape. Since the bandpass calibrator will not be used again, this is not strictly necessary, but is useful to check its calibrated flux density for example. We use a minimum signal-to-noise of 2 here as it seems to be a good compromise for using good data without rejecting too many solutions (minsnr=2).  
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, we will solve for them separately. In particular, if the phase changes significantly over a scan time, the amplitude would be decorrelated, if the un-corrected phase were averaged over this timescale. Note that we re-solve for the gain solutions of the bandpass calibrator, so we can derive new solutions that are corrected for the bandpass shape. Since the bandpass calibrator will not be used again, this is not strictly necessary, but is useful to check its calibrated flux density for example. We use a minimum signal-to-noise of 2 here as it seems to be a good compromise for using good data without rejecting too many solutions (minsnr = 2).  


<source lang="python">
<source lang="python">
Line 562: Line 546:
Here solint='int' coupled with calmode='p' will derive a single phase solution for each 10 second integration. 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.
Here solint='int' coupled with calmode='p' will derive a single phase solution for each 10 second integration. 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.


Note that quite a few solutions are rejected due to SNR<2 (printed to terminal). For the most part it  
Note that quite a few 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.
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 look at the phase solution, and note the obvious scatter within a scan time.
Now look at the phase solution, and note the obvious scatter within a scan time:


<source lang="python">
<source lang="python">
Line 575: Line 559:
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,''' for gaincal, the default of the combine parameter is combine='''' ''''.
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,''' for gaincal, the default of the combine parameter is combine='''' ''''.
[[Image:allcal_phaseinf2.png|thumb|Plot of phase solutions on a scan time.]]
[[Image:allcal_phaseinf2.png|thumb|Plot of phase solutions on a scan time.]]
<source lang="python">
<source lang="python">
# In CASA
# In CASA
Line 611: Line 596:
# In CASA       
# In CASA       
plotcal(caltable='amp.gcal',xaxis='time',yaxis='phase',
plotcal(caltable='amp.gcal',xaxis='time',yaxis='phase',
         iteration='antenna',subplot=331)
         iteration='antenna',subplot=331,plotrange=[-1,-1,-20,20])
</source>
</source>


Indeed, both antenna ea12 (all times) and ea23 (first 1/3 of observation) show particularly large residual phase noise.
Note that we have restricted the plot range so that all antennas are plotted with the same scale (instead of autoscaled; the "-1,-1" values for the x-axis allow autoscaling to take place still).  We can see that both antenna ea12 (all times) and ea23 (first 1/3 of observation) show particularly large residual phase noise.
[[Image:allcal_amp.png|thumb|Plot of amplitude solutions on a scan time]]
[[Image:allcal_amp.png|thumb|Plot of amplitude solutions on a scan time]]
<source lang="python">
<source lang="python">
# In CASA
# In CASA
plotcal(caltable='amp.gcal',xaxis='time',yaxis='amp',
plotcal(caltable='amp.gcal',xaxis='time',yaxis='amp',
         iteration='antenna',subplot=331)
         iteration='antenna',subplot=331,plotrange=[-1,-1,0,1])
</source>
</source>


Line 636: Line 621:


<pre style="background-color: #fffacd;">
<pre style="background-color: #fffacd;">
Flux density for J0954+1743 in SpW=0 is: 0.237135 +/- 0.000858511 (SNR = 276.216, nAnt= 19)
Flux density for J0954+1743 in SpW=0 is: 0.23717 +/- 0.0312543 (SNR = 7.58838, N= 38)
Flux density for J0954+1743 in SpW=1 is: 0.247597 +/- 0.00063498 (SNR = 389.928, nAnt= 19)
Flux density for J0954+1743 in SpW=1 is: 0.240683 +/- 0.0302257 (SNR = 7.96288, N= 38)
Flux density for J1229+0203 in SpW=0 is: 26.508 +/- 0 (SNR = inf, nAnt= 19)
Flux density for J1229+0203 in SpW=0 is: 28.2803 +/- 0.245907 (SNR = 115.004, N= 38)
Flux density for J1229+0203 in SpW=1 is: 26.2517 +/- 0 (SNR = inf, nAnt= 19)
Flux density for J1229+0203 in SpW=1 is: 28.4019 +/- 0.257634 (SNR = 110.241, N= 38)
 
</pre>
</pre>


<!--
Obviously, the signal-to-noise for J1229+0203 can't be infinity! This is just an indication that there is only one scan for this source, and we derived a scan based amplitude solution, so there is no variation to calculate.  
Obviously, the signal-to-noise for J1229+0203 can't be infinity! This is just an indication that there is only one scan for this source, and we derived a scan based amplitude solution, so there is no variation to calculate.  
-->


Next, check that the flux.cal table looks as expected.
Next, check that the flux.cal table looks as expected:


<source lang="python">
<source lang="python">
Line 652: Line 638:
         iteration='antenna',subplot=331)
         iteration='antenna',subplot=331)
</source>
</source>
Note the odd behavior for antennas ea07, ea12, and ea23.


==Applycal and Inspect==
==Applycal and Inspect==
Line 657: Line 645:
Now we apply the calibration to each source, according to which tables are appropriate, and which source should be used to do that particular calibration. For the calibrators, all bandpass solutions come from the bandpass calibrator (id=5), and the phase and amplitude calibration comes from their own solutions.  
Now we apply the calibration to each source, according to which tables are appropriate, and which source should be used to do that particular calibration. For the calibrators, all bandpass solutions come from the bandpass calibrator (id=5), and the phase and amplitude calibration comes from their own solutions.  


'''Note:''' In applycal we set calwt=F. 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. During the VLA era, meaningful weights were available for each visibility. However, EVLA is 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=T when the input weights are themselves not meaningful, especially for self-calibration on resolved sources (your flux calibrator and target for example). In a few months EVLA data will again have meaningful weights and the default calwt=T will likely again be the best option.
'''Note:''' In applycal we set calwt=F. 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. During the VLA era, meaningful weights were available for each visibility. However, at the time of this observation, the EVLA 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=T when the input weights are themselves not meaningful, especially for self-calibration on resolved sources (your flux calibrator and target, for example).  
 
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.


<source lang="python">
<source lang="python">
Line 738: Line 728:
</source>
</source>


it may be a 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 initial and second pointing. Recall that these are antennas we became suspicious of while inspecting the calibration solutions.
It may be a 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 initial and second pointing. Recall that these are antennas we became suspicious of while inspecting the calibration solutions.


[[Image:target_uvdist.png|thumb|IRC+10216 as a function of uv-distance.]]
[[Image:target_uvdist.png|thumb|IRC+10216 as a function of uv-distance.]]
Now lets see how the target looks. Because the target has resolved structure, its best to look at it as
Now let's see how the target looks. Because the target has resolved structure, its best to look at it as
a function of uvdistance. We'll go ahead and exclude the three antennas we already know have problems.
a function of uvdistance. We'll go ahead and exclude the three antennas we already know have problems.


Line 751: Line 741:
</source>
</source>


you can see that the spikes
The color indicates that the spikes are caused by a single antenna. Use, zoom, mark, and locate to see which one.
are caused by a single antenna. Use, zoom, mark, and locate to see which one.
Also look at spw=1.
Also look at spw=1.


Turns out to be ea28; to confirm, replot with antenna=!ea28, and
The bad antenna turns out to be ea28; to confirm, replot with antenna=!ea28:


<source lang="python">
<source lang="python">
Line 764: Line 753:
</source>
</source>


To see if it's restricted to a certain time
Looks much better!  To see if it's restricted to a certain time, do:


<source lang="python">
<source lang="python">
Line 772: Line 761:
       avgchannel='64',spw='0:4~60',antenna='ea28', coloraxis='antenna1')
       avgchannel='64',spw='0:4~60',antenna='ea28', coloraxis='antenna1')
</source>
</source>


Baselines with ea28 clearly show issues until about two-thirds of the way through the observation.  
Baselines with ea28 clearly show issues until about two-thirds of the way through the observation.  
Plot another distant antenna to compare. We will go ahead and flag it all, since its hanging far out on the north
Plot another distant antenna to compare. We will go ahead and flag it all, since it's hanging far out on the north
arm by itself.
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
The additional data we've identified as bad need to be flagged, and then all the calibration steps will need to be run again.
again.


<source lang="python">
<source lang="python">
Line 801: Line 788:
         gaintable=['antpos.cal'],
         gaintable=['antpos.cal'],
         opacity=[0.0410, 0.0408],gaincurve=T)
         opacity=[0.0410, 0.0408],gaincurve=T)
</source>
#
 
<source lang="python">
# In CASA
bandpass(vis='day2_TDEM0003_10s_norx',caltable='bandpass_redo.bcal',
bandpass(vis='day2_TDEM0003_10s_norx',caltable='bandpass_redo.bcal',
         field='5',
         field='5',
Line 810: Line 794:
         gaintable=['antpos.cal','bpphase_redo.gcal'],
         gaintable=['antpos.cal','bpphase_redo.gcal'],
         opacity=[0.0410, 0.0408],gaincurve=T)
         opacity=[0.0410, 0.0408],gaincurve=T)
</source>
#
 
<source lang="python">
# In CASA
gaincal(vis='day2_TDEM0003_10s_norx',caltable='intphase_redo.gcal',
gaincal(vis='day2_TDEM0003_10s_norx',caltable='intphase_redo.gcal',
         field='2,5,7',spw='0~1:4~60',
         field='2,5,7',spw='0~1:4~60',
Line 819: Line 800:
         gaintable=['antpos.cal','bandpass_redo.bcal'],
         gaintable=['antpos.cal','bandpass_redo.bcal'],
         opacity=[0.0410, 0.0408],gaincurve=T)
         opacity=[0.0410, 0.0408],gaincurve=T)
</source>
#
 
<source lang="python">
# In CASA
gaincal(vis='day2_TDEM0003_10s_norx',caltable='scanphase_redo.gcal',
gaincal(vis='day2_TDEM0003_10s_norx',caltable='scanphase_redo.gcal',
         field='2,5,7',spw='0~1:4~60',
         field='2,5,7',spw='0~1:4~60',
Line 828: Line 806:
         gaintable=['antpos.cal','bandpass_redo.bcal'],
         gaintable=['antpos.cal','bandpass_redo.bcal'],
         opacity=[0.0410, 0.0408],gaincurve=T)
         opacity=[0.0410, 0.0408],gaincurve=T)
</source>
#
 
<source lang="python">
# In CASA
gaincal(vis='day2_TDEM0003_10s_norx',caltable='amp_redo.gcal',
gaincal(vis='day2_TDEM0003_10s_norx',caltable='amp_redo.gcal',
         field='2,5,7',spw='0~1:4~60',
         field='2,5,7',spw='0~1:4~60',
Line 837: Line 812:
         gaintable=['antpos.cal','bandpass_redo.bcal','intphase_redo.gcal'],
         gaintable=['antpos.cal','bandpass_redo.bcal','intphase_redo.gcal'],
         opacity=[0.0410, 0.0408],gaincurve=T)
         opacity=[0.0410, 0.0408],gaincurve=T)
</source>
#
<source lang="python">
# In CASA     
fluxscale(vis='day2_TDEM0003_10s_norx',caltable='amp_redo.gcal',
fluxscale(vis='day2_TDEM0003_10s_norx',caltable='amp_redo.gcal',
           fluxtable='flux_redo.cal',reference='7')
           fluxtable='flux_redo.cal',reference='7')
Line 846: Line 818:


<pre style="background-color: #fffacd;">
<pre style="background-color: #fffacd;">
Flux density for J0954+1743 in SpW=0 is: 0.235345 +/- 0.000879422
Flux density for J0954+1743 in SpW=0 is: 0.265112 +/- 0.0286296 (SNR = 9.26005, N= 32)
    (SNR = 267.613, nAnt= 16)
Flux density for J0954+1743 in SpW=1 is: 0.265065 +/- 0.0292281 (SNR = 9.06882, N= 32)
Flux density for J0954+1743 in SpW=1 is: 0.241996 +/- 0.000930228
Flux density for J1229+0203 in SpW=0 is: 28.4587 +/- 0.286533 (SNR = 99.3209, N= 32)
    (SNR = 260.147, nAnt= 16)
Flux density for J1229+0203 in SpW=1 is: 28.335 +/- 0.29411 (SNR = 96.3414, N= 32)
Flux density for J1229+0203 in SpW=0 is: 25.2479 +/- 0  
    (SNR = inf, nAnt= 16)
Flux density for J1229+0203 in SpW=1 is: 24.9907 +/- 0
    (SNR = inf, nAnt= 16)
 
Flux density for J0954+1743 in SpW=0 is: 0.247052 +/- 0.000946345 (SNR = 261.059, nAnt= 16)
Flux density for J0954+1743 in SpW=1 is: 0.254038 +/- 0.00097531 (SNR = 260.47, nAnt= 16)
Flux density for J1229+0203 in SpW=0 is: 26.5079 +/- 0 (SNR = inf, nAnt= 16)
Flux density for J1229+0203 in SpW=1 is: 26.2335 +/- 0 (SNR = inf, nAnt= 16)
</pre>
</pre>


Line 873: Line 836:
         gainfield=['','5','2','2'],
         gainfield=['','5','2','2'],
         opacity=[0.0410, 0.0408],gaincurve=T,calwt=F)
         opacity=[0.0410, 0.0408],gaincurve=T,calwt=F)
</source>
#
 
<source lang="python">
# In CASA
applycal(vis='day2_TDEM0003_10s_norx',field='5',
applycal(vis='day2_TDEM0003_10s_norx',field='5',
         gaintable=['antpos.cal','bandpass_redo.bcal','intphase_redo.gcal','flux_redo.cal'],
         gaintable=['antpos.cal','bandpass_redo.bcal','intphase_redo.gcal','flux_redo.cal'],
         gainfield=['','5','5','5'],
         gainfield=['','5','5','5'],
         opacity=[0.0410, 0.0408],gaincurve=T,calwt=F)
         opacity=[0.0410, 0.0408],gaincurve=T,calwt=F)
</source>
#
 
<source lang="python">
# In CASA
applycal(vis='day2_TDEM0003_10s_norx',field='7',
applycal(vis='day2_TDEM0003_10s_norx',field='7',
         gaintable=['antpos.cal','bandpass_redo.bcal','intphase_redo.gcal','flux_redo.cal'],
         gaintable=['antpos.cal','bandpass_redo.bcal','intphase_redo.gcal','flux_redo.cal'],
         gainfield=['','5','7','7'],
         gainfield=['','5','7','7'],
         opacity=[0.0410, 0.0408],gaincurve=T,calwt=F)
         opacity=[0.0410, 0.0408],gaincurve=T,calwt=F)
</source>
#
 
[[Image:gaincal_corrflag.png|thumb|Gain calibrator after further flagging and recalibration]]
[[Image:target_corrflag.png|thumb|IRC+10216 after further flagging and recalibration (after selecting colorize by spw).]]
 
<source lang="python">
# In CASA
applycal(vis='day2_TDEM0003_10s_norx',field='3',
applycal(vis='day2_TDEM0003_10s_norx',field='3',
         gaintable=['antpos.cal','bandpass_redo.bcal','scanphase_redo.gcal','flux_redo.cal'],
         gaintable=['antpos.cal','bandpass_redo.bcal','scanphase_redo.gcal','flux_redo.cal'],
Line 901: Line 852:
         opacity=[0.0410, 0.0408],gaincurve=T,calwt=F)
         opacity=[0.0410, 0.0408],gaincurve=T,calwt=F)
</source>
</source>
[[Image:gaincal_corrflag.png|thumb|Gain calibrator after further flagging and recalibration]]
[[Image:target_corrflag.png|thumb|IRC+10216 after further flagging and recalibration (after selecting colorize by spw).]]


Now you can inspect the calibrated data again. Except for random scatter things look pretty good.
Now you can inspect the calibrated data again. Except for random scatter things look pretty good.
Line 930: Line 884:
split(vis='day2_TDEM0003_10s_norx',outputvis='J0954',
split(vis='day2_TDEM0003_10s_norx',outputvis='J0954',
       field='2')
       field='2')
</source>
#
 
<source lang="python">
# In CASA
split(vis='day2_TDEM0003_10s_norx',outputvis='IRC10216',
split(vis='day2_TDEM0003_10s_norx',outputvis='IRC10216',
       field='3')
       field='3')
Line 943: Line 894:
# In CASA
# In CASA
clearcal(vis='J0954')
clearcal(vis='J0954')
</source>
#
 
<source lang="python">
# In CASA
clearcal(vis='IRC10216')
clearcal(vis='IRC10216')
</source>
</source>
Line 956: Line 904:


--[[User:Cbrogan|Crystal Brogan]]     
--[[User:Cbrogan|Crystal Brogan]]     
--additions:[[User:Jott| Juergen Ott]]
--additions:[[User:Jott| Juergen Ott]], [[User:Mkrauss| Miriam Krauss]], [[User:Dshepher| Deb Shepherd]]
 
{{Checked 3.3.0}}

Latest revision as of 18:47, 17 May 2012


This is the first part of the EVLA IRC+10216 data reduction tutorial. CASA 3.4 version


This is the first part of the EVLA IRC+10216 data reduction tutorial. IRC+10216 Part 2: continuum subtraction, imaging, selfcal

Overview

VLT V-band image of IRC+10216 showing dust rings out to a radius of 90" by Leão et al. (2006, A&A, 455, 187).

This tutorial describes the data reduction for two spectral lines observed toward the AGB star IRC+10216. In this EVLA "OSRO1" mode observation one subband was observed in each of two basebands, with the subbands centered on the HC3N and SiS lines near 36 GHz. The raw data were loaded into CASA with importevla, where zero and shadowed data were flagged (note that this is done by the NRAO Data Archive). Then the data were split, so we could average from the native 1-second integrations to 10 seconds, select only antennas with Ka-band receivers, and select only spectral windows (called spw in CASA) with Ka-band data. This produces a significantly smaller dataset for processing.

IRC+10216 (CW Leo) is the brightest star in the sky at 5 microns but only 16th magnitude visually. It was discovered during the first survey of the infrared sky, carried out by Bob Leighton and Gerry Neugebauer in 1965. An Asymptotic Giant Branch star, it 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 emission; the continuum emission seen at radio wavelengths probes the actual 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 material drops, the chemistry freezes. But the molecules continue their long coast outward into the Galaxy, and as the shell thins ultraviolet light from the ambient galactic radiation field penetrates and initiates a 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. The result is the centrally condensed emission we see in the present observations. Interestingly, in the next lower transition at 18 GHz, the line shape is much different from what we see here. At the extreme velocities in the profile, very bright narrow emission is seen which has been interpreted as maser emission. The interested student can find EVLA observations of this line in the archive.

HC3N, a much more complex species, is created by the photochemistry which becomes active as atoms and pieces of molecules destroyed by ultraviolet radiation undergo the next phase of chemistry in the shell. HC3N has many vibrational modes which may be excited in addition to its rotational modes. 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, but during its brief existence its rise to abundance in the envelope results in a ring of emission, which is what is observed in this image made with the EVLA. A recent model by Cordiner & Millar (2009, ApJ, 697, 68) describes a new chemical model for the shell, which also takes into account the variation of mass loss by the star. They show that in addition to purely chemical effects, local gas and dust density peaks play a role in shaping the observed emission.

Getting the data

The post-split averaged data can be downloaded from http://casa.nrao.edu/Data/EVLA/IRC10216/day2_TDEM0003_10s_norx.tar.gz (data size: 1.1GB)

Once the download is complete, unzip and unpack the file:

# in a terminal, outside of CASA:
tar -xzvf day2_TDEM0003_10s_norx.tar.gz

How to Use This casaguide

Inputs from one of the clean commands from this tutorial

There are a number of possible ways to run CASA. Many aspects are described in Getting_Started_in_CASA. You should review this page if you are new to CASA. In brief you can run CASA interactively by looking at the inputs to tasks with inp taskname (example: inp clean), setting the parameters one by one (example: selectdata=T) as you desire and then run go. After setting parameters one by one in a task and then looking at the inputs again, you will notice that the parameters that have been set to something other than their defaults are blue. If you have mistyped any parameters, they will be red and must be fixed for the task to run correctly. You can get more detailed help on any task by typing help taskname (example: help clean). Once a task is run you can get the same parameters back by running tget taskname (example: tget clean); subsequent runs will overwrite the previous tget file.

The second way to run CASA is to provide task function calls. This tutorial is made up of such calls, which were developed by looking at the inputs for each task and deciding what needed to be changed from default values. For task function calls, only parameters that you want to be different from their defaults need to be set. A series of task function calls can be combined together into a script, and run with execfile('scriptname.py'). It is possible to extract a script containing all the CASA task function calls in this and other casaguides using the method described at the Extracting_scripts_from_these_tutorials page.

If you are a relative novice or just new to CASA it is strongly recommended to work through this tutorial by cutting and pasting the task function calls provided below after you have read all the associated explanations. Work at your own pace, look at the inputs to the tasks to see what other options exist, and read the help files. Later, when you are more comfortable, you might try to extract the script, modify it for your purposes, and begin to reduce other data.

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

For all EVLA observations, the operators keep an observing log. You can look at the observing logs at the observing log website by entering a date range covering only 26 Apr 2010, then selecting "TDEM0003". Pertinent information from this observation is repeated below:

INFORMATION FROM OBSERVING LOG:
Date of the observation: 26-April-2010
There are no Ka-band receivers on ea11, ea13, ea14, ea16, ea17, ea18, ea26  
Antenna ea06 is out of the array
Antenna 15 has a focus timeout
Antennas ea10, ea12, ea22 do not have good baseline positions

As mentioned in the log, 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, we need to insert these corrections by hand using gencal. The resulting table will need to be supplied as a "prior" calibration to all subsequent calibration steps. The corrections can be ascertained from the EVLA/VLA Baseline Corrections page. Be sure to carefully read to the bottom of that webpage to see how to calculate the additive corrections. The current case is simple as there is only a single correction for each antenna. In the future, CASA will include an automated lookup of the corrections.

# In CASA
gencal(vis='day2_TDEM0003_10s_norx',caltable='antpos.cal',
       caltype='antpos',
       antenna='ea12,ea22',
       parameter=[-0.0072,0.0045,-0.0017, -0.0220,0.0040,-0.0190])

Please note: if you are reducing VLA data taken before March 1, 2010, you need to set caltype='antposvla'.

There are two other types of "priors" that we will need to apply at each calibration stage described below:

  1. Gaincurve -- the gaincurve describes how each antenna behaves as a function of elevation, for each receiver band. Currently only gaincurves for the VLA/EVLA are available (see [1] for the incorporated models). This option should not be used with any other telescopes.
  2. Opacity correction -- the opacity of the observation can be computed from a seasonal model and/or weather station information. We are planning to have a task available for this information. At the moment, the CASA_EVLA_Scripts page hosts a script and a contributed task to display the weather information and to calculate the zenith opacities for each spectral window. After the zenith opacities are derived, they will be recomputed for the correct elevation of the data automatically using [math]\displaystyle{ e^{(-\csc[el]\tau_z)} }[/math] in gaincal, applycal, bandpass etc.

Now we want to plot the opacity of the atmosphere at 36 GHz at the time this observation was taken so it can be corrected for in subsequent calibrations. There is not currently a task in CASA to do this however, an external task, plotWX.py, has been written that will do this. In the following we will download the contributed task, and bring it into CASA. We will then use the results to determine the opacity of the atmosphere at the time these observations were taken.

Download these two files:

File:Task plotWX.py and File:PlotWX.xml. The first file contains the actual python code and the second the interface and inline help for CASA. Download the files and note that for obscure reasons the wiki will always capitalize the uploaded files abut we need to change this. In a regular terminal:

# in a terminal, outside of CASA:
wget http://casaguides.nrao.edu/images/e/e4/PlotWX.xml
mv PlotWX.xml plotWX.xml
wget http://casaguides.nrao.edu/images/2/24/Task_plotWX.py
mv Task_plotWX.py task_plotWX.py

Now, we go back to CASA and build the task itself with buildmytask (it is a shell command and not a CASA command, but we will run it from within CASA; it will create a few *.pyc and *cli files as well as mytasks.py). And then add it into your CASA session by running mytasks.py:

# In CASA
!buildmytasks
execfile('mytasks.py')

Now your brand-new plotWX task is available:

CASA <7>: inp plotWX
--------> inp(plotWX)
#  plotWX :: Plot elements of the weather table for a given MS
vis                 =         ''        #  MS name
seasonal_weight     =        0.5        #  weight of the seasonal model
doPlot              =       True        #  set this to True to create a plot
async               =      False        #  If true the taskname must be started using plotWX(...)

seasonal_weight is a parameter that combines global seasonal weather conditions (measured over many years) with the weather information from the weather station obtained during the observations. 0.5 will combine both with equal weight.

plotWX weather table figure

We will be running plotWX in a way that will assign the opacity list (one entry for each spw in ascending order) to the variable myTau:

# In CASA
myTau = plotWX(vis='day2_TDEM0003_10s_norx', doPlot=T)

With this input, the task

  • returns the opacity values to the logger:
wrote weather figure: day2_TDEM0003_10s_norx.plotWX.png
SPW : Frequency (GHz) : Zenith opacity (nepers)
0  :   36.387  :  0.041
1  :   36.305  :  0.041
  • creates a file "day2_TDEM0003_10s_norx.plotWX.png" with the elevation of the sun, the wind speed and direction, the temperature, 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 assigns the myTau variable to the list of opacities per spectral window:
# In CASA
myTau
[0.041019209411983663, 0.040779609355637694]

The calibration task opacity parameters can now take either the variable directly opacity = myTau or the list of opacities like opacity = [0.0410, 0.0408] or, if a single value is good enough to correct for all spectral windows, one can simple use opacity = 0.041.

Initial Inspection and Flagging

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:

# In CASA
listobs(vis='day2_TDEM0003_10s_norx')

Below we have cut and pasted the most relevant output from the logger.

Fields: 4
  ID   Code Name         RA            Decl           Epoch   SrcId nVis   
  2    D    J0954+1743   09:54:56.8236 +17.43.31.2224 J2000   2     65326  
  3    NONE IRC+10216    09:47:57.3820 +13.16.40.6600 J2000   3     208242 
  5    F    J1229+0203   12:29:06.6997 +02.03.08.5982 J2000   5     10836  
  7    E    J1331+3030   13:31:08.2880 +30.30.32.9589 J2000   7     5814   
   (nVis = Total number of time/baseline visibilities per field) 
Spectral Windows:  (2 unique spectral windows and 1 unique polarization setups)
  SpwID  #Chans Frame Ch1(MHz)    ChanWid(kHz)TotBW(kHz)  Ref(MHz)    Corrs           
  0          64 TOPO  36387.2295  125         8000        36387.2295  RR  RL  LR  LL  
  1          64 TOPO  36304.542   125         8000        36304.542   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.         
  0    ea01  W09       25.0 m   -107.37.25.2  +33.53.51.0  
  1    ea02  E02       25.0 m   -107.37.04.4  +33.54.01.1  
  2    ea03  E09       25.0 m   -107.36.45.1  +33.53.53.6  
  3    ea04  W01       25.0 m   -107.37.05.9  +33.54.00.5  
  4    ea05  W08       25.0 m   -107.37.21.6  +33.53.53.0  
  5    ea07  N06       25.0 m   -107.37.06.9  +33.54.10.3  
  6    ea08  N01       25.0 m   -107.37.06.0  +33.54.01.8  
  7    ea09  E06       25.0 m   -107.36.55.6  +33.53.57.7  
  8    ea12  E08       25.0 m   -107.36.48.9  +33.53.55.1  
  9    ea15  W06       25.0 m   -107.37.15.6  +33.53.56.4  
  10   ea19  W04       25.0 m   -107.37.10.8  +33.53.59.1  
  11   ea20  N05       25.0 m   -107.37.06.7  +33.54.08.0  
  12   ea21  E01       25.0 m   -107.37.05.7  +33.53.59.2  
  13   ea22  N04       25.0 m   -107.37.06.5  +33.54.06.1  
  14   ea23  E07       25.0 m   -107.36.52.4  +33.53.56.5  
  15   ea24  W05       25.0 m   -107.37.13.0  +33.53.57.8  
  16   ea25  N02       25.0 m   -107.37.06.2  +33.54.03.5  
  17   ea27  E03       25.0 m   -107.37.02.8  +33.54.00.5  
  18   ea28  N08       25.0 m   -107.37.07.5  +33.54.15.8  

In addition to source names, antenna names, etc. CASA indexes each of these parameters starting from 0. To specify, e.g., a source field, one can use the index or the name -- both are accepted.

Note that 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 EVLA Observing Tool 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.

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
Antenna locations from running plotants

Look at a graphical plot of the antenna locations and save hardcopy in case you want it later. This will be useful for picking a reference antenna -- typically a good choice is an antenna close to the center of the array. Unless it shows problems after inspection of the data, we provisionally choose ea02.

Elevation as a function of time (after selecting colorize by field).
# In CASA
plotants(vis='day2_TDEM0003_10s_norx',figfile='ant_locations.png')

Next, let's look at the elevation as a function of time for all sources. It's not the case for these data, but if the elevation is very low (usually at start or end of track) you may want to flag. Also, how near in elevation your flux calibrator is to your target will impact your ultimate absolute flux calibration accuracy.

Unfortunately, the target and flux calibrator are not particularly well-matched for this observation, as you can show by plotting the elevation for each source (each sources has a different colors). We will be using data in spectral window 0, channels 4 to 60 for this plot as given by the spw = '0:4~60' parameter; the general CASA selection syntax is described in the Measurement Selection Syntax Document:

# In CASA
plotms(vis='day2_TDEM0003_10s_norx',
       xaxis='time',yaxis='elevation',correlation='RR,LL',
       avgchannel='64',spw='0:4~60', coloraxis='field')

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!)

Result of plotms
Zooming in and marking region (hatched box)

Next, let's look at all the source amplitudes as a function of time:

# In CASA
plotms(vis='day2_TDEM0003_10s_norx',
       xaxis='time',yaxis='amp',correlation='RR,LL',
       avgchannel='64',spw='0:4~60', coloraxis='field')

Now zoom in on the region very near zero amplitude for sources J0954+1743 and IRC+10216. To zoom, select the Zoom tool Casaplotms-zoom.png in lower left corner of the plotms GUI, then you can left click to draw a box. Look for the low values (you may want to zoom a few times to really see the suspect points clearly). Now use the Mark Region MarkRegionsButton.png and Locate Casaplotms-locate-tool.png buttons (located along the bottom of the GUI) to see which antenna is causing problems. The output is be shown in the logger. Since all the "located" baselines include ea12, this is the responsible antenna.

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 "force reload" box in the lower left, or do Shift+Plot.

* If the plotms tool hangs during a plot, try clicking the cancel button on the 
load progress GUI, and/or if you see a "table locked" message try typing 
"clearstat" on the CASA command line.

* Occasionally, plotms will get into a strange state that you cannot clear from inside. 
We are working on these issues, but for now, when all else fails, exit from casapy and 
restart.  

Now click the clear region Casaplotms-clear-region.png button, and 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 next pointing run: 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 timerange, 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' for antenna, which removes ea07 from the plot (in CASA selection, "!" deselects).

We can set up a flagging command to get both bad antennas for the appropriate time and spw:

# In CASA
flagdata(vis='day2_TDEM0003_10s_norx',
         field=['2,3','2,3'],
         spw=['','1'],
         antenna=['ea12','ea07'],
         timerange=['03:41:00~04:10:00','03:21:40~04:10:00'])

flagdata works by spanning up a matrix. The first entries in each list must be taken as one flagging command, as well as the second entries etc. Lists within lists are fine. In the above example, the first flagging command is issued for fields 2 and 3 for all spws and within the 03:41:00~04:10:00 timerange. A second command is again for the fields 2 and 3 but for spw 1 only and for the second timerange in the list '03:21:40~04:10:00'.

Note that because the chosen timerange is limited to fields 2 and 3, the field parameter is not really needed; however, flagdata will run fastest if you put as many constraints as possible (since it will traverse less data).

Now remove the !ea07 from antenna 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!).

Amplitude vs. uv-distance for IRC+10216, both spw

Let's look more closely at IRC+10216:

# In 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')

You can see a that there are some noisy high points. But now try:

# In 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')

Now, we can see that most of the high points on IRC+10216 are due to large scale emission on short baselines, 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.

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, one typically used a point source (constant flux) model for the flux calibrator, possibly with a uvrange cutoff if necessary. More recently for the VLA/EVLA, 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).

Currently, CASA contains models for the most common calibrators at the most common frequencies but not yet the new EVLA frequency bands S and Ka. These will be added as soon as they become available. One may check the availability of calibration models in setjy:

# In CASA
setjy(vis='day2_TDEM0003_10s_norx',listmodimages=T)

The terminal will now show the available models, e.g. 3C286_C.im, 3C48_K.im etc. (setjy will search in the working directory for images that may contain models, as well as in a CASA directory where known calibrator models are stored.)

As mentioned above, the Ka band does not yet have a model incorporated at this time as the full configuration cycle has not ended yet. However, models for Ka-band that are good for D, C, and B configurations are available from the EVLA calibrator model webpage for the new frequency bands.

From this page, please download the file 36.68 GHz CASA image of 3C286, this is the direct link. To do so, you may run

# in a terminal, outside of CASA:
wget http://www.aoc.nrao.edu/evlacal/3C286_36.68G.model.tgz
tar xzvf 3C286_36.68G.model.tgz

to create the model 3C286_36.68G.model on disk.

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.

# In CASA
setjy(vis='day2_TDEM0003_10s_norx',field='7',spw='0~1',
      modimage='3C286_36.68G.model')
The logger output for each spw is:
J1331+3030  rawspwid=  0  [I=1.7762, Q=0, U=0, V=0] Jy, (Perley-Butler 2010)
J1331+3030  rawspwid=  1  [I=1.7794, Q=0, U=0, V=0] Jy, (Perley-Butler 2010)

The absolute fluxes for the frequencies have now been determined, and one can proceed to the bandpass and complex gain calibrations.

Bandpass

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 as it seems like a good candidate for the reference antenna.

# In CASA
plotms(vis='day2_TDEM0003_10s_norx',field='5',
       xaxis='channel',yaxis='phase',correlation='RR',
       avgtime='1e8',spw='0:4~60',antenna='ea02&ea23')
Phase as a function of channel for ea02 (after Custom and upping "Style" to 3.)

If you find the points difficult to see, go to the "Display" tab, set "Unflagged Points Symbol" to "Custom", and increase the "Style" number, then replot. The phase variation is modest ~10 degrees. Now expand to all baselines that include ea02, and add an extra dimension of color:

# In 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')
Phase as a function of time for all baselines with antenna ea02 (after Custom and upping "Style" to 3.)

From this you can see that the phase variation across the bandpass is modest. Next check LL, and spw=1, both correlations. Also check other antennas if you like.

Now look at the phase as a function of time:

# In CASA
plotms(vis='day2_TDEM0003_10s_norx',field='5',
       xaxis='time',yaxis='phase',correlation='RR',
       avgchannel='64',spw='0:4~60',antenna='ea02&ea23')

Expand to all antennas with ea02:

# In 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')

You can see that the phase variations are smooth, but do vary significantly over the 5 minutes of observation -- in most cases by a few 10s of degrees. Zoom in to see this better if you want.

The conclusion from this investigation is that you 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, you can average over several channels to increase the signal to noise of the phase vs. time solution. If the phase variation as a function of channel is larger you may need to use 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.

Since the bandpass calibrator is quite strong we do the phase-only solution on the integration time of 10 seconds (solint='int').

Remember that previously we determined the myTau variable as a list of opacities. It's repeated here but may still be in your CASA session.

Phase only calibration before bandpass. The 4 lines are both polarizations in both spw, unfortunately two of them get the same color green.
# In CASA
myTau=[0.041019209411983566, 0.040779609355637569]
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'],
        opacity=myTau,gaincurve=T)

Plot the solutions:

# In CASA
plotcal(caltable='bpphase.gcal',xaxis='time',yaxis='phase',
        iteration='antenna',subplot=331,plotrange=[0,0,-180,180])

These solutions will appear in the CASA plotter gui. If you closed it after plotting the antennas 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.

Now we can apply this phase solution on the fly while determining the bandpass solutions on the timescale of the bandpass calibrator scan (solint='inf').

We also use the opacity list now instead of myTau -- but both options will work.

# In CASA
bandpass(vis='day2_TDEM0003_10s_norx',caltable='bandpass.bcal',field='5',
        refant='ea02',solint='inf',solnorm=T,
        gaintable=['antpos.cal','bpphase.gcal'],
        opacity=[0.041, 0.0408],gaincurve=T)

A few words about solint and combine:

The use of solint='inf' in bandpass will derive one bandpass solution for the whole J1229+0203 scan. Note that if there had been two observations of the bandpass calibrator (for example), this command would have combined the data from both scans to form one bandpass solution, because the default of the combine parameter for bandpass is combine='scan'. To solve for one bandpass for each bandpass calibrator scan you would also need to include 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 combine='field' (combine='spw'), the latter two are not generally good ideas for bandpass solutions.

Plot the solutions, amplitude and phase:

Amplitude Bandpass solutions
Phase Bandpass solutions
# In CASA
plotcal(caltable='bandpass.bcal',xaxis='chan',yaxis='amp',
        iteration='antenna',subplot=331)
# In CASA
plotcal(caltable='bandpass.bcal',xaxis='chan',yaxis='phase',
        iteration='antenna',subplot=331)

Note the scale for ea02 (our reference antenna), and that the phases on ea12 look noiser than on other antennas. This jitter could indicate bad pointing since ea12 had just come back in the array.

This step isn't necessary from a calibration perspective, but if you want to go ahead and check the bandpass calibration on the bandpass calibrator you can run applycal here. In the future we hope to plot corrected data on-the-fly without this applycal step. Subsequent applycals will overwrite this one, so no need to worry.

Phase as a function of channel, plotting the corrected data (after Custom and upping "Style" to 3.)
applycal(vis='day2_TDEM0003_10s_norx',field='5',
        gaintable=['antpos.cal','bandpass.bcal'],
        gainfield=['','5'],
        opacity=[0.0410, 0.0408],gaincurve=T,calwt=F)

Similar to flagdata, applycal works like a matrix. The first entries in the 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. In the above example, 'antpos.cal' from any field is applied to source '5', and the 'bandpass.cal' that was obtained for field '5' (the bandpass observation) is also applied to field '5'. Again, lists within the lists are fine.

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')

Note that the phase and amplitude as a function of channel are very flat now.

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, we will solve for them separately. In particular, if the phase changes significantly over a scan time, the amplitude would be decorrelated, if the un-corrected phase were averaged over this timescale. Note that we re-solve for the gain solutions of the bandpass calibrator, so we can derive new solutions that are corrected for the bandpass shape. Since the bandpass calibrator will not be used again, this is not strictly necessary, but is useful to check its calibrated flux density for example. We use a minimum signal-to-noise of 2 here as it seems to be a good compromise for using good data without rejecting too many solutions (minsnr = 2).

# In CASA
gaincal(vis='day2_TDEM0003_10s_norx',caltable='intphase.gcal',
        field='2,5,7',spw='0~1:4~60',
        refant='ea02',calmode='p',solint='int',minsnr=2.0,
        gaintable=['antpos.cal','bandpass.bcal'],
        opacity=[0.0410, 0.0408],gaincurve=T)
Plot of phase solutions on an integration time.

Here solint='int' coupled with calmode='p' will derive a single phase solution for each 10 second integration. 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.

Note that quite a few 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 look at the phase solution, and note the obvious scatter within a scan time:

# In CASA
plotcal(caltable='intphase.gcal',xaxis='time',yaxis='phase',
        iteration='antenna',subplot=331,plotrange=[0,0,-180,180])

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, for gaincal, the default of the combine parameter is combine=' '.

Plot of phase solutions on a scan time.
# In CASA
gaincal(vis='day2_TDEM0003_10s_norx',caltable='scanphase.gcal',
        field='2,5,7',spw='0~1:4~60',
        refant='ea02',calmode='p',solint='inf',minsnr=2.0,
        gaintable=['antpos.cal','bandpass.bcal'],
        opacity=[0.0410, 0.0408],gaincurve=T)
# In CASA
plotcal(caltable='scanphase.gcal',xaxis='time',yaxis='phase',
        iteration='antenna',subplot=331,plotrange=[0,0,-180,180])

Note that there are no failed solutions here because of the added S/N afforded by the longer solint. Alternatively, instead of making a separate phase solution for application to the target, one can also run smoothcal to smooth the solutions derived on the integration time.

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

# In CASA
gaincal(vis='day2_TDEM0003_10s_norx',caltable='amp.gcal',
        field='2,5,7',spw='0~1:4~60',
        refant='ea02',calmode='ap',solint='inf',minsnr=2.0,
        gaintable=['antpos.cal','bandpass.bcal','intphase.gcal'],
        opacity=[0.0410, 0.0408],gaincurve=T)
Plot of residual phase solutions on a scan time

Now let's look at the resulting phase solutions. 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 you see scatter of more than a few degrees here, you should consider going back and looking for more data to flag, particularly bad timeranges etc.

# In CASA       
plotcal(caltable='amp.gcal',xaxis='time',yaxis='phase',
        iteration='antenna',subplot=331,plotrange=[-1,-1,-20,20])

Note that we have restricted the plot range so that all antennas are plotted with the same scale (instead of autoscaled; the "-1,-1" values for the x-axis allow autoscaling to take place still). We can see that both antenna ea12 (all times) and ea23 (first 1/3 of observation) show particularly large residual phase noise.

Plot of amplitude solutions on a scan time
# In CASA
plotcal(caltable='amp.gcal',xaxis='time',yaxis='amp',
        iteration='antenna',subplot=331,plotrange=[-1,-1,0,1])

Note that the amplitude solutions for ea12 are very low; this is another indication that this antenna is dubious.

Next we use the flux calibrator (whose flux density was set in setjy above) to derive the flux of the other calibrators. Note that the flux table REPLACES the amp.gcal in terms of future application of the calibration to the data, i.e. the flux table contains both the amp.gcal and flux scaling. Unlike the gain calibration steps, this is not an incremental table.

# In CASA
fluxscale(vis='day2_TDEM0003_10s_norx',caltable='amp.gcal',
          fluxtable='flux.cal',reference='7')
Plot of flux corrected amplitude solutions.

It is a good idea to note down for your records the derived flux densities:

Flux density for J0954+1743 in SpW=0 is: 0.23717 +/- 0.0312543 (SNR = 7.58838, N= 38)
Flux density for J0954+1743 in SpW=1 is: 0.240683 +/- 0.0302257 (SNR = 7.96288, N= 38)
Flux density for J1229+0203 in SpW=0 is: 28.2803 +/- 0.245907 (SNR = 115.004, N= 38)
Flux density for J1229+0203 in SpW=1 is: 28.4019 +/- 0.257634 (SNR = 110.241, N= 38)


Next, check that the flux.cal table looks as expected:

# In CASA
plotcal(caltable='flux.cal',xaxis='time',yaxis='amp',
        iteration='antenna',subplot=331)

Note the odd behavior for antennas ea07, ea12, and ea23.

Applycal and Inspect

Now we apply the calibration to each source, according to which tables are appropriate, and which source should be used to do that particular calibration. For the calibrators, all bandpass solutions come from the bandpass calibrator (id=5), and the phase and amplitude calibration comes from their own solutions.

Note: In applycal we set calwt=F. 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. During the VLA era, meaningful weights were available for each visibility. However, at the time of this observation, the EVLA 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=T when the input weights are themselves not meaningful, especially for self-calibration on resolved sources (your flux calibrator and target, for example).

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.

# In CASA
# for the gain/phase calibrator (field '2'):
applycal(vis='day2_TDEM0003_10s_norx',field='2',
        gaintable=['antpos.cal','bandpass.bcal','intphase.gcal','flux.cal'],
        gainfield=['','5','2','2'],
        opacity=[0.0410, 0.0408],gaincurve=T,calwt=F)
# In CASA
# for the bandpass calibrator (field '5'):
applycal(vis='day2_TDEM0003_10s_norx',field='5',
        gaintable=['antpos.cal','bandpass.bcal','intphase.gcal','flux.cal'],
        gainfield=['','5','5','5'],
        opacity=[0.0410, 0.0408],gaincurve=T,calwt=F)
# In CASA
# for the flux calibrator (field '7'):
applycal(vis='day2_TDEM0003_10s_norx',field='7',
        gaintable=['antpos.cal','bandpass.bcal','intphase.gcal','flux.cal'],
        gainfield=['','5','7','7'],
        opacity=[0.0410, 0.0408],gaincurve=T,calwt=F)

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

# In CASA
# for the target source IRC10216 (field '3'):
applycal(vis='day2_TDEM0003_10s_norx',field='3',
        gaintable=['antpos.cal','bandpass.bcal','scanphase.gcal','flux.cal'],
        gainfield=['','5','2','2'],
        opacity=[0.0410, 0.0408],gaincurve=T,calwt=F)

Now inspect the corrected data:

Plot of calibrated amplitudes over time.
# In CASA
plotms(vis='day2_TDEM0003_10s_norx',field='5',ydatacolumn='corrected',
       xaxis='time',yaxis='amp',correlation='RR,LL',
       avgchannel='64',spw='0:4~60',antenna='', coloraxis='antenna1')

This plot shows some data deviating from the average amplitudes. Use methods described above to mark a region for a small number of deviant data points, and click "Locate". You will find that ea12 is responsible.

# In CASA
plotms(vis='day2_TDEM0003_10s_norx',field='2',ydatacolumn='corrected',
       xaxis='time',yaxis='amp',correlation='RR,LL',
       avgchannel='64',spw='0:4~60',antenna='', coloraxis='antenna2')

Here we see some problems, with high points. Mark some regions and locate in 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 you find is that ea07 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 wasn't. 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 you'll see that ea07 is actually pretty noisy throughout.

Plot of antenna ea12 by itself
Plot of antenna ea23 by itself

From the locate we also find that ea12 and ea23 show some high points; to see this, replot baselines using each of them alone:

plotms(vis='day2_TDEM0003_10s_norx',field='2',ydatacolumn='corrected',
       xaxis='time',yaxis='amp',correlation='RR,LL',
       avgchannel='64',spw='0:4~60',antenna='ea12', coloraxis='antenna2')
plotms(vis='day2_TDEM0003_10s_norx',field='2',ydatacolumn='corrected',
       xaxis='time',yaxis='amp',correlation='RR,LL',
       avgchannel='64',spw='0:4~60',antenna='ea23', coloraxis='antenna2')

It may be a 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 initial and second pointing. Recall that these are antennas we became suspicious of while inspecting the calibration solutions.

IRC+10216 as a function of uv-distance.

Now let's see how the target looks. Because the target has resolved structure, its best to look at it as a function of uvdistance. We'll go ahead and exclude the three antennas we already know have problems.

# In 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', coloraxis='antenna2')

The color indicates that the spikes are caused by a single antenna. Use, zoom, mark, and locate to see which one. Also look at spw=1.

The bad antenna turns out to be ea28; to confirm, replot with antenna=!ea28:

# In 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')

Looks much better! To see if it's restricted to a certain time, do:

# In 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')

Baselines with ea28 clearly show issues until about two-thirds of the way through the observation. Plot another distant antenna to compare. We will go ahead and 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.

# In CASA
flagdata(vis='day2_TDEM0003_10s_norx',
         field=['',''],
         spw=['',''],
         antenna=['ea07,ea12,ea28','ea07,ea23'],
         timerange=['','03:21:40~04:10:00'])

Redo Calibration after more Flagging

After flagging, you'll need to repeat the calibration 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.

# In CASA
gaincal(vis='day2_TDEM0003_10s_norx',caltable='bpphase_redo.gcal',
        field='5',spw='0~1:20~40',
        refant='ea02',calmode='p',solint='int',minsnr=2.0,
        gaintable=['antpos.cal'],
        opacity=[0.0410, 0.0408],gaincurve=T)
#
bandpass(vis='day2_TDEM0003_10s_norx',caltable='bandpass_redo.bcal',
        field='5',
        refant='ea02',solint='inf',solnorm=T,
        gaintable=['antpos.cal','bpphase_redo.gcal'],
        opacity=[0.0410, 0.0408],gaincurve=T)
#
gaincal(vis='day2_TDEM0003_10s_norx',caltable='intphase_redo.gcal',
        field='2,5,7',spw='0~1:4~60',
        refant='ea02',calmode='p',solint='int',minsnr=2.0,
        gaintable=['antpos.cal','bandpass_redo.bcal'],
        opacity=[0.0410, 0.0408],gaincurve=T)
#
gaincal(vis='day2_TDEM0003_10s_norx',caltable='scanphase_redo.gcal',
        field='2,5,7',spw='0~1:4~60',
        refant='ea02',calmode='p',solint='inf',minsnr=2.0,
        gaintable=['antpos.cal','bandpass_redo.bcal'],
        opacity=[0.0410, 0.0408],gaincurve=T)
#
gaincal(vis='day2_TDEM0003_10s_norx',caltable='amp_redo.gcal',
        field='2,5,7',spw='0~1:4~60',
        refant='ea02',calmode='ap',solint='inf',minsnr=2.0,
        gaintable=['antpos.cal','bandpass_redo.bcal','intphase_redo.gcal'],
        opacity=[0.0410, 0.0408],gaincurve=T)
#
fluxscale(vis='day2_TDEM0003_10s_norx',caltable='amp_redo.gcal',
          fluxtable='flux_redo.cal',reference='7')
Flux density for J0954+1743 in SpW=0 is: 0.265112 +/- 0.0286296 (SNR = 9.26005, N= 32)
Flux density for J0954+1743 in SpW=1 is: 0.265065 +/- 0.0292281 (SNR = 9.06882, N= 32)
Flux density for J1229+0203 in SpW=0 is: 28.4587 +/- 0.286533 (SNR = 99.3209, N= 32)
Flux density for J1229+0203 in SpW=1 is: 28.335 +/- 0.29411 (SNR = 96.3414, N= 32)

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

# In CASA
applycal(vis='day2_TDEM0003_10s_norx',field='2',
        gaintable=['antpos.cal','bandpass_redo.bcal','intphase_redo.gcal','flux_redo.cal'],
        gainfield=['','5','2','2'],
        opacity=[0.0410, 0.0408],gaincurve=T,calwt=F)
#
applycal(vis='day2_TDEM0003_10s_norx',field='5',
        gaintable=['antpos.cal','bandpass_redo.bcal','intphase_redo.gcal','flux_redo.cal'],
        gainfield=['','5','5','5'],
        opacity=[0.0410, 0.0408],gaincurve=T,calwt=F)
#
applycal(vis='day2_TDEM0003_10s_norx',field='7',
        gaintable=['antpos.cal','bandpass_redo.bcal','intphase_redo.gcal','flux_redo.cal'],
        gainfield=['','5','7','7'],
        opacity=[0.0410, 0.0408],gaincurve=T,calwt=F)
#
applycal(vis='day2_TDEM0003_10s_norx',field='3',
        gaintable=['antpos.cal','bandpass_redo.bcal','scanphase_redo.gcal','flux_redo.cal'],
        gainfield=['','5','2','2'],
        opacity=[0.0410, 0.0408],gaincurve=T,calwt=F)
Gain calibrator after further flagging and recalibration
IRC+10216 after further flagging and recalibration (after selecting colorize by spw).

Now you can inspect the calibrated data again. Except for random scatter things look pretty good.

# In CASA
plotms(vis='day2_TDEM0003_10s_norx',field='2',ydatacolumn='corrected',
       xaxis='time',yaxis='amp',correlation='RR,LL',
       avgchannel='64',spw='0:4~60',antenna='', coloraxis='antenna2')

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.

# In CASA
plotms(vis='day2_TDEM0003_10s_norx',field='3',ydatacolumn='corrected',
       xaxis='uvdist',yaxis='amp',correlation='RR,LL',
       avgchannel='64',spw='0~1:4~60',antenna='', coloraxis='spw')

Split

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

Here, we split off the data for the phase calibrator and the target:

# In CASA
split(vis='day2_TDEM0003_10s_norx',outputvis='J0954',
      field='2')
#
split(vis='day2_TDEM0003_10s_norx',outputvis='IRC10216',
      field='3')

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

# In CASA
clearcal(vis='J0954')
#
clearcal(vis='IRC10216')

This concludes the calibration phase of the data reductions. The tutorial continues with continuum subtraction, imaging, and image analysis in EVLA high frequency Spectral Line tutorial - IRC+10216 part2.

CASAguides

--Crystal Brogan --additions: Juergen Ott, Miriam Krauss, Deb Shepherd

Last checked on CASA Version 3.3.0.