SMA CO Line Data 3.1: Difference between revisions

From CASA Guides
Jump to navigationJump to search
No edit summary
 
(66 intermediate revisions by 5 users not shown)
Line 1: Line 1:
== Overview ==
[http://www.cfa.harvard.edu/sma/general/index.html SMA Intro]


This script describes the data reduction of an SMA six-pointing mosaic of the infrared dark cloud G19.3+0.07 with 12CO(2-1) in the USB from Brogan et al. in prep. Please do not use these data for scientific purposes. The data were filled in Miriad, a Tsys correction was applied, and were then written out in uvfits format. In this tutorial, CASA is used to read the uvfits files, and flag, calibrate, continuum subtract, and image the data. Plotting for datasets with many spw such as SMA data ([[plotxy]]) can be slow. Where possible the new plotter [[plotms]] has been used.
[[File:g19.3CO.png|thumb| Integrated intensity SMA 12CO 2-1 map of G19.3+0.07 overlaid on a Spitzer image]]


This tutorial is based on the SMA CO Line Data script found in the CASA Scripts and Data page [http://casa.nrao.edu/casascripts.shtml here]. You will need the tarball of the uvfits data, found [http://casa.nrao.edu/Data/SMA/spw_fits_files.tar here]. An integrated intensity map of the SMA data overlaid on a Spitzer image can be found [http://casa.nrao.edu/Images/g19_COint.png here].
'''NOTE: This script works with CASA 3.1.0. For an updated script that works with CASA 3.2.0, please go to [[SMA CO Line Data 3.2]].'''
 
This script describes the data reduction of an SMA six-pointing mosaic of 12CO 2-1 emission in the infrared dark cloud G19.3+0.07. The CO line was placed in the USB. These data are from Brogan et al. in prep.; please do not use them for scientific purposes. This tutorial is based on the SMA CO Line Data script found in the [http://casa.nrao.edu/casascripts.shtml CASA Scripts and Data] page.
 
The data were previously filled in [http://www.atnf.csiro.au/computing/software/miriad/ Miriad]. A Tsys correction was applied, and the data were then written out in uvfits format. In this tutorial, CASA is used to read the uvfits files, and flag, calibrate, continuum subtract and image the data. Plotting for datasets with many spw such as SMA data ({{plotxy}}) can be slow. Where possible the new plotter {{plotms}} has been used.
 
You will need the tarball of the uvfits data, which can be downloaded [http://casa.nrao.edu/Data/SMA/spw_fits_files.tar here].
 
<pre style="background-color: #FFFF00;">
This tutorial uses the plotxy task in some places rather than plotms. Plotxy will eventually
be phased out, and this tutorial will be updated. In the meantime, this tutorial will work with
CASA Version 3.1.0.
</pre>


== Creating the measurement set from uvfits data ==
== Creating the measurement set from uvfits data ==
Line 17: Line 29:
</source>
</source>


Read in each uvfits file using the [[importuvfits]] filler (need to put upper range to one more than have, to satisfy ipython peculiarity). We create a list called 'myfiles' of the filenames that can be used in a later call.
Read in each uvfits file using the {{importuvfits}} filler (you need to set the upper range of the for loop to one more than you actually have, to satisfy ipython peculiarity). We create a list of the uvfits filenames in ''myfiles'' that will be used to create a single measurement set (ms).


<source lang="python">
<source lang="python">
# In CASA
os.system('rm -rf spw_ms_files')
os.system('rm -rf spw_ms_files')
os.system('mkdir spw_ms_files')
os.system('mkdir spw_ms_files')
myfiles = []
myfiles = []
for i in range(1,25):
for i in range(1,25):
        msfile = "spw_ms_files/g19_d2usb.spw"+str(i)+".ms"
    msfile = "spw_ms_files/g19_d2usb.spw"+str(i)+".ms"
        importuvfits(fitsfile="spw_fits_files/g19_d2usb.spw"+str(i),
    importuvfits(fitsfile="spw_fits_files/g19_d2usb.spw"+str(i),
                    vis=msfile)
                vis=msfile)
        myfiles.append(msfile)
    myfiles.append(msfile)
</source>
</source>


Create a single ms from the filenames listed in 'myfile' using [[concat]]:  
[[File:listobs.png|thumb| The beginning of the output from listobs. ]]
 
Using {{concat}}, create a single ms from the filenames listed in ''myfiles'':  


<source lang="python">
<source lang="python">
# In CASA
concat(vis=myfiles,concatvis='g19_d2usb.ms',timesort=True)
concat(vis=myfiles,concatvis='g19_d2usb.ms',timesort=True)
</source>
</source>


Finally, initialize the scratch columns using [[clearcal]] and use [[listobs]] to look at the properties in the concatenated file:  
Finally, initialize the scratch columns using {{clearcal}} and use {{listobs}} to look at the contents of the measurement set (which appear in the CASA Log Messages window):  


<source lang="python">
<source lang="python">
# In CASA
clearcal(vis='g19_d2usb.ms')
clearcal(vis='g19_d2usb.ms')
listobs(vis='g19_d2usb.ms',verbose=False)
listobs(vis='g19_d2usb.ms',verbose=False)
</source>
</source>


<pre>
<pre style="background-color: #ffe4b5;">
       0        jupiter      13:01:25.48      -05.18.53.00  J2000
##########################################
       1        1743-038      17:43:58.85      -03.50.04.62  J2000
##### Begin Task: listobs            #####
       2        g19.3a        18:25:58.70      -12.03.57.80  J2000
listobs::::casa
       3        g19.3b        18:25:57.40      -12.04.18.50  J2000
================================================================================
       4        g19.3c        18:25:56.09      -12.04.28.20  J2000
MeasurementSet Name:  /export/data_2/rfriesen/casa_tutorials/test_for_stable/sma_co/sma_test/g19_d2usb.ms      MS Version 2
       5        g19.3d        18:25:54.60      -12.04.23.10  J2000
================================================================================
       6        1908-201      19:11:09.65      -20.06.55.11  J2000
Observer: SmaUser    Project: 
       7        g19.3e        18:25:54.80      -12.04.43.80  J2000
Observation: SMA(8 antennas)
       8        g19.3f        18:25:53.30      -12.04.51.00  J2000
Data records: 571536      Total integration time = 33508 seconds
       9        3c454.3      22:53:57.74      +16.08.53.56  J2000
Observed from  19-Aug-2005/04:23:41.9  to  19-Aug-2005/13:42:09.9 (UTC)
       10        0359+509      03:59:29.75      +50.57.50.15  J2000
Fields: 14
       11        hip19167      04:06:35.04      +50.21.04.54  J2000
      ID  Code Name          RA              Decl          Epoch  SrcId
       12        polaris      02:31:48.70      +89.15.50.72  J2000
       0        jupiter      13:01:25.48      -05.18.53.00  J2000   0
       13        uranus        22:44:40.51      -08.50.23.77  J2000
       1        1743-038      17:43:58.85      -03.50.04.62  J2000   1
       2        g19.3a        18:25:58.70      -12.03.57.80  J2000   2
       3        g19.3b        18:25:57.40      -12.04.18.50  J2000   3
       4        g19.3c        18:25:56.09      -12.04.28.20  J2000   4
       5        g19.3d        18:25:54.60      -12.04.23.10  J2000   5
       6        1908-201      19:11:09.65      -20.06.55.11  J2000   6
       7        g19.3e        18:25:54.80      -12.04.43.80  J2000   7
       8        g19.3f        18:25:53.30      -12.04.51.00  J2000   8
       9        3c454.3      22:53:57.74      +16.08.53.56  J2000   9
       10        0359+509      03:59:29.75      +50.57.50.15  J2000   10
       11        hip19167      04:06:35.04      +50.21.04.54  J2000   11
       12        polaris      02:31:48.70      +89.15.50.72  J2000   12
       13        uranus        22:44:40.51      -08.50.23.77  J2000   13
      (nVis = Total number of time/baseline visibilities per field)
Spectral Windows:(24 unique spectral windows and 1 unique polarization setups)
Spectral Windows:(24 unique spectral windows and 1 unique polarization setups)
       SpwID  #Chans Frame Ch1(MHz)    ChanWid(kHz)TotBW(kHz)  Ref(MHz)    Corrs
       SpwID  #Chans Frame Ch1(MHz)    ChanWid(kHz)TotBW(kHz)  Ref(MHz)    Corrs
Line 88: Line 116:
</pre>
</pre>


Brief description of the purposes of the sources in this project. () gives source id number.
The sources contained in this measurement set are:
 
<pre>
Bandpass calibrator: 3c454.3 (9)
Source Name        Description                                    ID
Gain calibrators:    1743-038 (1) and 1908-201 (6)
====================================================================
Flux calibrators:    3c454.3 (9) and Uranus (13)
3c454.3           Bandpass and flux calibrator                    9
Science target:      g19 (2,3,4,5,7,8) six-pointing mosaic
1743-038           Gain calibrator                                1
1908-201           Gain calibrator                                6
Uranus            Flux calibrator                                13
g19.7              Science target (6-pointing mosaic)    2,3,4,5,7,8
====================================================================
</pre>
Other sources in the ms were not part of this program.
Other sources in the ms were not part of this program.


== Initial inspection and flagging ==
== Initial inspection and flagging ==


Use [[plotms]] to inspect the data. Screengrab here.  
[[File:plotms-noflags.png|thumb| Amplitude vs. time for the calibrator data, pre-flagging.]]
Use {{plotms}} to inspect the data.  


<source lang="python">
<source lang="python">
# In CASA
clearstat()
clearstat()
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='1~9,13',
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='1~9,13',
Line 106: Line 141:
</source>
</source>


There are a few bad data points in this ms. You can use the locate features in [[plotms]] to find the bad data points on field 8 and antenna ids 1-2 and 0-2. Flagging can be done in [[plotms]], or in a script using the [[flagdata]] command. Here we will do all flagging within the script.  
There are a few bad data points in this ms. You can use the locate features in {{plotms}} to find the bad data points on field 8 and antenna ids 1-2 and 0-2. Flagging can be done in {{plotms}}, or in a script using the {{flagdata}} command. Here we will do all flagging within the script.  


<div style="background-color: #dddddd;">
<div style="background-color: #dddddd;">
'''Tip:''' locate gives antenna ids, but [[flagdata]] looks first at the name column. Since SMA antenna names are numbers, this can be confusing. The corresponding ''names'' are antennas 2-3 and 1-3, which you can check with [[listobs]].
'''Tip:''' locate gives antenna ids, but {{flagdata}} looks first at the name column. Since SMA antenna names are numbers, this can be confusing. The corresponding ''names'' are antennas 2-3 and 1-3, which you can check with {{listobs}}.
</div>
</div>
[[File:plotms-flags.png|thumb| Amplitude vs. time for the calibrator data, after flagging bad data points.]]


<source lang="python">
<source lang="python">
# In CASA
flagdata(vis='g19_d2usb.ms',mode='manualflag',field='8',spw='20~23',
flagdata(vis='g19_d2usb.ms',mode='manualflag',field='8',spw='20~23',
timerange='07:52:45~07:52:47',antenna='2&3;1&3')
timerange='07:52:45~07:52:47',antenna='2&3;1&3')
</source>
</source>


Plotting again shows that the bad data are gone: (screen grab)
Plotting again shows that the bad data are gone:


<source lang="python">
<source lang="python">
# In CASA
clearstat()
clearstat()
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='8',
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='8',
Line 126: Line 165:


<div style="background-color: #dddddd;">
<div style="background-color: #dddddd;">
'''Tip:''' to confirm the data were flagged in [[plotms]], you need to change selection in plot (for example change field and then change back), and then change back to original to clear the cache and force a re-read of the data.  
'''Tip:''' to confirm the data were flagged in {{plotms}}, you need to change selection in plot (for example change field and then change back), and then change back to original to clear the cache and force a re-read of the data.  
</div>
</div>


Make an inital amplitude vs. channel plot of the bandpass calibrator 3c454.3, averaging all the data together in time. Here, we are using [[plotxy]] because we would like to have multiple subplots (not yet functional in [[plotms]]). We set the [[plotxy]] parameter '''crossscan=T''' because there are two scans on 3c454.3. Note the rolloff of 24 individual spw. Try zooming to see better. (Screengrab here)
[[File:ampVfreq_colorspw.png|thumb| Amplitude vs. frequency for the bandpass calibrator 3c454.3.]]


<div style="background-color: #dddddd;">
Make an inital amplitude vs. frequency plot of the bandpass calibrator 3c454.3, averaging the data together in time. Here, we colorize the points by spw in the {{plotms}} call using the '''coloraxis''' parameter. You can also do this under the 'Display' tab in the plotms window.  
'''Tip:''' after you run [[plotxy]], run clearstat() before trying to run [[plotms]] or you will get a table lock.
</div>


<source lang="python">
<source lang="python">
# In CASA
clearstat()
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp",datacolumn="data",
plotms(vis="g19_d2usb.ms",xaxis="freq",yaxis="amp",field="9",
      iteration="baseline",
       averagedata=True,avgscan=True,avgtime='500',avgbaseline=True,coloraxis="spw")  
      antenna="",spw="",field="9",
       averagemode="vector",width="1",timebin="all",
      crossscans=T,crossbls=False,stackspw=F,
      subplot=331,plotcolor="darkcyan",
      overplot=False,showflags=False,interactive=True,figfile="")
</source>
</source>


As you can see in this plot, SMA spw overlap by a fair amount. Here we will flag the first seven channels on either side of each spw to aid calibration:
As you can see, SMA spw overlap by a fair amount. Here we will flag the first seven channels on either side of each spw to aid calibration:


<source lang="python">
<source lang="python">
# In CASA
default('flagdata')
default('flagdata')
flagdata(vis='g19_d2usb.ms', mode='manualflag',spw='0~23:0~6;121~127')
flagdata(vis='g19_d2usb.ms', mode='manualflag',spw='0~23:0~6;121~127')
</source>
</source>


On many baselines the first integration of every scan is low. These can also be flagged with [[flagdata]]. The integration time of these data is 30.9s. Note that more recent SMA data rarely show this problem.
[[File:plotxy-ants.png|thumb| Antenna positions for these data.]]


On many baselines the first integration of every scan is low. This can be seen with {{plotms}}.
<source lang="python">
<source lang="python">
# In CASA
clearstat()
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",
      field='1,6,9',avgchannel='3072',avgspw=T,coloraxis="field")
</source>
Here, we have colored the points by field to see that the low amplitudes are found in the first integrations of many scans.  These can also be flagged with {{flagdata}}. The integration time of these data is 30.9s. Note that more recent SMA data rarely show this problem. 
[[File:Plotms_amp_vs_time_before_quack.png|thumb| Amplitude vs. Time before flagging first integrations.]]
<source lang="python">
# In CASA
flagdata(vis="g19_d2usb.ms",mode="quack",field='',
flagdata(vis="g19_d2usb.ms",mode="quack",field='',
spw="0~23",quackinterval=40.0)
spw="0~23",quackinterval=40.0)
</source>
</source>


[[plotxy] Plotxy] can be used to plot the antenna positions. For these observations, antenna 8 was in the barn for repair.
We can use {{plotxy}} to plot the antenna positions. For these observations, antenna 8 was in the barn for repair.


<source lang="python">
<source lang="python">
# In CASA
clearstat()
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="x")
plotxy(vis="g19_d2usb.ms",xaxis="x")
</source>
</source>
<div style="background-color: #dddddd;">
'''Tip:''' after you run {{plotxy}}, run clearstat() before trying to run {{plotms}} or you will get a table lock.
</div>


== Bandpass Calibration ==
== Bandpass Calibration ==


Next plot phase vs. channel. Note the large discrepancy in phase for spw=4-7 on some baselines, this will be important later.
[[File:plotxy-pha_vs_chan.png|thumb| Phase vs. channel for the bandpass calibrator 3c454.3.]]
 
Next plot phase vs. channel. Here, we use plotms to iterate over the SMA baselines. While multiple subplots are not yet operational within plotms, you can view the individual phase vs. channel plots using the arrow buttons at the bottom right of the {{plotms}} window. By setting the '''xaxis''' parameter to 'frequency', we can see how the spectral windows overlap. Note the large discrepancy in phase for spw=4-7 on some baselines - this will be important later.
 
<source lang="python">
clearstat()
plotms(vis="g19_d2usb.ms",xaxis="frequency",yaxis="phase",
      ydatacolumn="data",field="9",iteraxis="baseline",
      averagedata=True,avgtime="99999",colorize=True,coloraxis="spw")
</source>
 
We can also use the slower {{plotxy}} to show multiple subplots in one panel. We set the {{plotxy}} parameter '''crossscan=T''' because there are two scans on 3c454.3. 


<source lang="python">
<source lang="python">
# In CASA
clearstat()
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="phase",datacolumn="data",
plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="phase",datacolumn="data",
Line 181: Line 244:
       overplot=False,showflags=False,interactive=True,figfile="")
       overplot=False,showflags=False,interactive=True,figfile="")
</source>
</source>
[[File:plotms-amp_vs_time.png|thumb| Amplitude vs. time for the bandpass calibrator 3c454.3.]]


Now plot amplitude and phase as a function of time.
Now plot amplitude and phase as a function of time.


<source lang="python">
<source lang="python">
# In CASA
clearstat()
clearstat()
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='9',
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='9',
Line 190: Line 256:
</source>
</source>


The gap in the 3C454.3 observation is due to a pause to determine pointing solutions.  The flags below catch times just before and after when the amplitudes and phase are off. This flagging could also have been accomplished interactively in [[plotms]].
The gap in the 3C454.3 observation is due to a pause to determine pointing solutions.  The flags below catch times just before and after when the amplitudes and phase are off. This flagging could also have been accomplished interactively in {{plotms}}.


<source lang="python">
<source lang="python">
# In CASA
default('flagdata')
default('flagdata')
flagdata(vis="g19_d2usb.ms",field="9",timerange='12:16:58~12:19:00')
flagdata(vis="g19_d2usb.ms",field="9",timerange='12:16:58~12:19:00')
Line 198: Line 265:
</source>
</source>


Here because [[plotms]] cannot do multipanel iteration yet we use the slower [[plotxy]].
[[File:plotxy-pha_vs_time.png|thumb| Phase vs. time for the bandpass calibrator 3c454.3.]]
 
Here because {{plotms}} cannot do multipanel iteration yet we use the slower {{plotxy}}.


<source lang="python">
<source lang="python">
# In CASA
clearstat()
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="phase",datacolumn="data",
plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="phase",datacolumn="data",
Line 210: Line 280:
</source>
</source>


Start by determining an initial phase only solution for 3c454.3 using [[gaincal]] to take out major phase variations with time. For maximum sensitivity, we average all the spw together, with the exception of spw 4-7 which show large phase offsets from the other spw in phase vs. channel.   
[[File:plotcal-pha_vs_time.png|thumb| Phase solutions, all spw averaged together (minus spw 4-7).]]
 
Start by determining an initial phase only solution for 3c454.3 using {{gaincal}} to take out major phase variations with time. For maximum sensitivity, we average all the spw together, with the exception of spw 4-7 which show large phase offsets from the other spw in phase vs. channel.   


<source lang="python">
<source lang="python">
# In CASA
gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.pcal',
gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.pcal',
field='9', spw='0~3,8~23', gaintype='G', minsnr=2.0,
field='9', spw='0~3,8~23', gaintype='G', minsnr=2.0,
refant='3', calmode='p',solint='int', combine='spw')
refant='3', calmode='p',solint='int', combine='spw')
clearstat()
clearstat()
plotcal(caltable="g19_d2usb.ms.pcal",xaxis="time",yaxis="phase",
plotcal(caltable="g19_d2usb.ms.pcal",xaxis="time",yaxis="phase",
Line 224: Line 296:
</source>
</source>


Next we do a bandpass calibration. There are several possible solution types. A regular 'B' type solution can work well if you have lots of S/N on your bandpass calibrator. You may wish to make this table and compare to 'BPOLY' solution further down. '''elaborate'''
Next we do a bandpass calibration. There are several possible solution types. A regular 'B' type solution can work well if you have lots of S/N on your bandpass calibrator. You may wish to make this table and compare to 'BPOLY' solution further down.  


<source lang="python">
<source lang="python">
# In CASA
bandpass(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.bpcal',
bandpass(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.bpcal',
  field='9', spw='0~23', bandtype='B', solint='inf', combine='scan',  
  field='9', spw='0~23', bandtype='B', solint='inf', combine='scan',  
  refant='3', gaincurve=False, opacity=0.0,
  refant='3', gaincurve=False, opacity=0.0,
  gaintable='g19_d2usb.ms.pcal',spwmap=[0])  
  gaintable='g19_d2usb.ms.pcal',spwmap=[0])  
plotcal(caltable="g19_d2usb.ms.bpcal",xaxis="freq",yaxis="amp",
plotcal(caltable="g19_d2usb.ms.bpcal",xaxis="freq",yaxis="amp",
field="9",antenna="",spw="",timerange="",subplot=311,
field="9",antenna="",spw="",timerange="",subplot=311,
Line 237: Line 309:
plotsymbol="o",plotcolor="blue",showgui=True,figfile="")  
plotsymbol="o",plotcolor="blue",showgui=True,figfile="")  
</source>
</source>
</div>


A polynomial bandpass solution is more sensitive than simply fitting a polynomial to the channel based solutions 'B' after the fact. For the mm/submm, where the S/N is often low on the bandpass calibrator, 'BPOLY' is often the best choice, but every dataset is different. You may especially need to play with the '''degamp''' and '''degphase''' parameters in [[bandpass]]. The parameter '''combine='scan'''' is needed to average together the two 3c454.3 scans. The '''spwmap=[0]''' is telling it to apply the phase solution that was found from the average of all spw to each spw before determining the bandpass solution.  
[[File:plotcal-amp_vs_freq.png|thumb| Bandpass solution for the bandpass calibrator 3c454.3.]]
 
A polynomial bandpass solution is more sensitive than simply fitting a polynomial to the channel based solutions 'B' after the fact. For the mm/submm, where the S/N is often low on the bandpass calibrator, 'BPOLY' is often the best choice, but every dataset is different. You may especially need to play with the '''degamp''' and '''degphase''' parameters in {{bandpass}}. The parameter '''combine='scan'''' is needed to average together the two 3c454.3 scans. The '''spwmap=[0]''' is telling it to apply the phase solution that was found from the average of all spw to each spw before determining the bandpass solution.  


<source lang="python">
<source lang="python">
# In CASA
bandpass(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.bpoly7',
bandpass(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.bpoly7',
field='9', spw='0~23', bandtype='BPOLY',  
field='9', spw='0~23', bandtype='BPOLY',  
Line 247: Line 321:
refant='3', degamp=7, degphase=7,
refant='3', degamp=7, degphase=7,
gaintable='g19_d2usb.ms.pcal',spwmap=[0])
gaintable='g19_d2usb.ms.pcal',spwmap=[0])
clearstat()
clearstat()
plotcal(caltable="g19_d2usb.ms.bpoly7",xaxis="freq",yaxis="amp",
plotcal(caltable="g19_d2usb.ms.bpoly7",xaxis="freq",yaxis="amp",
Line 257: Line 330:
Unfortunately, this kind of plot with many spw and colors takes a long time. Be patient for now, as it's being worked on. You can safely ignore where the solution gets large at edges of spw. This is showing where we flagged channels - this is also being worked on. Zoom in if necessary to see that inner portions are nice and flat.
Unfortunately, this kind of plot with many spw and colors takes a long time. Be patient for now, as it's being worked on. You can safely ignore where the solution gets large at edges of spw. This is showing where we flagged channels - this is also being worked on. Zoom in if necessary to see that inner portions are nice and flat.


The next step is to apply the bandpass solution while re-calculating with [[gaincal]] the phase only solution for 3c454.3 on a short time interval (here, we use the integration time), but now including all spw (including 4~7 since antenna based solutions will have taken out phase vs.spw differences).
[[File:plotcal-pha_vs_time-2.png|thumb| Phase solutions for the bandpass calibrator 3c454.3, including all spw.]]
 
The next step is to apply the bandpass solution while re-calculating with {{gaincal}} the phase only solution for 3c454.3 on a short time interval (here, we use the integration time), but now including all spw (including 4~7 since antenna based solutions will have taken out phase vs.spw differences).


<source lang="python">
<source lang="python">
# In CASA
gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.pcal2',
gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.pcal2',
field='9', spw='0~23', gaintype='G', minsnr=2.0,
field='9', spw='0~23', gaintype='G', minsnr=2.0,
refant='3', calmode='p',solint='int', combine='spw',
refant='3', calmode='p',solint='int', combine='spw',
gaintable=['g19_d2usb.ms.bpoly7'])
gaintable=['g19_d2usb.ms.bpoly7'])
clearstat()
clearstat()
plotcal(caltable="g19_d2usb.ms.pcal2",xaxis="time",yaxis="phase",
plotcal(caltable="g19_d2usb.ms.pcal2",xaxis="time",yaxis="phase",
Line 272: Line 347:
</source>
</source>


Next we apply the antenna based bandpass, and phase only solution to solve for amplitude solutions over a longer '''solint''' to increase the S/N. The '''spwmap=[[0],[]]''' tells [[gaincal]] that the pcal solution is independent of spw while the bpoly7 solution is per spw. '''spwmap''' order must be in the same order as the tables in '''gaintable'''.
[[File:plotcal-amp_vs_time-new.png|thumb| Amplitude solutions for the bandpass calibrator 3c454.3, including all spw.]]
 
Next we apply the antenna based bandpass, and phase only solution to solve for amplitude solutions over a longer '''solint''' to increase the S/N. The '''spwmap=[[0],[]]''' tells {{gaincal}} that the pcal solution is independent of spw while the bpoly7 solution is per spw. '''spwmap''' order must be in the same order as the tables in '''gaintable'''.


<source lang="python">
<source lang="python">
# In CASA
gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.apcal',
gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.apcal',
field='9', spw='0~23', gaintype='G', minsnr=2.0,
field='9', spw='0~23', gaintype='G', minsnr=2.0,
Line 280: Line 358:
gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.bpoly7'],
gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.bpoly7'],
spwmap=[[0],[]])
spwmap=[[0],[]])
clearstat()
clearstat()
plotcal(caltable="g19_d2usb.ms.apcal",xaxis="time",yaxis="amp",
plotcal(caltable="g19_d2usb.ms.apcal",xaxis="time",yaxis="amp",
Line 288: Line 365:
</source>
</source>


Although in previous steps solutions were applied "on-the-fly", the data were not actually changed. Running [[applycal]] will fill the corrected data column in the ms. Again be careful with specification of '''gaintable''', and keeping the same order in '''spwmap''' and '''gainfield'''. This [[applycal]] is not necessary for later steps, but allows plotting of data calibration to this point.
Although in previous steps solutions were applied "on-the-fly", the data were not actually changed. Running {{applycal}} will fill the corrected data column in the ms. Again be careful with specification of '''gaintable''', and keeping the same order in '''spwmap''' and '''gainfield'''. This {{applycal}} is not necessary for later steps, but allows plotting of data calibration to this point.
 
[[File:plotxy-amp_vs_chan.png|thumb| Gain-corrected 3c454.3 data. Notice the few residual baseline-based problems.]]


<source lang="python">
<source lang="python">
# In CASA
applycal(vis='g19_d2usb.ms',spw='0~23', field='9',
applycal(vis='g19_d2usb.ms',spw='0~23', field='9',
gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.apcal',
gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.apcal',
'g19_d2usb.ms.bpoly7'],
'g19_d2usb.ms.bpoly7'],
spwmap=[[0],[0],[]],gainfield=['9','9','9'])
spwmap=[[0],[0],[]],gainfield=['9','9','9'])
clearstat()
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp",
plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp",
Line 307: Line 386:
</source>
</source>


As this [[plotxy]] shows, after antenna based bandpass, SMA data often show residual baseline based problems as a function of spw (i.e. a few abnormally low amplitudes for a few spw on a few baselines). To fix this problem we need to do a baseline based amplitude only solution per spw, but we'll need to remove an overall normalization factor first by combining all spw (this option will be incorporated into [[blcal]] one day so that only one run is necessary).
As this {{plotxy}} shows, after antenna based bandpass, SMA data often show residual baseline based problems as a function of spw (i.e. a few abnormally low amplitudes for a few spw on a few baselines). To fix this problem we need to do a baseline based amplitude only solution per spw, but we'll need to remove an overall normalization factor first by combining all spw (this option will be incorporated into {{blcal}} one day so that only one run is necessary).  Note the significant change on baseline 1-6 in the following call to {{plotcal}}.


<source lang="python">
<source lang="python">
# In CASA
blcal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.blcal1',
blcal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.blcal1',
field='9', spw='0~23',solint='inf', combine='spw,scan',  
field='9', spw='0~23',solint='inf', combine='spw,scan',  
Line 315: Line 395:
'g19_d2usb.ms.bpoly7'],  
'g19_d2usb.ms.bpoly7'],  
calmode='a',spwmap=[[0],[0],[]])
calmode='a',spwmap=[[0],[0],[]])
blcal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.blcal2',
blcal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.blcal2',
field='9', spw='0~23',solint='inf', combine='scan',  
field='9', spw='0~23',solint='inf', combine='scan',  
Line 321: Line 400:
'g19_d2usb.ms.blcal1','g19_d2usb.ms.bpoly7',],  
'g19_d2usb.ms.blcal1','g19_d2usb.ms.bpoly7',],  
calmode='a',spwmap=[[0],[0],[0],[]])
calmode='a',spwmap=[[0],[0],[0],[]])
clearstat()
clearstat()
plotcal(caltable="g19_d2usb.ms.blcal1",xaxis="freq",yaxis="amp",
plotcal(caltable="g19_d2usb.ms.blcal2",xaxis="freq",yaxis="amp",
field="9",antenna="",spw="",timerange="",subplot=511,
field="9",antenna="",spw="",timerange="",subplot=511,
overplot=False,clearpanel="Auto",iteration="antenna",
overplot=False,clearpanel="Auto",iteration="antenna",
Line 329: Line 407:
</source>
</source>


Now we can do a new [[applycal]] including ''only'' the second blcal2 solution. Note that [[applycal]] overwrites the corrected column, so there is no need to [[clearcal]]. Like the antenna based bandpass the [[blcal]] solution is independent for each spw, so we set '''spwmap = []'''.
Now we can do a new {{applycal}} including ''only'' the second blcal2 solution. Note that {{applycal}} overwrites the corrected column, so there is no need to {{clearcal}}. Like the antenna based bandpass the {{blcal}} solution is independent for each spw, so we set '''spwmap = []'''.
 
[[File:plotxy-amp_vs_chan-2.png|thumb| Gain-corrected 3c454.3 data, with residual problems fixed.]]


<source lang="python">
<source lang="python">
# In CASA
applycal(vis='g19_d2usb.ms',spw='0~23', field='9',
applycal(vis='g19_d2usb.ms',spw='0~23', field='9',
gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.apcal',
gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.apcal',
Line 341: Line 422:


<source lang="python">
<source lang="python">
# In CASA
clearstat()
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp",
plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp",
Line 354: Line 436:
== Inspection of gain calibrators and setting the flux scale ==
== Inspection of gain calibrators and setting the flux scale ==


Apply antenna and baseline based bandpass solutions to the gain calibrators with [[applycal]]:
Apply antenna and baseline based bandpass solutions to the gain calibrators with {{applycal}}:


<source lang="python">
<source lang="python">
# In CASA
clearstat()
clearstat()
applycal(vis='g19_d2usb.ms',spw='0~23', field='1,6',
applycal(vis='g19_d2usb.ms',spw='0~23', field='1,6',
Line 362: Line 445:
spwmap=[[],[]],gainfield=['9','9'])
spwmap=[[],[]],gainfield=['9','9'])
</source>
</source>
[[File:plotxy-amp_vs_chan-3.png|thumb| Corrected amplitude vs. channel for 1743-038.]]


Have a look at the corrected data for both calibrators with channel vs. amplitude:
Have a look at the corrected data for both calibrators with channel vs. amplitude:


<source lang="python">
<source lang="python">
# In CASA
# First field 1:
# First field 1:
clearstat()
clearstat()
Line 384: Line 470:
       overplot=False,showflags=False,interactive=True,figfile="")
       overplot=False,showflags=False,interactive=True,figfile="")
</source>
</source>
[[File:plotxy-amp_vs_time.png|thumb| Corrected amplitude vs. time for 1743-038.]]
[[File:plotxy-pha_vs_time-corr1.png|thumb| Corrected phase vs. time for 1743-038.]]


Have a look at the corrected data for both calibrators in phase and amplitude vs. time. For the moment, unfortunately, '''multicolor='field'''' doesn't work, so this is easiest to do separately.
Have a look at the corrected data for both calibrators in phase and amplitude vs. time. For the moment, unfortunately, '''multicolor='field'''' doesn't work, so this is easiest to do separately.


<source lang="python">
<source lang="python">
# In CASA
# First, amplitude vs. time for field 1:  
# First, amplitude vs. time for field 1:  
clearstat()
clearstat()
Line 425: Line 516:
</source>
</source>


At the time of these observations, the calibrator J1911 had the most recent and relatively stable monitor data (SMA website), so we will use it for absolute flux calibration. In future tools for using resolved planets will be implemented. ** How to set up nice table format?
At the time of these observations, the calibrator J1911 had the most recent and relatively stable monitor data (SMA website), so we will use it for absolute flux calibration. In future tools for using resolved planets will be implemented.  
 
<pre>
1911-201
1911-201
1mm  17 Jul 2005 11:38  SMA      274.5    1.88 +/-  0.13    mgurwell   
1mm  17 Jul 2005 11:38  SMA      274.5    1.88 +/-  0.13    mgurwell   
1mm  05 Sep 2005 08:10  SMA      221.4    1.99 +/-  0.12    mgurwell   
1mm  05 Sep 2005 08:10  SMA      221.4    1.99 +/-  0.12    mgurwell   
</pre>


<source lang="python">
<source lang="python">
# In CASA
setjy(vis='g19_d2usb.ms', field='6', spw='0~23', fluxdensity=[2.0,0.,0.,0.])
setjy(vis='g19_d2usb.ms', field='6', spw='0~23', fluxdensity=[2.0,0.,0.,0.])
</source>
</source>
Line 439: Line 534:


<source lang="python">
<source lang="python">
# In CASA
gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.allpcal',
gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.allpcal',
field='1,6,9', spw='0~23', gaintype='G', minsnr=2.0,
field='1,6,9', spw='0~23', gaintype='G', minsnr=2.0,
Line 449: Line 545:


<source lang="python">
<source lang="python">
# In CASA
gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.allpcalscan',
gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.allpcalscan',
field='1,6,9', spw='0~23', gaintype='G', minsnr=2.0,
field='1,6,9', spw='0~23', gaintype='G', minsnr=2.0,
Line 455: Line 552:
spwmap=[[],[]])
spwmap=[[],[]])
</source>
</source>
[[File:plotcal-pha_vs_time-3.png|thumb| Phase solutions for 1743-038 and 1908-201.]]


Plot the phase solutions of the phase calibrators. Can repeat for both short and long solution intervals.
Plot the phase solutions of the phase calibrators. Can repeat for both short and long solution intervals.


<source lang="python">
<source lang="python">
# In CASA
clearstat()
clearstat()
plotcal(caltable="g19_d2usb.ms.allpcalscan",xaxis="time",yaxis="phase",
plotcal(caltable="g19_d2usb.ms.allpcalscan",xaxis="time",yaxis="phase",
Line 469: Line 569:


<source lang="python">
<source lang="python">
# In CASA
gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.allapcal',
gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.allapcal',
field='1,6,9', spw='0~23', gaintype='G', minsnr=2.0,
field='1,6,9', spw='0~23', gaintype='G', minsnr=2.0,
Line 475: Line 576:
'g19_d2usb.ms.blcal2'],spwmap=[[0],[],[]])
'g19_d2usb.ms.blcal2'],spwmap=[[0],[],[]])
</source>
</source>
[[File:plotcal-amp_vs_time-2.png|thumb| Amplitude solutions for 1743-038 and 1908-201.]]


Plot the amplitude solutions:
Plot the amplitude solutions:


<source lang="python">
<source lang="python">
# In CASA
clearstat()
clearstat()
plotcal(caltable="g19_d2usb.ms.allapcal",xaxis="time",yaxis="amp",
plotcal(caltable="g19_d2usb.ms.allapcal",xaxis="time",yaxis="amp",
Line 486: Line 590:
</source>
</source>


Derive the absolute flux calibration using [[fluxscale]] (based on [[setjy]] above):
Derive the absolute flux calibration using {{fluxscale}} (based on {{setjy}} above):


<source lang="python">
<source lang="python">
# In CASA
fluxscale(vis='g19_d2usb.ms',caltable='g19_d2usb.ms.allapcal',
fluxscale(vis='g19_d2usb.ms',caltable='g19_d2usb.ms.allapcal',
fluxtable='g19_d2usb.ms.fluxcal',reference='6')
fluxtable='g19_d2usb.ms.fluxcal',reference='6')
Line 496: Line 601:


The final application of calibration for each source is done independently below. Flux densities from SMA flux monitoring are provided as checks of the absolute flux calibration.
The final application of calibration for each source is done independently below. Flux densities from SMA flux monitoring are provided as checks of the absolute flux calibration.
[[File:plotms-amp_vs_time-2.png|thumb| Corrected amplitude vs. time for 3C454.3.]]


<source lang="python">
<source lang="python">
# In CASA
applycal(vis='g19_d2usb.ms',spw='0~23', field='9',
applycal(vis='g19_d2usb.ms',spw='0~23', field='9',
gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.fluxcal',
gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.fluxcal',
'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]],
'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]],
gainfield=['9','9','9','9'])
gainfield=['9','9','9','9'])
clearstat()
clearstat()
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='9',
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='9',
Line 512: Line 619:
1mm  19 Aug 2005 16:01  SMA      225.3  27.36 +/-  1.37    mgurwell  
1mm  19 Aug 2005 16:01  SMA      225.3  27.36 +/-  1.37    mgurwell  
</pre>
</pre>
[[File:plotms-amp_vs_time-3.png|thumb| Corrected amplitude vs. time for 1743-038.]]


<source lang="python">
<source lang="python">
# In CASA
applycal(vis='g19_d2usb.ms',spw='0~23', field='1',
applycal(vis='g19_d2usb.ms',spw='0~23', field='1',
gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.fluxcal',
gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.fluxcal',
'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]],
'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]],
gainfield=['1','1','9','9'])
gainfield=['1','1','9','9'])
clearstat()
clearstat()
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='1',
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='1',
Line 528: Line 637:
1mm  26 Aug 2005 06:05  SMA      225.6    1.71 +/-  0.09    mgurwell  
1mm  26 Aug 2005 06:05  SMA      225.6    1.71 +/-  0.09    mgurwell  
</pre>
</pre>
[[File:plotms-amp_vs_time-4.png|thumb| Corrected amplitude vs. time for 1908-201.]]


<source lang="python">
<source lang="python">
# In CASA
applycal(vis='g19_d2usb.ms',spw='0~23', field='6',
applycal(vis='g19_d2usb.ms',spw='0~23', field='6',
gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.fluxcal',
gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.fluxcal',
'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]],
'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]],
gainfield=['6','6','9','9'])
gainfield=['6','6','9','9'])
clearstat()
clearstat()
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='6',
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='6',
Line 546: Line 657:
</pre>
</pre>


Since Uranus is resolved we apply the solutions from 3C454.3 to achieve at lest some reduction in instrumental gain variations.
[[File:plotms-amp_vs_uvdist-uranus.png|thumb| Corrected amplitude vs. uvdistance for Uranus.]]
 
Since Uranus is resolved we apply the solutions from 3C454.3 to achieve at least some reduction in instrumental gain variations.


<source lang="python">
<source lang="python">
# In CASA
applycal(vis='g19_d2usb.ms',spw='0~23', field='13',
applycal(vis='g19_d2usb.ms',spw='0~23', field='13',
gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.fluxcal',
gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.fluxcal',
'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]],
'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]],
gainfield=['9','9','9','9'])
gainfield=['9','9','9','9'])
clearstat()
clearstat()
plotms(vis="g19_d2usb.ms",xaxis="uvdist",yaxis="amp",field='13',
plotms(vis="g19_d2usb.ms",xaxis="uvdist",yaxis="amp",field='13',
Line 562: Line 675:


<source lang="python">
<source lang="python">
# In CASA
applycal(vis='g19_d2usb.ms',spw='0~23', field='2,3,4,5,7,8',
applycal(vis='g19_d2usb.ms',spw='0~23', field='2,3,4,5,7,8',
gaintable=['g19_d2usb.ms.allpcalscan','g19_d2usb.ms.fluxcal',
gaintable=['g19_d2usb.ms.allpcalscan','g19_d2usb.ms.fluxcal',
Line 571: Line 685:


<source lang="python">
<source lang="python">
# In CASA
split(vis="g19_d2usb.ms",outputvis='g19_d2usb_targets.ms',
split(vis="g19_d2usb.ms",outputvis='g19_d2usb_targets.ms',
       field="2,3,4,5,7,8")
       field="2,3,4,5,7,8")
</source>
To reinitialize the scratch columns for use by later tasks, we need to run {{clearcal}} for the new dataset. Then run {{listobs}} to check the target data names.
<source lang="python">
# In CASA
clearcal(vis='g19_d2usb_targets.ms')


listobs(vis="g19_d2usb_targets.ms",verbose=F)
listobs(vis="g19_d2usb_targets.ms",verbose=F)
</source>
</source>


<div style="background-color: #dddddd;">
<pre>
<pre>
ID  Code Name          Right Ascension  Declination  Epoch
ID  Code Name          Right Ascension  Declination  Epoch
Line 587: Line 708:
5        g19.3f        18:25:53.30      -12.04.51.00  J2000
5        g19.3f        18:25:53.30      -12.04.51.00  J2000
</pre>
</pre>
<div style="background-color: #dddddd;">
'''Tip:''' After splitting, the "corrected" data column in the original vis file will become the "data" and "corrected" columns and the field IDs are renumbered.
</div>
</div>


Note that after splitting, the "corrected" data column in the original vis file will become the "data" and "corrected" columns and the field IDs are renumbered.
[[File:plotms-amp_vs_freq.png|thumb| Amplitude vs. frequency for the target data.]]


<source lang="python">
<source lang="python">
# In CASA
clearstat()
clearstat()
plotms(vis='g19_d2usb_targets.ms',xaxis="frequency",yaxis="amp",field='0',
plotms(vis='g19_d2usb_targets.ms',xaxis="frequency",yaxis="amp",field='0',
Line 601: Line 726:
</source>
</source>


Data at times > 11:04:00 show lots of scatter: the elevation of these data are very low and are worth flagging. One can [[plotxy]] with '''xaxis="time"''', '''yaxis="elevation"''' to see this.
Data at times > 11:04:00 show lots of scatter: the elevation of these data are very low and are worth flagging. One can {{plotms}} with '''xaxis="time"''', '''yaxis="elevation"''' to see this.
For now to replot you must change the selection:
<source lang="python">
# In CASA
clearstat()
plotms(vis='g19_d2usb_targets.ms',xaxis="time",yaxis="elevation",field='')
</source>
[[File:el_time_ppt.png|thumb| Elevation vs. Time for the target data.]]


<source lang="python">
<source lang="python">
# In CASA
flagdata(vis="g19_d2usb_targets.ms",field="",timerange='>11:04:00')
flagdata(vis="g19_d2usb_targets.ms",field="",timerange='>11:04:00')
</source>
</source>


Note that because the Doppler tracked frequency was 13CO (+ 10GHz in the USB), the velocity diplayed in in [[plotxy]] will not be correct unless you set the '''restfreq''' to the line of interest. [[plotms] Plotms] cannot do this yet.
[[File:plotxy-amp_vs_vel.png|thumb| Amplitude vs. velocity for the target data.]]
 
Note that because the Doppler tracked frequency was 13CO (+ 10GHz in the USB), the velocity diplayed in in {{plotxy}} will not be correct unless you set the '''restfreq''' to the line of interest. {{plotms}} cannot do this yet.


CO(2-1) restfreq=230.53797 GHz
CO(2-1) restfreq=230.53797 GHz


<source lang="python">
<source lang="python">
# In CASA
clearstat()
clearstat()
plotxy(vis="g19_d2usb_targets.ms",xaxis="velocity",yaxis="amp",
plotxy(vis="g19_d2usb_targets.ms",xaxis="velocity",yaxis="amp",
Line 628: Line 762:
== Continuum imaging and subtraction ==
== Continuum imaging and subtraction ==


Use field 0 to pick line free spw. You need to make a clean mask. The strongest source is in the upper NE of the mosaic and there is a much weaker source to the SW corner.
Next we want to make a continuum-free line dataset. We use field 0 to pick the line free spw. You need to make a clean mask. The strongest source is in the upper NE of the mosaic and there is a much weaker source to the SW corner.


<source lang="python">
<source lang="python">
# In CASA
clearstat()
clearstat()
plotxy(vis="g19_d2usb_targets.ms",xaxis="channel",yaxis="amp",
plotxy(vis="g19_d2usb_targets.ms",xaxis="channel",yaxis="amp",
Line 642: Line 777:
The line is confined to spws 13~15.
The line is confined to spws 13~15.


Copy original ms in case something goes wrong, currently [[uvcontsub2]] will change the input model and corrected columns so its best to make a back up. This behavior will change in future.
Copy the original ms in case something goes wrong. Currently {{uvcontsub2}} will change the input model and corrected columns, so it's best to make a back up. This behavior will change in future.


<source lang="python">
<source lang="python">
# In CASA
os.system('cp -r g19_d2usb_targets.ms g19_d2usb_targets_line.ms')
os.system('cp -r g19_d2usb_targets.ms g19_d2usb_targets_line.ms')
</source>
</source>


Now use new task [[uvcontsub2]] to subtract the continuum. Unlike [[uvcontsub]], this task will subtract the continuum estimate for spw that are ''not'' in the fitspw list. In future releases [[uvcontsub2]] will likely replace [[uvcontsub]].
Now use new task {{uvcontsub2}} to subtract the continuum. Unlike {{uvcontsub}}, this task will subtract the continuum estimate for spw that are ''not'' in the fitspw list. In future releases {{uvcontsub2}} will likely replace {{uvcontsub}}. Setting '''fitspw='0~12,16~23'''' tells {{uvcontsub2}} where to fit the continuum emission.


<source lang="python">
<source lang="python">
# In CASA
uvcontsub2(vis='g19_d2usb_targets_line.ms',fitspw='0~12,16~23',
uvcontsub2(vis='g19_d2usb_targets_line.ms',fitspw='0~12,16~23',
  solint='int',combine='spw',want_cont=True,spw='0~23')
  solint='int',combine='spw',want_cont=True,spw='0~23')
</source>
</source>


The directory g19_d2usb_targets_line.ms.cont/ contains the estimated continuum, while g19_d2usb_targets_line.ms.contsub/ contains the continuum subtracted line data
The directory g19_d2usb_targets_line.ms.cont/ contains the estimated continuum, while g19_d2usb_targets_line.ms.contsub/ contains the continuum-subtracted line data.
 
Next we image the continuum data:


<div style="background-color: #dddddd;">
<div style="background-color: #dddddd;">
'''Tip:''' for sparse mosaics like this one, the subparameter '''ftmachine='ft'''' will typically give better results than the default '''ftmachine='mosaic''''. For Nyquist sampled mosaics, '''ftmachine='mosaic'''' is the best choice, and the only choice for heterogeneous arrays. Also note that in a mosaic, the rms noise changes as a function of position. Therefore, the rms noise that you calculate from the radiometer equation is only valid where the combined primary beam response of the mosaic is near one. [[clean] Clean] weights the decision of what to clean next by the primary beam response function (PBRF), so emission near PBRF=1 will be cleaned faster than elsewhere. Thus, if you have emission away from the "sweet spot(s)" you will need to clean to a '''threshold''' lower than you might expect from the anticipated rms noise. Development is on-going to provide other options which may work better for the case of strong emission near image boundaries.
'''Tip:''' for sparse mosaics like this one, the subparameter '''ftmachine='ft'''' will typically give better results than the default '''ftmachine='mosaic''''. For Nyquist sampled mosaics, '''ftmachine='mosaic'''' is the best choice, and the only choice for heterogeneous arrays. Also note that in a mosaic, the rms noise changes as a function of position. Therefore, the rms noise that you calculate from the radiometer equation is only valid where the combined primary beam response of the mosaic is near one. {{clean}} weights the decision of what to clean next by the primary beam response function (PBRF), so emission near PBRF=1 will be cleaned faster than elsewhere. Thus, if you have emission away from the "sweet spot(s)" you will need to clean to a '''threshold''' lower than you might expect from the anticipated rms noise. Development is on-going to provide other options which may work better for the case of strong emission near image boundaries.
</div>
</div>


<source lang="python">
<source lang="python">
# In CASA
clean(vis='g19_d2usb_targets_line.ms.cont',imagename='g19_d2usb_cont',
clean(vis='g19_d2usb_targets_line.ms.cont',imagename='g19_d2usb_cont',
       field='',spw='0~23',
       field='',spw='0~23',
Line 670: Line 810:
       phasecenter="J2000 18h25m56.09 -12d04m28.20",
       phasecenter="J2000 18h25m56.09 -12d04m28.20",
       pbcor=F,minpb=0.2)
       pbcor=F,minpb=0.2)
viewer('g19_d2usb_cont.image')
</source>
</source>


Here you can see one stronger continuum source to the NE (18:25:58.5 -12:03:58.9) and one weak one to the SW (18:25:52, -12:04:53). The weak one is only visible after one round of cleaning. Because the mosaic is sparse you should make a clean mask to get the best cleaned image.
Here you can see one stronger continuum source to the NE (18:25:58.5 -12:03:58.9) and one weak one to the SW (18:25:52, -12:04:53). The weak one is only visible after one round of cleaning. Because the mosaic is sparse you should make a clean mask to get the best cleaned image.


Make a dirty image of the continuum subtracted data, excluding the line spw's where the line is. You should see nothing but noise.
[[File:g19_contsub.png|thumb| Viewer display of the continuum-subtracted data, excluding line spws.]]
 
Make a dirty image of the continuum subtracted data, excluding the spws which contain line emission. You should see nothing but noise.


<source lang="python">
<source lang="python">
# In CASA
clean(vis='g19_d2usb_targets_line.ms.contsub',imagename='g19_d2usb_contsub',
clean(vis='g19_d2usb_targets_line.ms.contsub',imagename='g19_d2usb_contsub',
       field='',spw='0~12,16~23',
       field='',spw='0~12,16~23',
Line 689: Line 830:
</source>
</source>


You can look at the image with the [[viewer]]; you should see nothing because the continuum has been subtracted.
You can look at the image with the {{viewer}}; you should see nothing because the continuum has been subtracted.


<source lang="python">
<source lang="python">
# In CASA
viewer('g19_d2usb_contsub.image')
viewer('g19_d2usb_contsub.image')
</source>
</source>
Line 697: Line 839:
== Spectral line imaging ==
== Spectral line imaging ==


This image selects a channel range suitable to image the 12CO(2-1) line.  Notice use of '''restfreq''' parameter to get velocity scale correct.
[[File:clean-2.png|thumb| Viewer display of the {{clean}}-produced image at the 17 km/s channel.]]
 
This image selects a channel range suitable to image the 12CO(2-1) line.  Notice the use of the '''restfreq''' parameter to get the velocity scale correct.


<source lang="python">
<source lang="python">
# In CASA
clean(vis='g19_d2usb_targets_line.ms.contsub',imagename='g19_d2usb_coline',
clean(vis='g19_d2usb_targets_line.ms.contsub',imagename='g19_d2usb_coline',
       field='',spw='',
       field='',spw='',
Line 712: Line 857:
       pbcor=F,minpb=0.15)
       pbcor=F,minpb=0.15)
</source>
</source>
[[File:g19COm0.png|thumb| Viewer display of the moment 0 map.]]
[[File:g19COm1.png|thumb| Viewer display of the moment 1 map.]]


It is essential to run with '''interactive=T''' and make a clean mask. A clean mask is available in the original data directory spw_fits_files. To use the pre-made mask, run:  
It is essential to run with '''interactive=T''' and make a clean mask. A clean mask is available in the original data directory spw_fits_files. To use the pre-made mask, run:  


<source lang="python">
<source lang="python">
# In CASA
clean(vis='g19_d2usb_targets_line.ms.contsub',imagename='g19_d2usb_coline',
clean(vis='g19_d2usb_targets_line.ms.contsub',imagename='g19_d2usb_coline',
       field='',spw='',
       field='',spw='',
Line 727: Line 876:
       phasecenter="J2000 18h25m56.09 -12d04m28.20",
       phasecenter="J2000 18h25m56.09 -12d04m28.20",
       pbcor=F,minpb=0.15)
       pbcor=F,minpb=0.15)
<\source>
</source>


View the image produced:  
View the image produced:  


<source lang="python">
<source lang="python">
# In CASA
viewer('g19_d2usb_coline.image')
viewer('g19_d2usb_coline.image')
</source>
</source>
Line 738: Line 888:


<source lang="python">
<source lang="python">
# In CASA
# 0th moment map
# 0th moment map
immoments(imagename='g19_d2usb_coline.image',moments=[0],axis='spectral',
immoments(imagename='g19_d2usb_coline.image',moments=[0],axis='spectral',
Line 755: Line 906:
</source>
</source>


By clicking on "folder" icon in [[viewer]] window you can bring up the moment0 as a contour image. The strange appearence of the the mom1 map toward the NE source is probably due to more than one outflow being superposed.
By clicking on the "folder" icon in the {{viewer}} window you can bring up the moment0 as a contour image. The strange appearance of the the moment1 map toward the NE source is probably due to more than one outflow being superposed.
 
Note that one could and should correct the image for the primary beam response using {{immath}} to divide the image by the .flux image. One would then need to make a mask for the moment images in order to hide the ugly stuff at the edges.
 


Note that one could and should correct the image for the primary beam response using [[immath]] to divide the image by the .flux image. One would then need to make a mask for the moment images in order to hide the ugly stuff at the edges.
{{Checked 3.1.0}}

Latest revision as of 11:46, 18 May 2011

SMA Intro

Integrated intensity SMA 12CO 2-1 map of G19.3+0.07 overlaid on a Spitzer image

NOTE: This script works with CASA 3.1.0. For an updated script that works with CASA 3.2.0, please go to SMA CO Line Data 3.2.

This script describes the data reduction of an SMA six-pointing mosaic of 12CO 2-1 emission in the infrared dark cloud G19.3+0.07. The CO line was placed in the USB. These data are from Brogan et al. in prep.; please do not use them for scientific purposes. This tutorial is based on the SMA CO Line Data script found in the CASA Scripts and Data page.

The data were previously filled in Miriad. A Tsys correction was applied, and the data were then written out in uvfits format. In this tutorial, CASA is used to read the uvfits files, and flag, calibrate, continuum subtract and image the data. Plotting for datasets with many spw such as SMA data (plotxy) can be slow. Where possible the new plotter plotms has been used.

You will need the tarball of the uvfits data, which can be downloaded here.

This tutorial uses the plotxy task in some places rather than plotms. Plotxy will eventually 
be phased out, and this tutorial will be updated. In the meantime, this tutorial will work with 
CASA Version 3.1.0.

Creating the measurement set from uvfits data

First, remove any files from previous run-throughs and unpack the tarball:

# In CASA
# Remove any files from previous run-throughs
os.system('rm -rf spw_fits.files')
os.system('tar xvf spw_fits_files.tar')
os.system('rm -rf g19_d2usb*')

Read in each uvfits file using the importuvfits filler (you need to set the upper range of the for loop to one more than you actually have, to satisfy ipython peculiarity). We create a list of the uvfits filenames in myfiles that will be used to create a single measurement set (ms).

# In CASA
os.system('rm -rf spw_ms_files')
os.system('mkdir spw_ms_files')
myfiles = []
for i in range(1,25):
    msfile = "spw_ms_files/g19_d2usb.spw"+str(i)+".ms"
    importuvfits(fitsfile="spw_fits_files/g19_d2usb.spw"+str(i),
                 vis=msfile)
    myfiles.append(msfile)
The beginning of the output from listobs.

Using concat, create a single ms from the filenames listed in myfiles:

# In CASA
concat(vis=myfiles,concatvis='g19_d2usb.ms',timesort=True)

Finally, initialize the scratch columns using clearcal and use listobs to look at the contents of the measurement set (which appear in the CASA Log Messages window):

# In CASA
clearcal(vis='g19_d2usb.ms')
listobs(vis='g19_d2usb.ms',verbose=False)
##########################################
##### Begin Task: listobs            #####
listobs::::casa
================================================================================
MeasurementSet Name:  /export/data_2/rfriesen/casa_tutorials/test_for_stable/sma_co/sma_test/g19_d2usb.ms      MS Version 2
================================================================================
Observer: SmaUser     Project:   
Observation: SMA(8 antennas)
Data records: 571536       Total integration time = 33508 seconds
Observed from   19-Aug-2005/04:23:41.9   to   19-Aug-2005/13:42:09.9 (UTC)
Fields: 14
      ID   Code Name          RA               Decl          Epoch   SrcId 
      0         jupiter       13:01:25.48      -05.18.53.00  J2000   0
      1         1743-038      17:43:58.85      -03.50.04.62  J2000   1
      2         g19.3a        18:25:58.70      -12.03.57.80  J2000   2
      3         g19.3b        18:25:57.40      -12.04.18.50  J2000   3
      4         g19.3c        18:25:56.09      -12.04.28.20  J2000   4
      5         g19.3d        18:25:54.60      -12.04.23.10  J2000   5
      6         1908-201      19:11:09.65      -20.06.55.11  J2000   6
      7         g19.3e        18:25:54.80      -12.04.43.80  J2000   7
      8         g19.3f        18:25:53.30      -12.04.51.00  J2000   8
      9         3c454.3       22:53:57.74      +16.08.53.56  J2000   9
      10        0359+509      03:59:29.75      +50.57.50.15  J2000   10
      11        hip19167      04:06:35.04      +50.21.04.54  J2000   11
      12        polaris       02:31:48.70      +89.15.50.72  J2000   12
      13        uranus        22:44:40.51      -08.50.23.77  J2000   13
       (nVis = Total number of time/baseline visibilities per field) 
Spectral Windows:(24 unique spectral windows and 1 unique polarization setups)
      SpwID  #Chans Frame Ch1(MHz)    ChanWid(kHz)TotBW(kHz)  Ref(MHz)    Corrs
      0         128 LSRK  229347.598  812.5       104000      229347.598  XX
      1         128 LSRK  229429.601  812.5       104000      229429.601  XX
      2         128 LSRK  229504.79   812.5       104000      229504.79   XX
      3         128 LSRK  229586.792  812.5       104000      229586.792  XX
      4         128 LSRK  229675.607  812.5       104000      229675.607  XX
      5         128 LSRK  229757.609  812.5       104000      229757.609  XX
      6         128 LSRK  229832.798  812.5       104000      229832.798  XX
      7         128 LSRK  229914.8    812.5       104000      229914.8    XX
      8         128 LSRK  230003.615  812.5       104000      230003.615  XX
      9         128 LSRK  230085.617  812.5       104000      230085.617  XX
      10        128 LSRK  230160.806  812.5       104000      230160.806  XX
      11        128 LSRK  230242.808  812.5       104000      230242.808  XX
      12        128 LSRK  230331.623  812.5       104000      230331.623  XX
      13        128 LSRK  230413.625  812.5       104000      230413.625  XX
      14        128 LSRK  230488.814  812.5       104000      230488.814  XX
      15        128 LSRK  230570.817  812.5       104000      230570.817  XX
      16        128 LSRK  230659.631  812.5       104000      230659.631  XX
      17        128 LSRK  230741.634  812.5       104000      230741.634  XX
      18        128 LSRK  230816.823  812.5       104000      230816.823  XX
      19        128 LSRK  230898.825  812.5       104000      230898.825  XX
      20        128 LSRK  230987.639  812.5       104000      230987.639  XX
      21        128 LSRK  231069.641  812.5       104000      231069.641  XX
      22        128 LSRK  231144.831  812.5       104000      231144.831  XX
      23        128 LSRK  231226.833  812.5       104000      231226.833  XX

The sources contained in this measurement set are:

Source Name        Description                                    ID
====================================================================
3c454.3            Bandpass and flux calibrator                    9
1743-038           Gain calibrator                                 1
1908-201           Gain calibrator                                 6
Uranus             Flux calibrator                                13
g19.7              Science target (6-pointing mosaic)    2,3,4,5,7,8
====================================================================

Other sources in the ms were not part of this program.

Initial inspection and flagging

Amplitude vs. time for the calibrator data, pre-flagging.

Use plotms to inspect the data.

# In CASA
clearstat()
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='1~9,13',
       avgchannel='3072')

There are a few bad data points in this ms. You can use the locate features in plotms to find the bad data points on field 8 and antenna ids 1-2 and 0-2. Flagging can be done in plotms, or in a script using the flagdata command. Here we will do all flagging within the script.

Tip: locate gives antenna ids, but flagdata looks first at the name column. Since SMA antenna names are numbers, this can be confusing. The corresponding names are antennas 2-3 and 1-3, which you can check with listobs.

Amplitude vs. time for the calibrator data, after flagging bad data points.
# In CASA
flagdata(vis='g19_d2usb.ms',mode='manualflag',field='8',spw='20~23',
	 timerange='07:52:45~07:52:47',antenna='2&3;1&3')

Plotting again shows that the bad data are gone:

# In CASA
clearstat()
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='8',
       avgchannel='3072')

Tip: to confirm the data were flagged in plotms, you need to change selection in plot (for example change field and then change back), and then change back to original to clear the cache and force a re-read of the data.

Amplitude vs. frequency for the bandpass calibrator 3c454.3.

Make an inital amplitude vs. frequency plot of the bandpass calibrator 3c454.3, averaging the data together in time. Here, we colorize the points by spw in the plotms call using the coloraxis parameter. You can also do this under the 'Display' tab in the plotms window.

# In CASA
clearstat()
plotms(vis="g19_d2usb.ms",xaxis="freq",yaxis="amp",field="9",
       averagedata=True,avgscan=True,avgtime='500',avgbaseline=True,coloraxis="spw")

As you can see, SMA spw overlap by a fair amount. Here we will flag the first seven channels on either side of each spw to aid calibration:

# In CASA
default('flagdata')
flagdata(vis='g19_d2usb.ms', mode='manualflag',spw='0~23:0~6;121~127')
Antenna positions for these data.

On many baselines the first integration of every scan is low. This can be seen with plotms.

# In CASA
clearstat()
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",
       field='1,6,9',avgchannel='3072',avgspw=T,coloraxis="field")

Here, we have colored the points by field to see that the low amplitudes are found in the first integrations of many scans. These can also be flagged with flagdata. The integration time of these data is 30.9s. Note that more recent SMA data rarely show this problem.

Amplitude vs. Time before flagging first integrations.
# In CASA
flagdata(vis="g19_d2usb.ms",mode="quack",field='',
	 spw="0~23",quackinterval=40.0)

We can use plotxy to plot the antenna positions. For these observations, antenna 8 was in the barn for repair.

# In CASA
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="x")

Tip: after you run plotxy, run clearstat() before trying to run plotms or you will get a table lock.

Bandpass Calibration

Phase vs. channel for the bandpass calibrator 3c454.3.

Next plot phase vs. channel. Here, we use plotms to iterate over the SMA baselines. While multiple subplots are not yet operational within plotms, you can view the individual phase vs. channel plots using the arrow buttons at the bottom right of the plotms window. By setting the xaxis parameter to 'frequency', we can see how the spectral windows overlap. Note the large discrepancy in phase for spw=4-7 on some baselines - this will be important later.

clearstat()
plotms(vis="g19_d2usb.ms",xaxis="frequency",yaxis="phase",
       ydatacolumn="data",field="9",iteraxis="baseline",
       averagedata=True,avgtime="99999",colorize=True,coloraxis="spw")

We can also use the slower plotxy to show multiple subplots in one panel. We set the plotxy parameter crossscan=T because there are two scans on 3c454.3.

# In CASA
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="phase",datacolumn="data",
       iteration="baseline",
       antenna="",spw="",field="9",
       averagemode="vector",width="1",timebin="all",
       crossscans=T,crossbls=False,stackspw=F,
       subplot=331,plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")
Amplitude vs. time for the bandpass calibrator 3c454.3.

Now plot amplitude and phase as a function of time.

# In CASA
clearstat()
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='9',
       avgchannel='3072')

The gap in the 3C454.3 observation is due to a pause to determine pointing solutions. The flags below catch times just before and after when the amplitudes and phase are off. This flagging could also have been accomplished interactively in plotms.

# In CASA
default('flagdata')
flagdata(vis="g19_d2usb.ms",field="9",timerange='12:16:58~12:19:00')
flagdata(vis="g19_d2usb.ms",field="9",timerange='12:08:00~12:09:30')
Phase vs. time for the bandpass calibrator 3c454.3.

Here because plotms cannot do multipanel iteration yet we use the slower plotxy.

# In CASA
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="phase",datacolumn="data",
       iteration="baseline", antenna="",spw="0~23",field="9",
       averagemode="vector",width="allspw",timebin="0",
       subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")
Phase solutions, all spw averaged together (minus spw 4-7).

Start by determining an initial phase only solution for 3c454.3 using gaincal to take out major phase variations with time. For maximum sensitivity, we average all the spw together, with the exception of spw 4-7 which show large phase offsets from the other spw in phase vs. channel.

# In CASA
gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.pcal',
	field='9', spw='0~3,8~23', gaintype='G', minsnr=2.0,
	refant='3', calmode='p',solint='int', combine='spw')
clearstat()
plotcal(caltable="g19_d2usb.ms.pcal",xaxis="time",yaxis="phase",
	field="9",antenna="",spw="",timerange="",subplot=331,
	overplot=False,clearpanel="Auto",iteration="antennas",
	plotsymbol="o",plotcolor="blue",showgui=True,figfile="")

Next we do a bandpass calibration. There are several possible solution types. A regular 'B' type solution can work well if you have lots of S/N on your bandpass calibrator. You may wish to make this table and compare to 'BPOLY' solution further down.

# In CASA
bandpass(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.bpcal',
 	field='9', spw='0~23', bandtype='B', solint='inf', combine='scan', 
 	refant='3', gaincurve=False, opacity=0.0,
 	gaintable='g19_d2usb.ms.pcal',spwmap=[0]) 
plotcal(caltable="g19_d2usb.ms.bpcal",xaxis="freq",yaxis="amp",
	field="9",antenna="",spw="",timerange="",subplot=311,
	overplot=False,clearpanel="Auto",iteration="antenna",
	plotsymbol="o",plotcolor="blue",showgui=True,figfile="")
Bandpass solution for the bandpass calibrator 3c454.3.

A polynomial bandpass solution is more sensitive than simply fitting a polynomial to the channel based solutions 'B' after the fact. For the mm/submm, where the S/N is often low on the bandpass calibrator, 'BPOLY' is often the best choice, but every dataset is different. You may especially need to play with the degamp and degphase parameters in bandpass. The parameter combine='scan' is needed to average together the two 3c454.3 scans. The spwmap=[0] is telling it to apply the phase solution that was found from the average of all spw to each spw before determining the bandpass solution.

# In CASA
bandpass(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.bpoly7',
	field='9', spw='0~23', bandtype='BPOLY', 
	solint='inf', combine='scan', 
	refant='3', degamp=7, degphase=7,
	gaintable='g19_d2usb.ms.pcal',spwmap=[0])
clearstat()
plotcal(caltable="g19_d2usb.ms.bpoly7",xaxis="freq",yaxis="amp",
	field="9",antenna="",spw="",timerange="",subplot=311,
	overplot=False,clearpanel="Auto",iteration="antenna",
	plotsymbol="o",plotcolor="blue",showgui=True,figfile="")

Unfortunately, this kind of plot with many spw and colors takes a long time. Be patient for now, as it's being worked on. You can safely ignore where the solution gets large at edges of spw. This is showing where we flagged channels - this is also being worked on. Zoom in if necessary to see that inner portions are nice and flat.

Phase solutions for the bandpass calibrator 3c454.3, including all spw.

The next step is to apply the bandpass solution while re-calculating with gaincal the phase only solution for 3c454.3 on a short time interval (here, we use the integration time), but now including all spw (including 4~7 since antenna based solutions will have taken out phase vs.spw differences).

# In CASA
gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.pcal2',
	field='9', spw='0~23', gaintype='G', minsnr=2.0,
	refant='3', calmode='p',solint='int', combine='spw',
	gaintable=['g19_d2usb.ms.bpoly7'])
clearstat()
plotcal(caltable="g19_d2usb.ms.pcal2",xaxis="time",yaxis="phase",
	field="9",antenna="",spw="",timerange="",subplot=331,
	overplot=False,clearpanel="Auto",iteration="antennas",
	plotsymbol="o",plotcolor="blue",showgui=True,figfile="")
Amplitude solutions for the bandpass calibrator 3c454.3, including all spw.

Next we apply the antenna based bandpass, and phase only solution to solve for amplitude solutions over a longer solint to increase the S/N. The spwmap=[[0],[]] tells gaincal that the pcal solution is independent of spw while the bpoly7 solution is per spw. spwmap order must be in the same order as the tables in gaintable.

# In CASA
gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.apcal',
	field='9', spw='0~23', gaintype='G', minsnr=2.0,
	refant='3', calmode='a',solint='300s', combine='spw',
	gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.bpoly7'],
	spwmap=[[0],[]])
clearstat()
plotcal(caltable="g19_d2usb.ms.apcal",xaxis="time",yaxis="amp",
	field="9",antenna="",spw="",timerange="",subplot=331,
	overplot=False,clearpanel="Auto",iteration="antennas",
	plotsymbol="o",plotcolor="blue",showgui=True,figfile="")

Although in previous steps solutions were applied "on-the-fly", the data were not actually changed. Running applycal will fill the corrected data column in the ms. Again be careful with specification of gaintable, and keeping the same order in spwmap and gainfield. This applycal is not necessary for later steps, but allows plotting of data calibration to this point.

Gain-corrected 3c454.3 data. Notice the few residual baseline-based problems.
# In CASA
applycal(vis='g19_d2usb.ms',spw='0~23', field='9',
	gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.apcal',
	'g19_d2usb.ms.bpoly7'],
	spwmap=[[0],[0],[]],gainfield=['9','9','9'])
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp",
       datacolumn="corrected",iteration="baseline",
       antenna='',spw="0~23",field="9",
       averagemode="vector",width="1",timebin="all",crossscans=T,
       crossbls=False,stackspw=False,extendflag=F,
       extendchan="",extendspw="",extendant="",extendtime="",
       plotcolor='darkcyan',subplot=331,
       multicolor="none",figfile="")

As this plotxy shows, after antenna based bandpass, SMA data often show residual baseline based problems as a function of spw (i.e. a few abnormally low amplitudes for a few spw on a few baselines). To fix this problem we need to do a baseline based amplitude only solution per spw, but we'll need to remove an overall normalization factor first by combining all spw (this option will be incorporated into blcal one day so that only one run is necessary). Note the significant change on baseline 1-6 in the following call to plotcal.

# In CASA
blcal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.blcal1',
	field='9', spw='0~23',solint='inf', combine='spw,scan', 
	gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.apcal',
	'g19_d2usb.ms.bpoly7'], 
	calmode='a',spwmap=[[0],[0],[]])
blcal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.blcal2',
	field='9', spw='0~23',solint='inf', combine='scan', 
	gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.apcal',
	'g19_d2usb.ms.blcal1','g19_d2usb.ms.bpoly7',], 
	calmode='a',spwmap=[[0],[0],[0],[]])
clearstat()
plotcal(caltable="g19_d2usb.ms.blcal2",xaxis="freq",yaxis="amp",
	field="9",antenna="",spw="",timerange="",subplot=511,
	overplot=False,clearpanel="Auto",iteration="antenna",
	plotsymbol="o",plotcolor="blue",showgui=True,figfile="")

Now we can do a new applycal including only the second blcal2 solution. Note that applycal overwrites the corrected column, so there is no need to clearcal. Like the antenna based bandpass the blcal solution is independent for each spw, so we set spwmap = [].

Gain-corrected 3c454.3 data, with residual problems fixed.
# In CASA
applycal(vis='g19_d2usb.ms',spw='0~23', field='9',
	gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.apcal',
	'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],
	spwmap=[[0],[0],[],[]],gainfield=['9','9','9','9'])

Note how flat the full bandpass is after all this calibration:

# In CASA
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp",
       datacolumn="corrected",iteration="baseline",
       antenna='',spw="0~23",field="9",
       averagemode="vector",width="1",timebin="all",crossscans=T,
       crossbls=False,stackspw=False,extendflag=F,
       extendchan="",extendspw="",extendant="",extendtime="",
       plotcolor='darkcyan',subplot=331,
       multicolor="none",figfile="")

Inspection of gain calibrators and setting the flux scale

Apply antenna and baseline based bandpass solutions to the gain calibrators with applycal:

# In CASA
clearstat()
applycal(vis='g19_d2usb.ms',spw='0~23', field='1,6',
	gaintable=['g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],
	spwmap=[[],[]],gainfield=['9','9'])
Corrected amplitude vs. channel for 1743-038.

Have a look at the corrected data for both calibrators with channel vs. amplitude:

# In CASA
# First field 1:
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp",datacolumn="corrected",
       iteration="baseline", antenna="",spw="",field="1",
       averagemode="vector",width="1",timebin="all",
       subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")

# Next field 6:
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp",datacolumn="corrected",
       iteration="baseline", antenna="",spw="",field="6",
       averagemode="vector",width="1",timebin="all",
       subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")
Corrected amplitude vs. time for 1743-038.
Corrected phase vs. time for 1743-038.

Have a look at the corrected data for both calibrators in phase and amplitude vs. time. For the moment, unfortunately, multicolor='field' doesn't work, so this is easiest to do separately.

# In CASA
# First, amplitude vs. time for field 1: 
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",datacolumn="corrected",
       iteration="baseline", antenna="",spw="0~23",field="1",
       averagemode="vector",width="allspw",timebin="0",
       subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")

# Second, amplitude vs. time for field 6:
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",datacolumn="corrected",
       iteration="baseline", antenna="",spw="0~23",field="6",
       averagemode="vector",width="allspw",timebin="0",
       subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="red",
       overplot=F,showflags=False,interactive=True,figfile="")

# Third, phase vs. time for field 1:
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="phase",datacolumn="corrected",
       iteration="baseline", antenna="",spw="0~23",field="1",
       averagemode="vector",width="allspw",timebin="0",
       subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")

# Last, phase vs. time for field 6:
clearstat()
plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="phase",datacolumn="corrected",
       iteration="baseline", antenna="",spw="0~23",field="6",
       averagemode="vector",width="allspw",timebin="0",
       subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")

At the time of these observations, the calibrator J1911 had the most recent and relatively stable monitor data (SMA website), so we will use it for absolute flux calibration. In future tools for using resolved planets will be implemented.

1911-201
1mm   17 Jul 2005 11:38   SMA      274.5    1.88 +/-  0.13    mgurwell  
1mm   05 Sep 2005 08:10   SMA      221.4    1.99 +/-  0.12    mgurwell  
# In CASA
setjy(vis='g19_d2usb.ms', field='6', spw='0~23', fluxdensity=[2.0,0.,0.,0.])

Final gain and flux calibration for all calibrators

Next, we calculate two separate phase solutions from the data. To avoid decorrelation of the amplitude, the phase solutions fed to the amplitude calibration step should be at the shortest solint your S/N can support. Typically this will be the integration time.

# In CASA
gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.allpcal',
	field='1,6,9', spw='0~23', gaintype='G', minsnr=2.0,
	refant='3', calmode='p',solint='int', combine='spw',
	gaintable=['g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],
	spwmap=[[],[]])

Unless you chose to do some more sophisticated smoothing, the phase solutions applied to your target might as well be on the scan time of your phase calibrator(s). Without combine='scan', the solution will not cross scan boundaries even though solint='inf'. This is an easy way to "manually" interpolate.

# In CASA
gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.allpcalscan',
	field='1,6,9', spw='0~23', gaintype='G', minsnr=2.0,
	refant='3', calmode='p',solint='inf', combine='spw',
	gaintable=['g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],
	spwmap=[[],[]])
Phase solutions for 1743-038 and 1908-201.

Plot the phase solutions of the phase calibrators. Can repeat for both short and long solution intervals.

# In CASA
clearstat()
plotcal(caltable="g19_d2usb.ms.allpcalscan",xaxis="time",yaxis="phase",
	field="1,6",antenna="",spw="",timerange="",subplot=331,
	overplot=False,clearpanel="Auto",iteration="antennas",
	plotsymbol="o",plotcolor="blue",showgui=True,figfile="")

Amplitude calibration with longer solint and applying short solint phase solution.

# In CASA
gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.allapcal',
	field='1,6,9', spw='0~23', gaintype='G', minsnr=2.0,
	refant='3', calmode='ap',solint='300s', combine='spw',
	gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.bpoly7',
	'g19_d2usb.ms.blcal2'],spwmap=[[0],[],[]])
Amplitude solutions for 1743-038 and 1908-201.

Plot the amplitude solutions:

# In CASA
clearstat()
plotcal(caltable="g19_d2usb.ms.allapcal",xaxis="time",yaxis="amp",
	field="1,6",antenna="",spw="",timerange="",subplot=331,
	overplot=False,clearpanel="Auto",iteration="antennas",
	plotsymbol="o",plotcolor="blue",showgui=True,figfile="")

Derive the absolute flux calibration using fluxscale (based on setjy above):

# In CASA
fluxscale(vis='g19_d2usb.ms',caltable='g19_d2usb.ms.allapcal',
	fluxtable='g19_d2usb.ms.fluxcal',reference='6')

Apply final calibrations

The final application of calibration for each source is done independently below. Flux densities from SMA flux monitoring are provided as checks of the absolute flux calibration.

Corrected amplitude vs. time for 3C454.3.
# In CASA
applycal(vis='g19_d2usb.ms',spw='0~23', field='9',
	gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.fluxcal',
	'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]],
	gainfield=['9','9','9','9'])
clearstat()
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='9',
       avgchannel='3072',ydatacolumn='corrected')
3C454.3
1mm   19 Aug 2005 16:01   SMA      225.3   27.36 +/-  1.37    mgurwell 
Corrected amplitude vs. time for 1743-038.
# In CASA
applycal(vis='g19_d2usb.ms',spw='0~23', field='1',
	gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.fluxcal',
	'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]],
	gainfield=['1','1','9','9'])
clearstat()
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='1',
       avgchannel='3072',ydatacolumn='corrected')
J1743-038
1mm   26 Aug 2005 06:05   SMA      225.6    1.71 +/-  0.09    mgurwell 
Corrected amplitude vs. time for 1908-201.
# In CASA
applycal(vis='g19_d2usb.ms',spw='0~23', field='6',
	gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.fluxcal',
	'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]],
	gainfield=['6','6','9','9'])
clearstat()
plotms(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",field='6',
       avgchannel='3072',ydatacolumn='corrected')
J1911-201
1mm   17 Jul 2005 11:38   SMA      274.5    1.88 +/-  0.13    mgurwell  
1mm   05 Sep 2005 08:10   SMA      221.4    1.99 +/-  0.12    mgurwell  
Corrected amplitude vs. uvdistance for Uranus.

Since Uranus is resolved we apply the solutions from 3C454.3 to achieve at least some reduction in instrumental gain variations.

# In CASA
applycal(vis='g19_d2usb.ms',spw='0~23', field='13',
	gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.fluxcal',
	'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]],
	gainfield=['9','9','9','9'])
clearstat()
plotms(vis="g19_d2usb.ms",xaxis="uvdist",yaxis="amp",field='13',
       avgchannel='3072',ydatacolumn='corrected')

A UV-distance plot of the expected flux of Uranus can be obtained from the SMA webpage. In future we will implement tools to use planets for absolute flux calibration.

# In CASA
applycal(vis='g19_d2usb.ms',spw='0~23', field='2,3,4,5,7,8',
	gaintable=['g19_d2usb.ms.allpcalscan','g19_d2usb.ms.fluxcal',
	'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]],
	gainfield=['1,6','1,6','9','9'])

Split off target data for further inspection

# In CASA
split(vis="g19_d2usb.ms",outputvis='g19_d2usb_targets.ms',
      field="2,3,4,5,7,8")

To reinitialize the scratch columns for use by later tasks, we need to run clearcal for the new dataset. Then run listobs to check the target data names.

# In CASA
clearcal(vis='g19_d2usb_targets.ms')

listobs(vis="g19_d2usb_targets.ms",verbose=F)
ID   Code Name          Right Ascension  Declination   Epoch
0         g19.3a        18:25:58.70      -12.03.57.80  J2000
1         g19.3b        18:25:57.40      -12.04.18.50  J2000
2         g19.3c        18:25:56.09      -12.04.28.20  J2000
3         g19.3d        18:25:54.60      -12.04.23.10  J2000
4         g19.3e        18:25:54.80      -12.04.43.80  J2000
5         g19.3f        18:25:53.30      -12.04.51.00  J2000

Tip: After splitting, the "corrected" data column in the original vis file will become the "data" and "corrected" columns and the field IDs are renumbered.

Amplitude vs. frequency for the target data.
# In CASA
clearstat()
plotms(vis='g19_d2usb_targets.ms',xaxis="frequency",yaxis="amp",field='0',
       avgtime='1e7',avgscan=True,avgbaseline=True)

clearstat()
plotms(vis='g19_d2usb_targets.ms',xaxis="time",yaxis="amp",field='',
       avgchannel='3072',avgspw=True)

Data at times > 11:04:00 show lots of scatter: the elevation of these data are very low and are worth flagging. One can plotms with xaxis="time", yaxis="elevation" to see this.

# In CASA
clearstat()
plotms(vis='g19_d2usb_targets.ms',xaxis="time",yaxis="elevation",field='')
Elevation vs. Time for the target data.
# In CASA
flagdata(vis="g19_d2usb_targets.ms",field="",timerange='>11:04:00')
Amplitude vs. velocity for the target data.

Note that because the Doppler tracked frequency was 13CO (+ 10GHz in the USB), the velocity diplayed in in plotxy will not be correct unless you set the restfreq to the line of interest. plotms cannot do this yet.

CO(2-1) restfreq=230.53797 GHz

# In CASA
clearstat()
plotxy(vis="g19_d2usb_targets.ms",xaxis="velocity",yaxis="amp",
	datacolumn="data",iteration="",selectdata=True,
	antenna='',spw="0~23",field="0",timerange="",averagemode="vector",
        restfreq='230.53797GHz',
	width="1",timebin="all",crossscans=T,crossbls=T,
	stackspw=False,subplot=111,interactive=T,figfile="")

The other weak line detected in this sideband is not visible in a vector average plot.

A-CH3OH (8 -1 8 0 -- 7 +0 7 0) 229.75876 GHz

Continuum imaging and subtraction

Next we want to make a continuum-free line dataset. We use field 0 to pick the line free spw. You need to make a clean mask. The strongest source is in the upper NE of the mosaic and there is a much weaker source to the SW corner.

# In CASA
clearstat()
plotxy(vis="g19_d2usb_targets.ms",xaxis="channel",yaxis="amp",
	datacolumn="data",iteration="",selectdata=True,
	antenna='',spw="0~23",field="0",timerange="",averagemode="vector",
        restfreq='230.53797GHz',
	width="1",timebin="all",crossscans=T,crossbls=T,
	stackspw=False,subplot=111,interactive=T,figfile="")

The line is confined to spws 13~15.

Copy the original ms in case something goes wrong. Currently uvcontsub2 will change the input model and corrected columns, so it's best to make a back up. This behavior will change in future.

# In CASA
os.system('cp -r g19_d2usb_targets.ms g19_d2usb_targets_line.ms')

Now use new task uvcontsub2 to subtract the continuum. Unlike uvcontsub, this task will subtract the continuum estimate for spw that are not in the fitspw list. In future releases uvcontsub2 will likely replace uvcontsub. Setting fitspw='0~12,16~23' tells uvcontsub2 where to fit the continuum emission.

# In CASA
uvcontsub2(vis='g19_d2usb_targets_line.ms',fitspw='0~12,16~23',
	  solint='int',combine='spw',want_cont=True,spw='0~23')

The directory g19_d2usb_targets_line.ms.cont/ contains the estimated continuum, while g19_d2usb_targets_line.ms.contsub/ contains the continuum-subtracted line data.

Next we image the continuum data:

Tip: for sparse mosaics like this one, the subparameter ftmachine='ft' will typically give better results than the default ftmachine='mosaic'. For Nyquist sampled mosaics, ftmachine='mosaic' is the best choice, and the only choice for heterogeneous arrays. Also note that in a mosaic, the rms noise changes as a function of position. Therefore, the rms noise that you calculate from the radiometer equation is only valid where the combined primary beam response of the mosaic is near one. clean weights the decision of what to clean next by the primary beam response function (PBRF), so emission near PBRF=1 will be cleaned faster than elsewhere. Thus, if you have emission away from the "sweet spot(s)" you will need to clean to a threshold lower than you might expect from the anticipated rms noise. Development is on-going to provide other options which may work better for the case of strong emission near image boundaries.

# In CASA
clean(vis='g19_d2usb_targets_line.ms.cont',imagename='g19_d2usb_cont',
      field='',spw='0~23',
      mode='mfs',niter=200,gain=0.1,threshold=0.0,
      psfmode='clark',imagermode='mosaic',scaletype='SAULT',
      ftmachine='ft',
      interactive=T,imsize=500,cell="0.5arcsec",
      phasecenter="J2000 18h25m56.09 -12d04m28.20",
      pbcor=F,minpb=0.2)

Here you can see one stronger continuum source to the NE (18:25:58.5 -12:03:58.9) and one weak one to the SW (18:25:52, -12:04:53). The weak one is only visible after one round of cleaning. Because the mosaic is sparse you should make a clean mask to get the best cleaned image.

Viewer display of the continuum-subtracted data, excluding line spws.

Make a dirty image of the continuum subtracted data, excluding the spws which contain line emission. You should see nothing but noise.

# In CASA
clean(vis='g19_d2usb_targets_line.ms.contsub',imagename='g19_d2usb_contsub',
      field='',spw='0~12,16~23',
      mode='mfs',niter=0,gain=0.1,threshold=0.0,
      psfmode='clark',imagermode='mosaic',scaletype='SAULT',
      ftmachine='ft',
      interactive=F,imsize=500,cell="0.5arcsec",
      phasecenter="J2000 18h25m56.09 -12d04m28.20",
      pbcor=F,minpb=0.2)

You can look at the image with the viewer; you should see nothing because the continuum has been subtracted.

# In CASA
viewer('g19_d2usb_contsub.image')

Spectral line imaging

Viewer display of the clean-produced image at the 17 km/s channel.

This image selects a channel range suitable to image the 12CO(2-1) line. Notice the use of the restfreq parameter to get the velocity scale correct.

# In CASA
clean(vis='g19_d2usb_targets_line.ms.contsub',imagename='g19_d2usb_coline',
      field='',spw='',
      mode='velocity',start='-40km/s',nchan=50,width='3km/s',
      interpolation='linear',
      niter=1000,gain=0.1,threshold=0.2,
      psfmode='clark',imagermode='mosaic',scaletype='SAULT',
      ftmachine='ft',restfreq='230.53797GHz',
      interactive=T,npercycle=400,
      imsize=500,cell="0.5arcsec",
      phasecenter="J2000 18h25m56.09 -12d04m28.20",
      pbcor=F,minpb=0.15)
Viewer display of the moment 0 map.
Viewer display of the moment 1 map.

It is essential to run with interactive=T and make a clean mask. A clean mask is available in the original data directory spw_fits_files. To use the pre-made mask, run:

# In CASA
clean(vis='g19_d2usb_targets_line.ms.contsub',imagename='g19_d2usb_coline',
      field='',spw='',
      mode='velocity',start='-40km/s',nchan=50,width='3km/s',      
      niter=1000,gain=0.1,threshold=0.2,
      psfmode='clark',imagermode='mosaic',scaletype='SAULT',
      ftmachine='ft',restfreq='230.53797GHz',
      interactive=F,npercycle=400,
      imsize=500,cell="0.5arcsec",
      mask='spw_fits_files/coline.mask',
      phasecenter="J2000 18h25m56.09 -12d04m28.20",
      pbcor=F,minpb=0.15)

View the image produced:

# In CASA
viewer('g19_d2usb_coline.image')

Now make 0th and 1st moment maps.

# In CASA
# 0th moment map
immoments(imagename='g19_d2usb_coline.image',moments=[0],axis='spectral',
	  chans='12~39',outfile='g19_d2usb_coline.mom0')

viewer('g19_d2usb_coline.mom0')

# 1st moment map
immoments(imagename='g19_d2usb_coline.image',moments=[1],axis=3,
	  chans='12~39',outfile='g19_d2usb_coline.mom1',
	  includepix=[0,2])

immath(outfile='g19_d2usb_coline.mom1masked',mode='evalexpr',
       expr='"g19_d2usb_coline.mom1"["g19_d2usb_coline.mom0">9]') 

viewer('g19_d2usb_coline.mom1masked')

By clicking on the "folder" icon in the viewer window you can bring up the moment0 as a contour image. The strange appearance of the the moment1 map toward the NE source is probably due to more than one outflow being superposed.

Note that one could and should correct the image for the primary beam response using immath to divide the image by the .flux image. One would then need to make a mask for the moment images in order to hide the ugly stuff at the edges.


Last checked on CASA Version 3.1.0.