User talk:Preshanthj: Difference between revisions

From CASA Guides
Jump to navigationJump to search
← Older edit
Preshanthj (talk | contribs)
136 edits
Preshanthj (talk | contribs)
136 edits
No edit summary
 
(52 intermediate revisions by the same user not shown)
Line 1: Line 1:
Back to [[Main Page | '''CASA guides''']]
Back to [[Main Page | '''CASA guides''']]


By [https://safe.nrao.edu/wiki/bin/view/Main/HuibIntemaWebHome Huib Intema]   & Preshanth Jagannathan
By [https://www.aoc.nrao.edu/~pjaganna/Site/Home.html Preshanth Jagannathan] & [https://safe.nrao.edu/wiki/bin/view/Main/HuibIntemaWebHome Huib Intema]


{{Checked 4.3.1}}
{{Checked 4.3.1}}
Line 138: Line 138:


<source lang='python'>
<source lang='python'>
execfile( 'mytasks.py' )
execfile('mytasks.py')


fixlowband(vis='TSUB0001_pband.ms')
fixlowband(vis='TSUB0001_pband.ms')
Line 215: Line 215:


=== Ionospheric TEC Corrections ===
=== Ionospheric TEC Corrections ===
[[File:Pband_vla_tecu.png|200px|thumb|right]]


Low frequency observations are affected by the ionosphere. A delay in the signal path is introduced between the two polarization of light that varies both as a function of time and line of sight (direction dependent). The delay is proportional to the Total Electron Content (TEC) along the line of sight and is inversely proportional to the square of the frequency.GPS measurements at two different frequencies provides us with an estimate of the TEC per square metre. This correction has been implemented in CASA which we shall apply as a calibration table by means of the <i>gencal</i> task. The task requires a TEC map that we will generate utilizing casa recipes.
Low frequency observations are affected by the ionosphere. A delay in the signal path is introduced between the two polarization of light that varies both as a function of time and line of sight (direction dependent). The delay is proportional to the Total Electron Content (TEC) along the line of sight and is inversely proportional to the square of the frequency.GPS measurements at two different frequencies provides us with an estimate of the TEC per square metre. This correction has been implemented in CASA which we shall apply as a calibration table by means of the <i>gencal</i> task. The task requires a TEC map that we will generate utilizing casa recipes.
Line 221: Line 222:
from recipes import tec_maps
from recipes import tec_maps


tec_image, tec_rms_image = tec_maps.create(vis='TSUB0001_pband.ms',plot_vla_tec=T)
tec_image, tec_rms_image = tec_maps.create(vis='TSUB0001_pband.ms',doplot=True)


gencal(vis='TSUB0001_pband.ms',caltable='TSUB0001_pband.tecim',caltype='tecim',infile=tec_image)
gencal(vis='TSUB0001_pband.ms',caltable='TSUB0001_pband.tecim',caltype='tecim',infile=tec_image)
Line 246: Line 247:
== Dead and swapped antennas ==
== Dead and swapped antennas ==


Now it is time to have a first visual look at the uncalibrated visibility data. It is important to identify dead antennas / polarizations early on, so we can exclude them from further data processing (which is the most efficient data 'reduction' :-) ). A convenient way of doing this is through the <i>plotms</i> task on bright calibrator(s) (again, see the <i>listobs()</i> output), in our case 3C48.
Now it is time to have a first visual look at the uncalibrated visibility data. It is important to identify dead antennas / polarizations early on, so we can exclude them from further data processing (which is the most efficient data 'reduction' :-) ). A convenient way of doing this is through the <i>plotms</i> task on bright calibrator(s) (again, see the <i>listobs()</i> output), in our case 3C48.[[File:TSUB0001_pband_3C48_prebp.png|200px|thumb|right]]
<source lang='python'>
<source lang='python'>
plotms(vis='TSUB0001_pband.ms',xaxis='freq',yaxis='amp',antenna='ea01',correlation='XX,YY',
plotms(vis='TSUB0001_pband.ms',xaxis='freq',yaxis='amp',antenna='ea01',correlation='XX,YY', field='3C48',  
      coloraxis='spw',xlabel='Frequency',ylabel='Amplitude',plotfile='amp_freq_ea01_pband.png')
          plotrange=[0.2,0.5,0.0,100.0], coloraxis='spw',xlabel='Frequency',ylabel='Amplitude',iteraxis='baseline',
 
          plotfile='TSUB0001_pband_3C48_prebp.png')
</source>
</source>


This will bring forward a data inspection window, with (per scan, per baseline, per spectral window) the top row a set of intensity plots of the visibility amplitudes in all polarization products (XX, YY, XY, and YX) before flagging, and the bottom row the same data after applying some outlier flagging (marked in blue). The intensity plots do scale automatically to min-max, so any bright RFI pixel will skew the intensity range (hence the difference between the top and bottom row).
This will load the <i>plotms</i> window that has queued up all the baselines of antenna ea01 and is currently displaying the first amplitude vs frequency plot with the colors representing the different spectral windows. Note that we are only plotting the XX and YY correlation of the antennas. This is because we expect most of the power in the linear cross correlation products to be in these two correlations. Notice that some spectral windows in particular are badly affected due to rfi. This is the issue that we will deal with after identifying dead antennas. If you scroll through the baselines using the green forward buttons in the plotms you will notice that the amplitude is particularly low when you encounter ea10 and ea13. If we make a quick plot of the baseline made by ea10 and ea13 we can see that there is no power in any of the polarizations so these two antennas were dead and can be flagged out. We also notice that ea25 has low XX polarization but also notice the excess power in the XY polarization. This tells us that the polarization inputs have been wrongly labeled we have a local task that we built already and we will utilize that to swap the polarization labels for antenna ea25. In addition to the dead antennas we also have a dummy scan during slew on a source named 3C48D, we shall flag that scan by hand too.


In the case of 16 x 16 MHz spectral windows, we want to initially inspect a more central spectral window (spw 5 through 7) than the default one at the low edge of the band (spw = 0). Notice when clicking <i>Next SPW</i> a few times that several spectral windows are heavily affected by RFI, for instance spw 2 and 3, and 8 and 9. Also notice that not all RFI is captured; this is mainly because within a spectral window there is too little contrast between healthy and affected data. We will resolve that later.
<source lang='python'>
flagdata(vis='TSUB0001_pband.ms',mode='manual',antenna='ea10,ea13')
flagdata(vis='TSUB0001_pband.ms',mode='manual',field='3C48D')
swappol(vis='TSUB0001_pband.ms',antenna='ea25')
</source>


A healthy spectral window will have little flagged data, and show a clear difference between the parallel-hand polarizations XX,YY (bottom row left two panels) and the cross-hand polarizations XY,YX (bottom row right two panels), where the former shows a time-constant, smooth spectral behaviour, while the latter shows a noise-like behaviour. Cycling through the baselines (<i>Previous Baseline</i> / <i>Next Baseline</i> buttons) please note the following:
NOTE: <i>swappol()</i> is not part of the default CASA release, but was imported together with <i>fixlowband()</i>. It is a pure python implementation operating on visibilities, therefore it may take a while to finish. The way in which we run <i>swappol()</i> we overwrite the previous visibilities, so no new measurement set is created.
* Baselines containing antenna ea07 have a ripple in XX.
* Baselines containing antennas ea10 or ea13 contain only noise in all polarization products.
* Baselines containing antennas ea25 contain noise in XX and YY, signal in XY, and noise in YX
* Baselines containing antenna ea28 have a ripple in YY
What this means is that antennas ea10 and ea13 are dead and need to be flagged. We will also take this opportunity to flag the dummy scan 1 on '3C48D'. Antenna ea25 has swapped polarizations and a dead Y-polarization; for demonstration purposes we will unswap this antenna, but it becomes useless later on since CASA cannot form proper stokes I visibilities (which is different behaviour than in some other reduction software, e.g,. AIPS). The ripples in antennas ea07 (X-polarization) and ea28 (Y-polarization) will show up, and calibrate out, in the bandpass calibration performed below.
 
If we're done inspecting the data in <i>flagdata()</i>, press the <i>Quit</i> button (and not <i>Stop Display</i>). The next commands will perform the necessary flag operations, and unswap antenna ea25.


We also noticed sharp rfi peaks during our data examination, to prevent Gibbs ringing it is best to hanning smooth the data at this juncture before we proceed further with automatic flagging. [[File:TSUB0001_pband_3C48_prebp.png|200px|thumb|right]]
<source lang='python'>
<source lang='python'>
dead_antennas = 'ea10,ea13'
hanningsmooth(vis='TSUB0001_pband.ms',datacolumn='data')
flagdata( vis = pband_ms_path, mode = 'manual', antenna = dead_antennas, flagbackup = False )
 
dummy_scan = '1'
flagdata( vis = pband_ms_path, mode = 'manual', scan = dummy_scan, flagbackup = False )
 
swapped_antennas = 'ea25'
swappol( vis = pband_ms_path, antenna = swapped_antennas )
</source>
</source>


NOTE: <i>swappol()</i> is not part of the default CASA release, but was imported together with <i>fixlowband()</i>. It is a pure python implementation operating on visibilities, therefore it may take a while to finish. The way in which we run <i>swappol()</i> we overwrite the previous visibilities, so no new measurement set is created.
We should replot the primary calibrator to see the effect of hanning smoothing on the data and the rfi.  


You can run <i>flagdata()</i> again in <i>calculate</i> mode to inspect the results.
<source lang='python'>
plotms(vis='TSUB0001_pband.ms',xaxis='freq',yaxis='amp',antenna='ea01',correlation='XX,YY',
          field='3C48', plotrange=[0.2,0.5,0.0,100.0], coloraxis='spw',xlabel='Frequency',ylabel='Amplitude',
          iteraxis='baseline', plotfile='TSUB0001_pband_3C48_prebp_hanning.png')


Next, we will apply Hanning smoothing to suppress Gibbs ringing around bright RFI lines. Visibilities are overwritten; no new measurement set is created.
<source lang='python'>
hanningsmooth( vis = pband_ms_path, datacolumn = 'data' )
</source>
</source>


== Automatic flagging ==
== Automatic flagging ==


Before throwing away bad spectral windows by hand, we'll let the automated flaggers in CASA have a go at it. The main task for this is <i>flagdata()</i>. While working on your data, <i>flagdata()</i> produces an abundance of output, not all of which I find easy to understand. To help with this, we can get the flag status of our data before and after auto-flagging by running <i>flagdata()</i> in the <i>summary</i> mode. This first call provides the flag statistics of the flagging up to now, before any auto-flagging:
Before throwing away bad spectral windows by hand, we'll let the automated flaggers in CASA have a go at it. The main task for this is <i>flagdata()</i>. While working on your data, <i>flagdata()</i> produces an abundance of output, not all of which I find easy to understand. To help with this, we can get the flag status of our data before and after auto-flagging by running <i>flagdata()</i> in the <i>summary</i> mode. This first call provides the flag statistics of the flagging up to now, before any auto-flagging:
<source lang='python'>
<source lang='python'>
summary_1 = flagdata( vis = pband_ms_path, mode = 'summary' )
summary_1 = flagdata(vis='TSUB0001_pband.ms', mode='summary')
</source>
</source>


Line 305: Line 298:


<source lang='python'>
<source lang='python'>
flagdata( vis = pband_ms_path, field = '*', mode = 'tfcrop', datacolumn = 'data', timecutoff = 5., freqcutoff = 5., maxnpieces = 5,
flagdata(vis='TSUB0001_pband.ms', field='*', mode='tfcrop', datacolumn='data', timecutoff=4., freqcutoff=3., maxnpieces=5,
   action = 'apply', display = 'report', flagbackup = True )
   action='apply', display='report', flagbackup=True, combinescans=True, ntime='3600s', correlation='ABS_XY,ABS_YX')
flagdata( vis = pband_ms_path, field = '*', mode = 'tfcrop', datacolumn = 'data', timecutoff = 5., freqcutoff = 5., maxnpieces = 5,
 
   action = 'apply', display = 'report', flagbackup = False )
flagdata(vis='TSUB0001_pband.ms', field='*', mode='tfcrop', datacolumn='data', timecutoff=3., freqcutoff=3., maxnpieces=2,
   action='apply', display='report', flagbackup=True, combinescans=True, ntime='3600s', correlation='ABS_XY,ABS_YX')
 
flagdata(vis='TSUB0001_pband.ms', mode='extend')
</source>
</source>


Let's get another flag summary to see how much extra data got flagged:
Let's get another flag summary to see how much extra data got flagged:  
<source lang='python'>
<source lang='python'>
summary_2 = flagdata( vis = pband_ms_path, mode = 'summary' )
summary_2 = flagdata(vis='TSUB0001_pband.ms' , mode='summary')


axis = 'scan'
axis = 'scan'
Line 319: Line 315:
   old_stats = summary_1[ axis ][ value ]
   old_stats = summary_1[ axis ][ value ]
   print '%s %s: %5.1f percent flagged additionally' % ( axis, value, 100. * ( stats[ 'flagged' ] - old_stats[ 'flagged' ] ) / stats[ 'total' ] )
   print '%s %s: %5.1f percent flagged additionally' % ( axis, value, 100. * ( stats[ 'flagged' ] - old_stats[ 'flagged' ] ) / stats[ 'total' ] )
</source>


There is quite a bit of RFI that gets missed in RFI-rich spectral windows. The solution to this is to provide the flagging routines with more contrast between healthy and affected data. For this, we will combine the 16 spectral windows into one, and perform a preliminary bandpass calibration to take out the bandpass edge effects between spectral windows, and other intrinsic bandpass structure. Combining spectral windows can only be done with the experimental task <i>mstransform()</i>, which seems to work ok for this type of operation. A new measurement set will be created.
plotms(vis='TSUB0001_pband.ms',xaxis='freq',yaxis='amp',antenna='ea01',correlation='XX,YY',
<source lang='python'>
          field='3C48', plotrange=[0.2,0.5,0.0,100.0], coloraxis='spw',xlabel='Frequency',ylabel='Amplitude',
onespw_ms_path = project_name + '_PBAND_1SPW.ms'
          iteraxis='baseline', plotfile='TSUB0001_pband_3C48_prebp_hanning.png')
if os.path.isdir( pband_ms_path ):
  rmtables( onespw_ms_path )


mstransform( vis = pband_ms_path, outputvis = onespw_ms_path, datacolumn = 'data', combinespws = True )
</source>
</source>
[[File:TSUB0001_pband_3C48_prbp_tfcropped.png|200px|thumb|right]]
There is quite a bit of RFI that gets missed in RFI-rich spectral windows. The solution to this is to provide the flagging routines with more contrast between healthy and affected data. For this, we will perform a preliminary bandpass calibration to take out the bandpass shape. We will do some coarse preliminary calibration apply it to the calibrator before flagging for RFI once more.


Before we do more flagging, we will approximately remove the frequency structure across the observing band that is induced by our instrument (the VLA). This bandpass calibration, one per antenna and polarization, is derived using the primary calibrator 3C48. In the measurement set, all (uncalibrated) visibility data is located in a column named <i>DATA</i>. To perform any kind of calibration, we need to create and fill a <i>MODEL_DATA</i> column that matches the dimensions of the <i>DATA</i> column. Logically, the <i>MODEL_DATA</i> column contains the visibilities corresponding to whatever radio source model we choose to calibrate against. For our bandpass calibration against 3C48, we assume an unpolarized point source model with a frequency-dependent flux density as described by the Scaife & Heald (2012) low-frequency flux scale. Note that this operation roughly doubles the size of the measurement set.
<source lang='python'>
<source lang='python'>
fluxscale = setjy( vis = onespw_ms_path, field = bcal_field, scalebychan = True, standard = 'Scaife-Heald 2012', usescratch = True )
</source>


There will probably be some warnings about a time-dependent feed table in this and subsequent CASA task output logs. This is caused by a problem in the <i>FEED</i> sub-table created by <i>mstransform()</i>, but can be ignored. The logger reports the stokes I flux of 3C48, which should be in the range of 50-60 Jy. For P-band you really need calibrators of this magnitude. Any calibrator will be competing with typically 40-50 Jy of total apparent flux of other sources within the same wide field-of-view. Experience shows that observing a 20-30 Jy source like 3C286 for a few minutes is insufficient to derive good-quality bandpass calibrations. Better options are 3C48 and 3C147.
gaincal(vis='TSUB0001_pband.ms', caltable='TSUB0001_pband.G0', gaintype='G', calmode='p', solint='int', field='3C48',refant='ea09',
  gaintable=['TSUB0001_pband.antpos','TSUB0001_pband.rq','TSUB0001_pband.tecim'])
 
gaincal(vis='TSUB0001_pband.ms', caltable='TSUB0001_pband.K0', gaintype='K', solint='inf', field='3C48',refant='ea09',
  gaintable=['TSUB0001_pband.antpos','TSUB0001_pband.rq','TSUB0001_pband.tecim','TSUB0001_pband.G0'])


We set the minimum signal-to-noise (SNR) requirement to 0.1 to allow for as many solutions as possible, even for channels that have much RFI. Our aim here is not to create a clean bandpass calibration, but to have a first-order correction to perform more flagging.
bandpass(vis='TSUB0001_pband.ms', caltable='TSUB0001_pband.bp0', solint='inf', field='3C48',refant='ea09', minsnr=2.0,
<source lang='python'>
  gaintable=['TSUB0001_pband.antpos','TSUB0001_pband.rq','TSUB0001_pband.tecim','TSUB0001_pband.G0','TSUB0001_pband.K0'])
bcal_field = '3C48'
bcal_scan = '4'
reference_antenna = 'ea01'
bandpass_table_v0 = onespw_ms_path + '.B0'


if os.path.isdir( bandpass_table_v0 ):
applycal(vis='TSUB0001_pband.ms', field='*', applymode='calflagstrict',
   rmtables( bandpass_table_v0 )
   gaintable=['TSUB0001_pband.antpos','TSUB0001_pband.rq','TSUB0001_pband.tecim','TSUB0001_pband.G0','TSUB0001_pband.K0'] )


bandpass( vis = onespw_ms_path, caltable = bandpass_table_v0, field = bcal_field, selectdata = True, refant = reference_antenna,
  solint = 'inf', combine = 'scan,field', minsnr = 0.1, solnorm = False, parang = True )
</source>
</source>


The task <i>bandpass()</i> will report failed solutions for several channels, which is due to our previous flagging and because of heavily RFI-affected channels. The resulting bandpass calibration table can be inspected using <i>plotcal()</i>. It may be necessary to adjust the <i>plotrange</i> parameters ( <i>min_freq_MHz</i>, <i>max_freq_MHz</i>, <i>min_amp</i>, <i>max_amp</i> ) to get an optimal view for all antennas.
We will now flag the corrected data column that contains the coarsely calibrated visibilities which provide better contrast to the flagging algorithms to remove the RFI present. We will begin by running the <i>RFLAG</i> algorithm in the task <i>flagdata</i>. We will yet again in summary mode see how much more of the observed data was flagged.
<source lang='python'>
<source lang='python'>
plotcal( caltable = bandpass_table_v0, xaxis = 'freq', yaxis = 'amp', iteration = 'antenna', plotrange = [ 200.,500.,0.,0.1 ], markersize = 2.0 )
flagdata(vis='TSUB0001_pband.ms', field='*', mode='rflag', datacolumn='corrected', timedevscale=4., freqdevscale=3.,
</source>
  action='apply',  flagbackup=True, combinescans=True, ntime='3600s', maxnpieces=5)


Each plot shows the bandpass amplitudes for both polarizations X and Y. The Mark Region and Locate buttons can be used to identify the polarization, pol = 0 being X and pol = 1 being Y. You will notice that each antenna and polarization has a different overall gain level. Spectral windows are recognizable by the depressions in the bandpass plots, where spw = 2,9,13 are particularly affected by RFI. Several antennas / polarizations show significant ripples across frequency, which are mostly caused by cable connector problems (e.g., see antenna ea07). More information is provided in the Appendix below. As expected, antenna ea25 has one polarization that has a very low-amplitude bandpass plot.
flagdata(vis='TSUB0001_pband.ms', field='*', mode='rflag', datacolumn='corrected', timedevscale=4., freqdevscale=3.,
 
  action='apply',  flagbackup=True, combinescans=True, ntime='3600s', maxnpieces=5)
<source lang='python'>
summary_3 = flagdata(vis='TSUB0001_pband.ms' , mode='summary')
summary_3 = flagdata( vis = onespw_ms_path, mode = 'summary' )


axis = 'scan'
axis = 'scan'
for id, stats in summary_3[ axis ].iteritems():
for value, stats in summary_2[ axis ].iteritems():
   print '%s %s: %5.1f percent flagged' % ( axis, id, 100. * stats[ 'flagged' ] / stats[ 'total' ] )
  old_stats = summary_1[ axis ][ value ]
</source>
   print '%s %s: %5.1f percent flagged additionally' % ( axis, value, 100. * ( stats[ 'flagged' ] - old_stats[ 'flagged' ] ) / stats[ 'total' ] )


Although we could, we are not going to flag any bad bandpass calibration solutions in =plotcal()= (we'll do that later). The bandpass table, with all its obvious problems, is temporarily applied to the data to facilitate more flagging. This means that the visibility data (of all sources) in the <i>DATA</i> column gets corrected with the bandpass calibration table and copied into the <i>CORRECTED_DATA</i> column. Channels for which no solution was found will be flagged. Note that this operation, together with the <i>MODEL_DATA</i> column created by <i>setjy()</i>, roughly triples the size of the measurement set with respect to the original. This may take a while.
plotms(vis='TSUB0001_pband.ms',xaxis='freq',yaxis='amp',antenna='ea01',correlation='XX,YY',
 
          field='3C48', plotrange=[0.2,0.5,0.0,100.0], coloraxis='spw',xlabel='Frequency',ylabel='Amplitude',
<source lang='python'>
          iteraxis='baseline', plotfile='3C48_postbandpass_rflag.png', ydatacolumn='corrected')
applycal( vis = onespw_ms_path, field = '*', gaintable = [ bandpass_table_v0 ], interp = [ 'nearest,nearestflag' ],
  applymode = 'calflagstrict', flagbackup = True )
</source>
</source>
[[File:3c48_postbandpass_rflag.png|200px|thumb|right]]
Taking another visual look at the data using <i>plotms</i> we examine the effect of the automated flagging on the calibrator. We find that most spectral windows are RFI free and look great but for spectral windows 2, 9, 13, 15. We dig deeper into the cause of the RFI in these instance by using the locate tool in the <i>plotms</i> tool. It comes to our notice that most of the RFI in spw 13 comes from antenna ea04. So we can flag ea04 for spa 13. This still leaves 3 problematic spw's 2,9,15. These spectral windows have been designed to catch most of the RFI in the P-Band and in the event that they are terribly contaminated we can ignore these spw's and proceed with the actual calibration. Nonetheless a final pass of <i>flagdata</i> in RFLAG mode over the contaminated spw's does not hurt.


For convenience, we have made a backup of the flagged data status. Now that we have flattened out most of the bandpass structure, we run <i>flagdata()</i> in the <i>rflag</i> mode on the corrected data, which performs a more rigorous outlier detection. First we check the results before applying them. Remember that the data now only has one spectral window.
<source lang='python'>
<source lang='python'>
flagdata( vis = onespw_ms_path, field = '*', mode = 'rflag', datacolumn = 'corrected', winsize = 5,
flagdata(vis='TSUB0001_pband.ms', field='3C48', mode='rflag', datacolumn='corrected', timedevscale=4., freqdevscale=3.,
  action = 'calculate', display = 'both', flagbackup = False )
   action='apply', flagbackup=True, combinescans=True, ntime='3600s', maxnpieces=5, spw='2,9,15')
</source>
 
You can see that most of antenna ea25 gets flagged. This is likely because much of the dead polarization data gets flagged, and these flags get extended to the other polarization products. The task =flagdata()= in =rflag= mode does a good job in catching most of the obvious RFI without over-flagging, so we'll now run it twice to apply the flags.
 
<source lang='python'>
flagdata( vis = onespw_ms_path, field = '*', mode = 'rflag', datacolumn = 'corrected', winsize = 5,
   action = 'apply', display = 'report', flagbackup = False )
flagdata( vis = onespw_ms_path, field = '*', mode = 'rflag', datacolumn = 'corrected', winsize = 5,
  action = 'apply', display = 'report', flagbackup = False )
</source>
</source>


By getting the new flagging statistics, we can indeed confirm that ea25 is now also completely flagged.
==Flux Density Calibration==
The initial flagging has cleaned up most of the stray RFI across the band. This allows us to proceed with the actual calibration of the data. Before we get to the calibration tables its essential to do the flux density calibration of our calibrator 3C48. This is done by means of the <i>setjy</i> task. Before we run the task on hand we first clear the preliminary calibration that was carried out to enable better flagging. We do that but running the <i>clearcal</i> task in casa.


<source lang='python'>
<source lang='python'>
summary_4 = flagdata( vis = onespw_ms_path, mode = 'summary' )
clearcal(vis='TSUB0001_pband.ms')
 
axis = 'antenna'
for value, stats in summary_4[ axis ].iteritems():
  old_stats = summary_3[ axis ][ value ]
  print '%s %s: %5.1f percent flagged' % ( axis, value, 100. * ( stats[ 'flagged' ] ) / stats[ 'total' ] )
 
axis = 'scan'
for value, stats in summary_4[ axis ].iteritems():
  old_stats = summary_3[ axis ][ value ]
  print '%s %s: %5.1f percent flagged additionally' % ( axis, value, 100. * ( stats[ 'flagged' ] - old_stats[ 'flagged' ] ) / stats[ 'total' ] )
</source>


Now we have removed most of the RFI, we will start the actual calibration of the data. We go back to our original data in the <i>DATA</i> column by removing the the <i>CORRECTED_DATA</i> column. For this, we use the table toolkit.
setjy(vis='TSUB0001_pband.ms', standard='Scaife-Heald 2012', field='3C48')
<source lang='python'>
tb.open( onespw_ms_path, nomodify = False )
tb.removecols( 'CORRECTED_DATA' )
tb.close()
</source>
</source>
[...]


== Delay and bandpass calibration ==
== Delay and bandpass calibration ==


The initial flagging has cleaned the data of most of the RFI. This allows us to proceed with deriving the necessary calibrations. First we use the full bandwidth on a single scan on the primary calibrator to determine a single delay per antenna per polarization. This will determine a single (approximate) phase slope across frequency, mainly caused by propagation effects in the (time-variable) ionosphere, cable length differences and electronics in the signal paths from antenna feeds to correlator. Note that for a short (5-10 minutes), single scan on the calibrator we can get away with solving for a time-invariant delay per antenna per polarization. We will again use scan 4, which is the longest scan on 3C48.
First we use the full bandwidth on a single scan on the primary calibrator to determine a single delay per antenna per polarization. This will determine a single (approximate) phase slope across frequency, mainly caused by propagation effects in the (time-variable) ionosphere, cable length differences and electronics in the signal paths from antenna feeds to correlator. Note that for a short (5-10 minutes), single scan on the calibrator we can get away with solving for a time-invariant delay per antenna per polarization. We will again use scan 4, which is the longest scan on 3C48.


<source lang='python'>
<source lang='python'>
bcal_field = '3C48'
gaincal(vis='TSUB0001_pband.ms', caltable='TSUB0001_pband.K1', field='3C48', solint='inf', refant='ea09', gaintype='K', parang=True)
bcal_scan = '4'
reference_antenna = 'ea01'
delay_table_v1 = onespw_ms_path + '.K1'


if os.path.isdir( delay_table_v1 ):
plotcal(caltable='TSUB0001_pband.K1', xaxis='antenna', yaxis='delay', markersize=2.0)
  rmtables( delay_table_v1 )
 
gaincal( vis = onespw_ms_path, caltable = delay_table_v1, field = bcal_field, selectdata = True, scan = bcal_scan,
  solint = 'inf', refant = reference_antenna, minsnr = 3.0, gaintype = 'K', parang = True )
 
plotcal( caltable = delay_table_v1, xaxis = 'antenna', yaxis = 'delay', markersize = 2.0 )
</source>
</source>


Since we picked ea01 as our reference antenna, the delays for this antenna ID (=0) are arbitrarily set to zero. All other delays should be within 100 nanosec or so (which they should be for this data set). Larger values should be treated with suspicion, probably indicating a problem with the antenna / polarization. We will use the bandpass calibration (next) to verify this.
Since we picked ea09 as our reference antenna, the delays for this antenna ID (=0) are arbitrarily set to zero. All other delays should be within 30 nanosec or so (which they should be for this data set). Larger values should be treated with suspicion, probably indicating a problem with the antenna / polarization. We will use the bandpass calibration (next) to verify this.


Next is to determine the bandpass calibration. We use the same source and scan, and apply the delay calibration before solving for the bandpass. Note that we request a minimum SNR of 3, which will make the solve fail for the worst channels, but not all.
Next is to determine the bandpass calibration. We use the same source and scan, and apply the delay calibration before solving for the bandpass. Note that we request a minimum SNR of 3, which will make the solve fail for the worst channels, but not all.
<source lang='python'>
<source lang='python'>
bandpass_table_v1 = onespw_ms_path + '.B1'
if os.path.isdir( bandpass_table_v1 ):
  rmtables( bandpass_table_v1 )


bandpass( vis = onespw_ms_path, caltable = bandpass_table_v1, field = bcal_field, selectdata = True, scan = bcal_scan,
bandpass(vis='TSUB0001_pband.ms', caltable='TSUB0001_pband.B1', field='3C48', solint='inf', refant='ea09', minsnr=3.0, parang = True,
  solint = 'inf', refant = reference_antenna, minsnr = 3.0, solnorm = False, parang = True,
   gaintable=['TSUB0001_pband.K1'], interp=['nearest,nearestflag'])
   gaintable = [ delay_table_v1 ], interp = [ 'nearest,nearestflag' ] )
</source>
</source>


After solving, we inspect the bandpass calibration amplitudes and phases, and flag any obvious residual outliers. Outliers may be found above or below the average bandpass curves, and tend to arise in the same channel (frequency) ranges for all antennas. Don't spend more than 1 minute per antenna. And don't worry about missing some bad points; these bad channels will likely get flagged in a later stage anyway. Note that the bandpass amplitudes and phases may differ between polarizations of the same antenna. Also note that the bandpass phases across frequency should not have an overall gradient (this was removed by the delay calibration). Also also note that bandpass phases may wrap around from +/-180 to -/+180 degrees. If that is problematic, re-run the phase <i>plotcal()</i> with <i>plotrange = [ 200.,500.,0.,360. ]</i>. Also note that there are no bandpass solutions for antennas ea10, ea13 and ea25.
After solving, we inspect the bandpass calibration amplitudes and phases, and flag any obvious residual outliers. Outliers may be found above or below the average bandpass curves, and tend to arise in the same channel (frequency) ranges for all antennas. Don't spend more than 1 minute per antenna. And don't worry about missing some bad points; these bad channels will likely get flagged in a later stage anyway. Note that the bandpass amplitudes and phases may differ between polarizations of the same antenna. Also note that the bandpass phases across frequency should not have an overall gradient (this was removed by the delay calibration). Also also note that bandpass phases may wrap around from +/-180 to -/+180 degrees. If that is problematic, re-run the phase <i>plotcal()</i> with <i>plotrange = [ 200.,500.,0.,360. ]</i>. Also note that there are no bandpass solutions for antennas ea10, ea13 and ea25.
<source lang='python'>
<source lang='python'>
plotcal( caltable = bandpass_table_v1, xaxis = 'freq', yaxis = 'amp', iteration = 'antenna', markersize = 2.0 )
plotcal(caltable='TSUB0001_pband.B1', xaxis='freq', yaxis='amp', iteration='antenna', markersize=2.0)


plotcal( caltable = bandpass_table_v1, xaxis = 'freq', yaxis = 'phase', iteration = 'antenna', plotrange = [ 200.,500.,-180.,180. ], markersize = 2.0 )
plotcal(caltable='TSUB0001_pband.B1', xaxis='freq', yaxis='phase', iteration='antenna', plotrange=[ 200.,500.,-180.,180. ], markersize = 2.0)
</source>
</source>


Line 460: Line 408:


<source lang='python'>
<source lang='python'>
gain_table_v1 = onespw_ms_path + '.G1'
gaincal(vis='TSUB0001_pband.ms', caltable='TSUB0001_pband.G1', field='3C48', solint = 'int', refant = reference_antenna, minsnr = 3.0, gaintype = 'G',
 
calmode = 'ap', gaintable = [ delay_table_v1, bandpass_table_v1 ], interp = ['','','','nearest,nearestflag', 'nearest,nearestflag' ], parang = True )
if os.path.isdir( gain_table_v1 ):
  rmtables( gain_table_v1 )
 
gaincal( vis = onespw_ms_path, caltable = gain_table_v1, field = bcal_field, selectdata = True, scan = bcal_scan,
  solint = 'int', refant = reference_antenna, minsnr = 3.0, gaintype = 'G', calmode = 'ap',
  gaintable = [ delay_table_v1, bandpass_table_v1 ], interp = [ 'nearest,nearestflag', 'nearest,nearestflag' ], parang = True )
</source>
</source>


Line 473: Line 415:


<source lang='python'>
<source lang='python'>
plotcal( caltable = gain_table_v1, xaxis = 'time', yaxis = 'amp', iteration = 'antenna', markersize = 2.0 )
plotcal(caltable='TSUB0001_pband.ms', xaxis='time', yaxis='amp', iteration='antenna', markersize=2.0 )


plotcal( caltable = gain_table_v1, xaxis = 'time', yaxis = 'phase', iteration = 'antenna', plotrange = [ None,None,-180.,180. ], markersize = 2.0 )
plotcal(caltable='TSUB0001_pband.ms', xaxis='time', yaxis='phase', iteration='antenna', plotrange=[ None,None,-180.,180.], markersize=2.0)
</source>
</source>


Line 482: Line 424:
Next, we will prepare a smoothed and interpolated version of the gain calibration table, which will be applied later on to the target field data. This prevents flagging of target field data when one of the edges of a calibrator scan is flagged for one or more antennas. The result can be inspected with <i>plotcal</i>.
Next, we will prepare a smoothed and interpolated version of the gain calibration table, which will be applied later on to the target field data. This prevents flagging of target field data when one of the edges of a calibrator scan is flagged for one or more antennas. The result can be inspected with <i>plotcal</i>.
<source lang='python'>
<source lang='python'>
smooth_gain_table_v1 = onespw_ms_path + '.Gs1'


if os.path.isdir( smooth_gain_table_v1 ):
smoothcal(vis='TSUB0001_pband.ms', tablein='TSUB0001_pband.G1', caltable='TSUB0001_pband.Gs1', smoothtype='median', smoothtime = 60.*60.)
  rmtables( smooth_gain_table_v1 )


smoothcal( vis = onespw_ms_path, tablein = gain_table_v1, caltable = smooth_gain_table_v1, smoothtype = 'median', smoothtime = 60. * 60. )
plotcal(caltable='TSUB0001_pband.Gs1', xaxis='time', yaxis='amp', iteration='antenna', plotrange=[ None,None,0.8,1.2 ], markersize=2.0)


plotcal( caltable = smooth_gain_table_v1, xaxis = 'time', yaxis = 'amp', iteration = 'antenna', plotrange = [ None,None,0.8,1.2 ], markersize = 2.0 )
plotcal(caltable='TSUB0001_pband.Gs1', xaxis='time', yaxis='phase', iteration='antenna', plotrange=[ None,None,-10.,10. ], markersize=2.0)


plotcal( caltable = smooth_gain_table_v1, xaxis = 'time', yaxis = 'phase', iteration = 'antenna', plotrange = [ None,None,-10.,10. ], markersize = 2.0 )
</source>
</source>


Line 500: Line 439:
First we calibrate the cross delay, which is the difference in delay between the two polarizations of the reference antenna, in our case ea01. Remember that the regular delay calibration arbitrarily set the delays for both polarizations of the reference antenna to zero. In reality, there will be a small but noticeable physical path length difference between these two polarizations, for which we try to correct. This calibration uses the same scan on 3C48, assuming that 3C48 is unpolarized at P-band (which is a good approximation).
First we calibrate the cross delay, which is the difference in delay between the two polarizations of the reference antenna, in our case ea01. Remember that the regular delay calibration arbitrarily set the delays for both polarizations of the reference antenna to zero. In reality, there will be a small but noticeable physical path length difference between these two polarizations, for which we try to correct. This calibration uses the same scan on 3C48, assuming that 3C48 is unpolarized at P-band (which is a good approximation).
<source lang='python'>
<source lang='python'>
cross_delay_table_v1 = onespw_ms_path + '.Kc1'
gaincal(vis='TSUB0001_pband.ms', caltable='TSUB0001_pband.Kc1', field='3C48', solint='inf', refant='ea09', minsnr=3.0, gaintype='KCROSS', parang=True,
  gaintable=[ 'TSUB0001_pband.antpos','TSUB0001_pband.rq','TSUB0001_pband.tecim','TSUB0001_pband.K1', 'TSUB0001_pband.B1','TSUB0001_pband.Kc1'],
  interp=['','','','nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag'] )


if os.path.isdir( cross_delay_table_v1 ):
plotcal(caltable='TSUB0001_pband.Kc1', xaxis='antenna', yaxis='delay', markersize=2.0)
  rmtables( cross_delay_table_v1 )
 
gaincal( vis = onespw_ms_path, caltable = cross_delay_table_v1, field = bcal_field, selectdata = True, scan = bcal_scan,
  solint = 'inf', refant = reference_antenna, minsnr = 3.0, gaintype = 'KCROSS', parang = True,
  gaintable = [ delay_table_v1, bandpass_table_v1, gain_table_v1 ], interp = [ 'nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag' ] )
 
plotcal( caltable = cross_delay_table_v1, xaxis = 'antenna', yaxis = 'delay', markersize = 2.0 )
</source>
</source>


Line 517: Line 451:


<source lang='python'>
<source lang='python'>
leakage_table_v1 = onespw_ms_path + '.Df1'
if os.path.isdir( leakage_table_v1 ):
  rmtables( leakage_table_v1 )


polcal( vis = onespw_ms_path, caltable = leakage_table_v1, field = bcal_field, selectdata = True, scan = bcal_scan,
polcal(vis='TSUB0001_pband.ms', caltable='TSUB0001_pband.Df1', field='3C48', solint='inf', refant='ea09', minsnr=3.0, poltype = 'Df',
  solint = 'inf', preavg = 1., refant = reference_antenna, minsnr = 3.0, poltype = 'Df',
   gaintable=[ 'TSUB0001_pband.antpos','TSUB0001_pband.rq','TSUB0001_pband.tecim','TSUB0001_pband.K1', 'TSUB0001_pband.B1','TSUB0001_pband.Kc1'],
   gaintable = [ delay_table_v1, bandpass_table_v1, gain_table_v1, cross_delay_table_v1 ],
   interp=['','','','nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag' ] )
   interp = [ 'nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag' ] )


plotcal( caltable = leakage_table_v1, xaxis = 'freq', yaxis = 'amp', iteration = 'antenna', plotrange = [ 200.,500.,0.,1. ], markersize = 2.0 )
plotcal(caltable='TSUB0001_pband.Df1', xaxis='freq', yaxis='amp', iteration='antenna', plotrange=[ 200.,500.,0.,1. ], markersize=2.0)
</source>
</source>


Line 535: Line 464:
An alternative is to flag based on the real and imaginary parts, which you may or may not want to do.
An alternative is to flag based on the real and imaginary parts, which you may or may not want to do.
<source lang='python'>
<source lang='python'>
plotcal( caltable = leakage_table_v1, xaxis = 'freq', yaxis = 'real', iteration = 'antenna', plotrange = [ 200.,500.,-0.5,0.5 ], markersize = 2.0 )
plotcal(caltable='TSUB0001_pband.Df1', xaxis='freq', yaxis='real', iteration='antenna', plotrange=[ 200.,500.,-0.5,0.5 ], markersize=2.0)


plotcal( caltable = leakage_table_v1, xaxis = 'freq', yaxis = 'imag', iteration = 'antenna', plotrange = [ 200.,500.,-0.5,0.5 ], markersize = 2.0 )
plotcal(caltable='TSUB0001_pband.Df1', xaxis='freq', yaxis='imag', iteration='antenna', plotrange=[ 200.,500.,-0.5,0.5 ], markersize=2.0)
</source>
</source>


Line 547: Line 476:


<source lang='python'>
<source lang='python'>
bandpass_table_v1a = onespw_ms_path + '.B1a'
if os.path.isdir( bandpass_table_v1a ):
  rmtables( bandpass_table_v1a )


bandpass( vis = onespw_ms_path, caltable = bandpass_table_v1a, field = bcal_field, selectdata = True, scan = bcal_scan,
bandpass(vis='TSUB0001_pband.ms', caltable='TSUB0001_pband.B1', field='3C48', solint='inf', refant='ea09', minsnr=3.0, parang = True,
  solint = 'inf', refant = reference_antenna, minsnr = 3.0, solnorm = False, parang = True,
   gaintable=['TSUB0001_pband.antpos','TSUB0001_pband.rq','TSUB0001_pband.tecim','TSUB0001_pband.K1'], interp=['','','','nearest,nearestflag'])
   gaintable = [ delay_table_v1, bandpass_table_v1, gain_table_v1, cross_delay_table_v1, leakage_table_v1 ],
  interp = [ 'nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag' ] )


plotcal( caltable = bandpass_table_v1a, xaxis = 'freq', yaxis = 'amp', iteration = 'antenna', plotrange = [ 200.,500.,0.,2. ], markersize = 2.0 )
plotcal(caltable='TSUB0001_pband.B1' , xaxis='freq', yaxis='amp', iteration='antenna', plotrange=[ 200.,500.,0.,2. ], markersize=2.0)


plotcal( caltable = bandpass_table_v1a, xaxis = 'freq', yaxis = 'phase', iteration = 'antenna', plotrange = [ 200.,500.,-180.,180. ], markersize = 2.0 )
plotcal(caltable='TSUB0001_pband.B1', xaxis='freq', yaxis='phase', iteration='antenna', plotrange=[ 200.,500.,-180.,180. ], markersize=2.0)
</source>
</source>


Line 565: Line 488:


<source lang='python'>
<source lang='python'>
leakage_table_v1a = onespw_ms_path + '.Df1a'


if os.path.isdir( leakage_table_v1a ):
polcal(vis='TSUB0001_pband.ms', caltable='TSUB0001_pband.Df1', field='3C48', solint='inf', refant='ea09', minsnr=3.0, poltype = 'Df',
   rmtables( leakage_table_v1a )
  gaintable=[ 'TSUB0001_pband.antpos','TSUB0001_pband.rq','TSUB0001_pband.tecim','TSUB0001_pband.K1',
              'TSUB0001_pband.B1','TSUB0001_pband.Kc1','TSUB0001_pband.Df1'],
   interp=['','','','nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag' ] )


polcal( vis = onespw_ms_path, caltable = leakage_table_v1a, field = bcal_field, selectdata = True, scan = bcal_scan,
plotcal(caltable='TSUB0001_pband.Df1', xaxis='freq', yaxis='amp', iteration='antenna', plotrange=[ 200.,500.,None,None ], markersize=2.0)
  solint = 'inf', preavg = 1., refant = reference_antenna, minsnr = 3.0, poltype = 'Df',
  gaintable = [ delay_table_v1, bandpass_table_v1, bandpass_table_v1a, gain_table_v1, cross_delay_table_v1, leakage_table_v1 ],
  interp = [ 'nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag' ] )
 
plotcal( caltable = leakage_table_v1a, xaxis = 'freq', yaxis = 'amp', iteration = 'antenna', plotrange = [ 200.,500.,None,None ], markersize = 2.0 )
</source>
</source>


Line 583: Line 502:


<source lang='python'>
<source lang='python'>
applycal( vis = onespw_ms_path, field = '*', selectdata = True, parang = True, applymode = 'calflagstrict', flagbackup = True,
applycal(vis='TSUB0001_pband.ms', parang=True, applymode='calflagstrict', flagbackup=True, gaintable=['TSUB0001_pband.antpos','TSUB0001_pband.rq',
    gaintable = [ delay_table_v1, bandpass_table_v1, bandpass_table_v1a, smoothed_gain_table_v1, cross_delay_table_v1, leak_table_v1, leak_table_v1a ],
            'TSUB0001_pband.tecim','TSUB0001_pband.K1', 'TSUB0001_pband.B1','TSUB0001_pband.Gs1','TSUB0001_pband.Kc1','TSUB0001_pband.Df1'],
    interp = [ 'nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag',
        interp = ['','','','nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag','nearest,nearestflag'])
    'nearest,nearestflag', 'nearest,nearestflag' ] )
</source>
</source>


Line 592: Line 510:


<source lang='python'>
<source lang='python'>
target_name = '0313-192'
split(vis='TSUB0001_pband.ms', outputvis='TSUB0001_pband_target.ms', datacolumn='corrected', field='0313-192')
target_ms_path = target_name + '_PBAND_1SPW.ms'
 
if os.path.isdir( target_ms_path ):
  rmtables( target_ms_path )
 
split( vis = onespw_ms_path, outputvis = target_ms_path, datacolumn = 'corrected', field = target_name, keepflags = False )
</source>
</source>


Line 606: Line 518:


<source lang='python'>
<source lang='python'>
flagdata( vis = pband_ms_path_1spw, mode = 'rflag', field = '*', correlation = 'XY,YX', datacolumn = 'corrected', winsize = 5,
flagdata(vis='TSUB0001_pband_target.ms', mode = 'rflag', field = '*', correlation = 'ABS_XY,ABS_YX', datacolumn = 'corrected', winsize = 5,
   action = 'apply', display = 'report', flagbackup = False )
   action = 'apply', display = 'report', flagbackup = False )
</source>
==Imaging==
At this point, we're ready to make a first image of our target field. Imaging in CASA is done by means of the <i>CLEAN</i> task. The task implements algorithms for wide field imaging such as W-Projection (Cornwell et al. 2008 http://arxiv.org/abs/0807.4161) and Multi Term Multi Frequency Synthesis ( Rau et al. 2011 http://arxiv.org/abs/1106.2745) both of which we will utilize to make out initial image. Even though our source of concern the giant radio galaxy  lies in the centre of our field. The observations were carried out in B configuration of P Band resulting in an effective maximum resolution of approximately 5x5 arcsecs. To effectively model the source spectral index in the sky we will utilize the MTMFS algorithm and use the W-Projection algorithm to make a wide field image (The full beam at P Band is about 3 degrees in diameter). So going by these requirements we can now compute the required cell and image size. To Nyqvist sample the beam we need at least 3 pixels across the synthesized beam width, along the major and minor axes of the point spread function. Consulting the resolution guide of NRAO science page we can see that the expected HPBW for B - Configuration, P -Band Observation is 18.5 arcsec. So we can sample it well by using a 5 arcsec cellsize. Having done that if we decide to make a wide field image to account for all the point sources we will set the image size parameter to 4096 pixels. To enable the wide-field algorithm we set the <i>gridmode='widefield'</i> this invokes the W-Projection algorithm, upon which we set the number of W-Projection planes to be 128. We set the imager mode as Multi Frequency Synthesis using the <i>mode='mfs'</i> parameter along with the number of taylor terms to be considered during imaging to be 2. This allows for the source spectral variation to be modeled by a second order polynomial. We launch the interactive clean process.
<source lang='python'>
clean(vis='TSUB0001_pband_target.ms', imagename='0313-0192_intial_clean', cell=['5.0arcsec','5.0arcsec'], imsize=[4096,4096], mode='mfs',
        nterms=2, gridmode='widefield', wprojplanes=128, stokes='I', niter=10000, spw='3~8,10~15', interactive=True,
        minpb=0.1, usescrarch=T, weighting='briggs', robust=0.0)
</source>
[[File:0313_radio_gal_clean.png|200px|thumb|right]]
The <i>CLEAN</i> command launches an interactive session after a 100 iterations of clean and produces a wide field map with the source at the center and a lot of bright sources far out in the field. As it is a snapshot image the bright sources have significant side lobes and so tight clean boxing helped. This can done in the interactive viewer interface that pops up. This topic has been covered extensively in the 3C391 imaging tutorial which can be found here (https://casaguides.nrao.edu/index.php?title=EVLA_Continuum_Tutorial_3C391#Initial_Imaging). If we proceed with interactive clean with subsequent steps to keep boxing out the strong sources that pop up in the image we finally see the extended emission from the target radio galaxy begin to emerge. Boxing and cleaning to ensure that the residuals of the boxed cleaning look like nice (~10000 clean iterations). We stop the interactive task and look at the final image it produced.
Note the image has some sources still showing strong side lobes and imaging artifacts around them. We expected this as we have carried out phase calibration so far only on our flux calibrator and have just transferred the solutions over to our target field. Since we used the <i>usescratch=T</i> the MODEL_DATA column in the measurement set now contains the initial image model which we will self calibrate against to produce a better image. Also in clean do notice that I the spw's utilized are the cleanest spectral windows that are totally RFI free.
==Self Calibration ==
We now proceed to compute gain phase solutions for our target field using the gain cal task as the first step in self-calibration.
<source lang='python'>
gaincal(vis='TSUB0001_pband_target.ms', caltable='TSUB0001_pband_target.ScG0', field='3C48', solint='inf', refant='ea09',
          spw='3~8,10~15',minsnr=3.0, gaintype='G', parang=True, calmode='p')
applycal(vis='TSUB0001_pband_target.ms', gaintable=['TSUB0001_pband.ScG0'], applymode='calflagstrict')
</source>
Having applied these gain solutions we will once again image the target measurement set which we now expect to have better gain solutions and consequently a better image. We do this by invoking the <i>CLEAN</i> command.
<source lang='python'>
clean(vis='TSUB0001_pband_target.ms', imagename='0313-0192_clean_sc0', cell=['5.0arcsec','5.0arcsec'], imsize=[4096,4096], mode='mfs',
        nterms=2, gridmode='widefield', wprojplanes=128, stokes='I', niter=10000, spw='3~8,10~15', interactive=True,
        minpb=0.1, usescrarch=T, weighting='briggs', robust=0.0)
</source>
</source>
[[File:0313_radio_gal_cselfcal0.png|200px|thumb|right]]


At this point, we're ready to make a first image of our target field, after which we will self-calibrate and re-image.
On boxing and cleaning we already notice that the imaging artifacts has reduced significantly. We also see that the target source appears to contain more structure and a greater amount of flux. Further self calibration iterations involving an amplitude & phase gain calibration, even an on target bandpass and leakage calibration are all possible steps that can be introduced into the self-cal imaging loop. These steps and introducing Multi Scale Clean for the extended source are possible improvements on the current strategy implemented.
[...]


== Appendix: Some P-band data issues you may want to know about ==
== Appendix: Some P-band data issues you may want to know about ==

Latest revision as of 03:29, 12 March 2016

Back to CASA guides

By Preshanth Jagannathan & Huib Intema

Last checked on CASA Version 4.3.1.


Overview

This webpage provides a basic description on how to use CASA to reduce data from the upper part of the new VLA low-band system, also known as P-band, covering roughly 220-480 MHz. The goal is to make a wide-field continuum stokes I image of a typical blank field using the full effective bandwidth.

Obtaining the raw data

For this guide we'll use some test data that was taken in B-configuration. To get a copy of the raw data, go to the NRAO archive query page, and enter the following search parameters (leave the rest on default): 

Telescopes: Jansky VLA (tickbox)
Project Code: TSUB0001
Telescope Config: B (tickbox)
Observing Bands: P (tickbox)

The search should return (at least) 19 rows of results. Enter your valid e-mail address, set the download format to SDM-BDF, check the box on the 17th row (data from Dec. 21, 2013, or "13-Dec-21 03:11:18") at the bottom of the page, and click Get My Data. On the next page, click Retrieve over internet. The next page should report that the data staging is in progress. Wait until you receive an e-mail reporting that your archive data is copied, which should take a few minutes. Download the data onto your computer running CASA.

Starting CASA

Start CASA by typing 

casa

on the (linux) command line. This should start a CASA interactive python (iPython) session, and open a separate log window. The CASA version is reported at startup, both in the python session and the log window:

CASA Version 4.3.1-REL (r32491)

Note that this guide has been written for CASA release 4.3.1, which is the current stable release at NRAO. Due to ongoing development of CASA, this tutorial may or may not work with other versions. You may want to confirm your version before proceeding.

Importing the raw data into CASA

We will begin by importing our data into the measurement sent set format (CASA standard) from the binary format (SDM-BDF) as downloaded from the archive. We do this by means of the importevla task. 

importevla(asdm='TSUB0001.sb28588128.eb28590840.56647.13284232639', vis ='TSUB0001_pband.ms', savecmds=T, outfile='importflags.txt')

This task calls the external tool asdm2MS to perform the conversion. Together with the visibility data, also the flag commands from the VLA online system are imported. These flags are not directly applied, but stored in the importflags.txt file and applied later.

Preliminary data inspection

listobs(vis='TSUB0001_pband.ms', verbose = True, listfile = 'TSUB001_pband_import.listobs')

Alternatively you can also output to the logger by running the command without the listfile argument. The sample listobs output shown has flag verbose=F for the sake of brevity. 

listobs(vis='TSUB0001_pband.ms',verbose=F)

================================================================================
           MeasurementSet Name:  /lustre/pjaganna/evla/P_Band/casa_guide/TSUB0001_pband.ms      MS Version 2
================================================================================
   Observer: Frazer Owen     Project: uid://evla/pdb/1695465  
Observation: EVLA(26 antennas)
Data records: 5382000       Total elapsed time = 2076 seconds
   Observed from   21-Dec-2013/03:11:20.0   to   21-Dec-2013/03:45:56.0 (UTC)

Fields: 3
  ID   Code Name                RA               Decl           Epoch   SrcId      nRows
  0    NONE 3C48D               01:37:41.299431 +33.09.35.13299 J2000   0         884000
  1    NONE 3C48                01:37:41.299431 +33.09.35.13299 J2000   1        2168400
  2    NONE 0313-192            03:15:52.039999 -19.06.44.59999 J2000   2        2329600
Spectral Windows:  (16 unique spectral windows and 1 unique polarization setups)
  SpwID  Name           #Chans   Frame   Ch0(MHz)  ChanWid(kHz)  TotBW(kHz) CtrFreq(MHz) BBC Num  Corrs          
  0      EVLA_P#A0C0#0     128   TOPO     222.000       125.000     16000.0    229.9375       12  RR  RL  LR  LL
  1      EVLA_P#A0C0#1     128   TOPO     238.000       125.000     16000.0    245.9375       12  RR  RL  LR  LL
  2      EVLA_P#A0C0#2     128   TOPO     254.000       125.000     16000.0    261.9375       12  RR  RL  LR  LL
  3      EVLA_P#A0C0#3     128   TOPO     270.000       125.000     16000.0    277.9375       12  RR  RL  LR  LL
  4      EVLA_P#A0C0#4     128   TOPO     286.000       125.000     16000.0    293.9375       12  RR  RL  LR  LL
  5      EVLA_P#A0C0#5     128   TOPO     302.000       125.000     16000.0    309.9375       12  RR  RL  LR  LL
  6      EVLA_P#A0C0#6     128   TOPO     318.000       125.000     16000.0    325.9375       12  RR  RL  LR  LL
  7      EVLA_P#A0C0#7     128   TOPO     334.000       125.000     16000.0    341.9375       12  RR  RL  LR  LL
  8      EVLA_P#A0C0#8     128   TOPO     350.000       125.000     16000.0    357.9375       12  RR  RL  LR  LL
  9      EVLA_P#A0C0#9     128   TOPO     366.000       125.000     16000.0    373.9375       12  RR  RL  LR  LL
  10     EVLA_P#A0C0#10    128   TOPO     382.000       125.000     16000.0    389.9375       12  RR  RL  LR  LL
  11     EVLA_P#A0C0#11    128   TOPO     398.000       125.000     16000.0    405.9375       12  RR  RL  LR  LL
  12     EVLA_P#A0C0#12    128   TOPO     414.000       125.000     16000.0    421.9375       12  RR  RL  LR  LL
  13     EVLA_P#A0C0#13    128   TOPO     430.000       125.000     16000.0    437.9375       12  RR  RL  LR  LL
  14     EVLA_P#A0C0#14    128   TOPO     446.000       125.000     16000.0    453.9375       12  RR  RL  LR  LL
  15     EVLA_P#A0C0#15    128   TOPO     462.000       125.000     16000.0    469.9375       12  RR  RL  LR  LL
Antennas: 26 'name'='station' 
   ID=   0-3: 'ea01'='N32', 'ea02'='N28', 'ea03'='E28', 'ea04'='E24', 
   ID=   4-7: 'ea05'='W08', 'ea07'='N12', 'ea08'='N16', 'ea09'='W28', 
   ID=  8-11: 'ea10'='E04', 'ea11'='W24', 'ea12'='E08', 'ea13'='W20', 
   ID= 12-15: 'ea14'='N24', 'ea15'='E20', 'ea16'='N08', 'ea18'='E36', 
   ID= 16-19: 'ea19'='W16', 'ea20'='N04', 'ea21'='E12', 'ea22'='N20', 
   ID= 20-23: 'ea23'='W36', 'ea24'='W04', 'ea25'='W32', 'ea26'='W12', 
   ID= 24-25: 'ea27'='N36', 'ea28'='E16'

You can inspect the output text file TSUB0001_20131221.listobs in your favorite text editor (e.g., gedit). When taking some time to familiarize yourself with the output format, you will see that:

  • The observation consists of 4 scans, namely [1] 3C48D (hardware setup), [2] 3C48 (primary calibrator), [3] 0313-192 (target), and [4] again 3C48 (primary calibrator).
  • The frequency coverage is 224 - 480 MHz, divided into 16 x 16 MHz spectral windows, each having 128 x 0.125 MHz channels
  • Visibilities are recorded every 2 seconds in full polarization (4 polarization products, (wrongly!) labelled RR,RL,LR,LL).
  • Of the 28 VLA antennas labelled ea01 - ea028, antennas ea06 and ea17 are not participating in this observation.
  • Source, antenna and spectral window IDs start at zero, but scan IDs start at one.


The array configuration can be inspected using: 

plotants(vis='TSUB0001_pband.ms', figfile='TSUB0001_pband_plotants.png')

As noted above, the polarizations of P-band in this data set are wrongly labelled as circular rather than linear. This may seem mostly harmless, but does make a difference for polarization calibration (which we will do later). The following (non-standard) task will check and fix this (and another related problem; see appendix on P-band data issues): 

execfile( '/lustre/hintema/cspam/mytasks.py' )

fixlowband(vis='TSUB0001_pband.ms')

Re-run the listobs() task to check if the fix was correctly applied.

Addendum: If you're not working on the NRAO DSOC network, you won't have access to the lustre file system. In that case, download this tarball in your CASA working directory, open a shell there, and type the following commands: 

tar -xzvf casa_vla_lowband.tar.gz

buildmytasks

Now go back to your CASA session and type: 

execfile('mytasks.py')

fixlowband(vis='TSUB0001_pband.ms')

The re-run listobs output will show the correct polarization. 

================================================================================
           MeasurementSet Name:  /lustre/pjaganna/evla/P_Band/casa_guide/TSUB0001_pband.ms      MS Version 2
================================================================================
   Observer: Frazer Owen     Project: uid://evla/pdb/1695465
Observation: EVLA(26 antennas)
Data records: 5382000       Total elapsed time = 2076 seconds
   Observed from   21-Dec-2013/03:11:20.0   to   21-Dec-2013/03:45:56.0 (UTC)

Fields: 3
  ID   Code Name                RA               Decl           Epoch   SrcId      nRows
  0    NONE 3C48D               01:37:41.299431 +33.09.35.13299 J2000   0         884000
  1    NONE 3C48                01:37:41.299431 +33.09.35.13299 J2000   1        2168400
  2    NONE 0313-192            03:15:52.039999 -19.06.44.59999 J2000   2        2329600
Spectral Windows:  (16 unique spectral windows and 1 unique polarization setups)
  SpwID  Name           #Chans   Frame   Ch0(MHz)  ChanWid(kHz)  TotBW(kHz) CtrFreq(MHz) BBC Num  Corrs
  0      EVLA_P#A0C0#0     128   TOPO     222.000       125.000     16000.0    229.9375       12  XX  XY  YX  YY
  1      EVLA_P#A0C0#1     128   TOPO     238.000       125.000     16000.0    245.9375       12  XX  XY  YX  YY
  2      EVLA_P#A0C0#2     128   TOPO     254.000       125.000     16000.0    261.9375       12  XX  XY  YX  YY
  3      EVLA_P#A0C0#3     128   TOPO     270.000       125.000     16000.0    277.9375       12  XX  XY  YX  YY
  4      EVLA_P#A0C0#4     128   TOPO     286.000       125.000     16000.0    293.9375       12  XX  XY  YX  YY
  5      EVLA_P#A0C0#5     128   TOPO     302.000       125.000     16000.0    309.9375       12  XX  XY  YX  YY
  6      EVLA_P#A0C0#6     128   TOPO     318.000       125.000     16000.0    325.9375       12  XX  XY  YX  YY
  7      EVLA_P#A0C0#7     128   TOPO     334.000       125.000     16000.0    341.9375       12  XX  XY  YX  YY
  8      EVLA_P#A0C0#8     128   TOPO     350.000       125.000     16000.0    357.9375       12  XX  XY  YX  YY
  9      EVLA_P#A0C0#9     128   TOPO     366.000       125.000     16000.0    373.9375       12  XX  XY  YX  YY
  10     EVLA_P#A0C0#10    128   TOPO     382.000       125.000     16000.0    389.9375       12  XX  XY  YX  YY
  11     EVLA_P#A0C0#11    128   TOPO     398.000       125.000     16000.0    405.9375       12  XX  XY  YX  YY
  12     EVLA_P#A0C0#12    128   TOPO     414.000       125.000     16000.0    421.9375       12  XX  XY  YX  YY
  13     EVLA_P#A0C0#13    128   TOPO     430.000       125.000     16000.0    437.9375       12  XX  XY  YX  YY
  14     EVLA_P#A0C0#14    128   TOPO     446.000       125.000     16000.0    453.9375       12  XX  XY  YX  YY
  15     EVLA_P#A0C0#15    128   TOPO     462.000       125.000     16000.0    469.9375       12  XX  XY  YX  YY
Antennas: 26 'name'='station'
   ID=   0-3: 'ea01'='N32', 'ea02'='N28', 'ea03'='E28', 'ea04'='E24',
   ID=   4-7: 'ea05'='W08', 'ea07'='N12', 'ea08'='N16', 'ea09'='W28',
   ID=  8-11: 'ea10'='E04', 'ea11'='W24', 'ea12'='E08', 'ea13'='W20',
   ID= 12-15: 'ea14'='N24', 'ea15'='E20', 'ea16'='N08', 'ea18'='E36',
   ID= 16-19: 'ea19'='W16', 'ea20'='N04', 'ea21'='E12', 'ea22'='N20',
   ID= 20-23: 'ea23'='W36', 'ea24'='W04', 'ea25'='W32', 'ea26'='W12',
   ID= 24-25: 'ea27'='N36', 'ea28'='E16'

Initial processing steps - Antenna Position, Requantizer Gain, Ionospheric Correction

When importing the data, we saved the flags from the online system to an ASCII text file. This gives us the opportunity to review the flag commands before applying them. For our data set, the flag file is called TSUB0001_20131221.flags. Please load this file into your favorite text editor. The bulk of the flag commands refer to times when the VLA is slewing (ANTENNA_NOT_ON_SOURCE) or when the movable secondary reflector of the VLA's Cassegrain system is not in place (SUBREFLECTOR_ERROR). The latter can cause antenna gain variations (amplitude and phase), so it is safest to apply all the flags. The final two rows in the flag file will remove visibilities that are pure zero, and flag antennas that are partly blocked by other antennas (shadowing; this occurs mostly in compact configurations when observing along a VLA arm). 

flagdata(vis='TSUB0001_pband.ms', mode='list', inpfile='importflags.txt', action='apply', reason='any', flagbackup=True)

Antenna Position Corrections

Now that we have applied the import flags we will correct for the antenna position offsets. Antenna positional errors translates to an error in the measured visibilities and need to be accounted for before we proceed with any of the other calibration steps. 

gencal(vis='TSUB0001_pband.ms',caltable='TSUB0001_pband.antpos',caltype='antpos')

The output shows that there are five antennas with positional offsets, this is the kind 


2015-06-17 05:07:03 INFO gencal	offsets for antenna ea03 : -0.00240  -0.00410   0.00340
2015-06-17 05:07:03 INFO gencal	offsets for antenna ea04 : -0.00100  -0.00110   0.00180
2015-06-17 05:07:03 INFO gencal	offsets for antenna ea05 : -0.00100   0.00100   0.00200
2015-06-17 05:07:03 INFO gencal	offsets for antenna ea15 : -0.00190  -0.00240   0.00180
2015-06-17 05:07:03 INFO gencal	offsets for antenna ea18 : -0.00180  -0.00590   0.00500

Ionospheric TEC Corrections

Low frequency observations are affected by the ionosphere. A delay in the signal path is introduced between the two polarization of light that varies both as a function of time and line of sight (direction dependent). The delay is proportional to the Total Electron Content (TEC) along the line of sight and is inversely proportional to the square of the frequency.GPS measurements at two different frequencies provides us with an estimate of the TEC per square metre. This correction has been implemented in CASA which we shall apply as a calibration table by means of the gencal task. The task requires a TEC map that we will generate utilizing casa recipes. 

from recipes import tec_maps

tec_image, tec_rms_image = tec_maps.create(vis='TSUB0001_pband.ms',doplot=True)

gencal(vis='TSUB0001_pband.ms',caltable='TSUB0001_pband.tecim',caltype='tecim',infile=tec_image)

A word of caution regarding the TEC map generation. The IGS website updates measurements only two weeks after the date of observation.

Requantizer Gains

The next step will correct the visibility amplitudes for the signal leveling that occurs at the inputs of the WIDAR correlator, the so-called requantizer gains. These levels (per antenna, per polarization, per spectral window) are stored with the measurement set (in the SYSPOWER sub-table). Currently, this step is not essential, since the levels get set only once at the start of an observation, and bandpass calibration will correct for this. But it will make your bandpass plots look better if you have multiple spectral windows. And, more importantly, in the future it may be that the levels will become time-variable during an observation, so we'd better be prepared.

The correction is done by means of the gencal task, where information in the SYSPOWER sub-table gets translated into a gain table. 

gencal(vis='TSUB0001_pband.ms/',caltype='rq',caltable='TSUB0001_pband.rq')


NOTE: The experimental task mstransform() combines split() and hanningsmooth() and several other data manipulation operations, but its output is not yet to be trusted for all operations.

Dead and swapped antennas

Now it is time to have a first visual look at the uncalibrated visibility data. It is important to identify dead antennas / polarizations early on, so we can exclude them from further data processing (which is the most efficient data 'reduction' :-) ). A convenient way of doing this is through the plotms task on bright calibrator(s) (again, see the listobs() output), in our case 3C48.

plotms(vis='TSUB0001_pband.ms',xaxis='freq',yaxis='amp',antenna='ea01',correlation='XX,YY', field='3C48', 
           plotrange=[0.2,0.5,0.0,100.0], coloraxis='spw',xlabel='Frequency',ylabel='Amplitude',iteraxis='baseline', 
           plotfile='TSUB0001_pband_3C48_prebp.png')

This will load the plotms window that has queued up all the baselines of antenna ea01 and is currently displaying the first amplitude vs frequency plot with the colors representing the different spectral windows. Note that we are only plotting the XX and YY correlation of the antennas. This is because we expect most of the power in the linear cross correlation products to be in these two correlations. Notice that some spectral windows in particular are badly affected due to rfi. This is the issue that we will deal with after identifying dead antennas. If you scroll through the baselines using the green forward buttons in the plotms you will notice that the amplitude is particularly low when you encounter ea10 and ea13. If we make a quick plot of the baseline made by ea10 and ea13 we can see that there is no power in any of the polarizations so these two antennas were dead and can be flagged out. We also notice that ea25 has low XX polarization but also notice the excess power in the XY polarization. This tells us that the polarization inputs have been wrongly labeled we have a local task that we built already and we will utilize that to swap the polarization labels for antenna ea25. In addition to the dead antennas we also have a dummy scan during slew on a source named 3C48D, we shall flag that scan by hand too. 

flagdata(vis='TSUB0001_pband.ms',mode='manual',antenna='ea10,ea13')
flagdata(vis='TSUB0001_pband.ms',mode='manual',field='3C48D')
swappol(vis='TSUB0001_pband.ms',antenna='ea25')

NOTE: swappol() is not part of the default CASA release, but was imported together with fixlowband(). It is a pure python implementation operating on visibilities, therefore it may take a while to finish. The way in which we run swappol() we overwrite the previous visibilities, so no new measurement set is created.

We also noticed sharp rfi peaks during our data examination, to prevent Gibbs ringing it is best to hanning smooth the data at this juncture before we proceed further with automatic flagging. 

hanningsmooth(vis='TSUB0001_pband.ms',datacolumn='data')

We should replot the primary calibrator to see the effect of hanning smoothing on the data and the rfi. 

plotms(vis='TSUB0001_pband.ms',xaxis='freq',yaxis='amp',antenna='ea01',correlation='XX,YY',
           field='3C48', plotrange=[0.2,0.5,0.0,100.0], coloraxis='spw',xlabel='Frequency',ylabel='Amplitude',
           iteraxis='baseline', plotfile='TSUB0001_pband_3C48_prebp_hanning.png')

Automatic flagging

Before throwing away bad spectral windows by hand, we'll let the automated flaggers in CASA have a go at it. The main task for this is flagdata(). While working on your data, flagdata() produces an abundance of output, not all of which I find easy to understand. To help with this, we can get the flag status of our data before and after auto-flagging by running flagdata() in the summary mode. This first call provides the flag statistics of the flagging up to now, before any auto-flagging: 

summary_1 = flagdata(vis='TSUB0001_pband.ms', mode='summary')

This returns a python dictionary with flagged versus total visibilities along various axes (antenna, scan, spw, field, correlation, etc.). For example, if we want to know the percentage flagged per scan, run the following (note that the scans may not appear in sorted order): 

axis = 'scan'
for id, stats in summary_1[ axis ].iteritems():
  print '%s %s: %5.1f percent flagged' % ( axis, id, 100. * stats[ 'flagged' ] / stats[ 'total' ] )

For this example, you will notice that scan 1 (the dummy scan on 3C48D) is 100 percent flagged, which is what we did during the initial flagging.

We will run flagdata() in the tfcrop mode, which will (per scan, per baseline, per spectral window, per polarization) look for visibility amplitude outliers. It uses a 5-piece polynomial in an attempt to remove any intrinsic bandpass (amplitude) structure (we did not calibrate for bandpass yet). We run the task twice to allow for slightly deeper flagging. For the first run, we tell CASA to make a backup of our visibility flag status up to here, giving us an option to restore them if we choose over-aggresive flagging parameters and consequently over-flag our data. The flag backup file name is reported in the log window, and can be found in the TSUB0001_20131221_PBAND.ms.flagversions directory. 

flagdata(vis='TSUB0001_pband.ms', field='*', mode='tfcrop', datacolumn='data', timecutoff=4., freqcutoff=3., maxnpieces=5,
  action='apply', display='report', flagbackup=True, combinescans=True, ntime='3600s', correlation='ABS_XY,ABS_YX')

flagdata(vis='TSUB0001_pband.ms', field='*', mode='tfcrop', datacolumn='data', timecutoff=3., freqcutoff=3., maxnpieces=2,
  action='apply', display='report', flagbackup=True, combinescans=True, ntime='3600s', correlation='ABS_XY,ABS_YX')

flagdata(vis='TSUB0001_pband.ms', mode='extend')

Let's get another flag summary to see how much extra data got flagged: 

summary_2 = flagdata(vis='TSUB0001_pband.ms' , mode='summary')

axis = 'scan'
for value, stats in summary_2[ axis ].iteritems():
  old_stats = summary_1[ axis ][ value ]
  print '%s %s: %5.1f percent flagged additionally' % ( axis, value, 100. * ( stats[ 'flagged' ] - old_stats[ 'flagged' ] ) / stats[ 'total' ] )

plotms(vis='TSUB0001_pband.ms',xaxis='freq',yaxis='amp',antenna='ea01',correlation='XX,YY',
           field='3C48', plotrange=[0.2,0.5,0.0,100.0], coloraxis='spw',xlabel='Frequency',ylabel='Amplitude',
           iteraxis='baseline', plotfile='TSUB0001_pband_3C48_prebp_hanning.png')

There is quite a bit of RFI that gets missed in RFI-rich spectral windows. The solution to this is to provide the flagging routines with more contrast between healthy and affected data. For this, we will perform a preliminary bandpass calibration to take out the bandpass shape. We will do some coarse preliminary calibration apply it to the calibrator before flagging for RFI once more. 

gaincal(vis='TSUB0001_pband.ms', caltable='TSUB0001_pband.G0', gaintype='G', calmode='p', solint='int', field='3C48',refant='ea09',
  gaintable=['TSUB0001_pband.antpos','TSUB0001_pband.rq','TSUB0001_pband.tecim'])

gaincal(vis='TSUB0001_pband.ms', caltable='TSUB0001_pband.K0', gaintype='K', solint='inf', field='3C48',refant='ea09',
  gaintable=['TSUB0001_pband.antpos','TSUB0001_pband.rq','TSUB0001_pband.tecim','TSUB0001_pband.G0'])

bandpass(vis='TSUB0001_pband.ms', caltable='TSUB0001_pband.bp0', solint='inf', field='3C48',refant='ea09', minsnr=2.0,
  gaintable=['TSUB0001_pband.antpos','TSUB0001_pband.rq','TSUB0001_pband.tecim','TSUB0001_pband.G0','TSUB0001_pband.K0'])

applycal(vis='TSUB0001_pband.ms', field='*', applymode='calflagstrict',
  gaintable=['TSUB0001_pband.antpos','TSUB0001_pband.rq','TSUB0001_pband.tecim','TSUB0001_pband.G0','TSUB0001_pband.K0'] )

We will now flag the corrected data column that contains the coarsely calibrated visibilities which provide better contrast to the flagging algorithms to remove the RFI present. We will begin by running the RFLAG algorithm in the task flagdata. We will yet again in summary mode see how much more of the observed data was flagged. 

flagdata(vis='TSUB0001_pband.ms', field='*', mode='rflag', datacolumn='corrected', timedevscale=4., freqdevscale=3.,
  action='apply',  flagbackup=True, combinescans=True, ntime='3600s', maxnpieces=5)

flagdata(vis='TSUB0001_pband.ms', field='*', mode='rflag', datacolumn='corrected', timedevscale=4., freqdevscale=3.,
  action='apply',  flagbackup=True, combinescans=True, ntime='3600s', maxnpieces=5)
summary_3 = flagdata(vis='TSUB0001_pband.ms' , mode='summary')

axis = 'scan'
for value, stats in summary_2[ axis ].iteritems():
  old_stats = summary_1[ axis ][ value ]
  print '%s %s: %5.1f percent flagged additionally' % ( axis, value, 100. * ( stats[ 'flagged' ] - old_stats[ 'flagged' ] ) / stats[ 'total' ] )

plotms(vis='TSUB0001_pband.ms',xaxis='freq',yaxis='amp',antenna='ea01',correlation='XX,YY',
           field='3C48', plotrange=[0.2,0.5,0.0,100.0], coloraxis='spw',xlabel='Frequency',ylabel='Amplitude',
           iteraxis='baseline', plotfile='3C48_postbandpass_rflag.png', ydatacolumn='corrected')

Taking another visual look at the data using plotms we examine the effect of the automated flagging on the calibrator. We find that most spectral windows are RFI free and look great but for spectral windows 2, 9, 13, 15. We dig deeper into the cause of the RFI in these instance by using the locate tool in the plotms tool. It comes to our notice that most of the RFI in spw 13 comes from antenna ea04. So we can flag ea04 for spa 13. This still leaves 3 problematic spw's 2,9,15. These spectral windows have been designed to catch most of the RFI in the P-Band and in the event that they are terribly contaminated we can ignore these spw's and proceed with the actual calibration. Nonetheless a final pass of flagdata in RFLAG mode over the contaminated spw's does not hurt. 

flagdata(vis='TSUB0001_pband.ms', field='3C48', mode='rflag', datacolumn='corrected', timedevscale=4., freqdevscale=3.,
  action='apply',  flagbackup=True, combinescans=True, ntime='3600s', maxnpieces=5, spw='2,9,15')

Flux Density Calibration

The initial flagging has cleaned up most of the stray RFI across the band. This allows us to proceed with the actual calibration of the data. Before we get to the calibration tables its essential to do the flux density calibration of our calibrator 3C48. This is done by means of the setjy task. Before we run the task on hand we first clear the preliminary calibration that was carried out to enable better flagging. We do that but running the clearcal task in casa. 

clearcal(vis='TSUB0001_pband.ms')

setjy(vis='TSUB0001_pband.ms', standard='Scaife-Heald 2012', field='3C48')

Delay and bandpass calibration

First we use the full bandwidth on a single scan on the primary calibrator to determine a single delay per antenna per polarization. This will determine a single (approximate) phase slope across frequency, mainly caused by propagation effects in the (time-variable) ionosphere, cable length differences and electronics in the signal paths from antenna feeds to correlator. Note that for a short (5-10 minutes), single scan on the calibrator we can get away with solving for a time-invariant delay per antenna per polarization. We will again use scan 4, which is the longest scan on 3C48. 

gaincal(vis='TSUB0001_pband.ms', caltable='TSUB0001_pband.K1', field='3C48', solint='inf', refant='ea09', gaintype='K', parang=True)

plotcal(caltable='TSUB0001_pband.K1', xaxis='antenna', yaxis='delay', markersize=2.0)

Since we picked ea09 as our reference antenna, the delays for this antenna ID (=0) are arbitrarily set to zero. All other delays should be within 30 nanosec or so (which they should be for this data set). Larger values should be treated with suspicion, probably indicating a problem with the antenna / polarization. We will use the bandpass calibration (next) to verify this.

Next is to determine the bandpass calibration. We use the same source and scan, and apply the delay calibration before solving for the bandpass. Note that we request a minimum SNR of 3, which will make the solve fail for the worst channels, but not all. 

bandpass(vis='TSUB0001_pband.ms', caltable='TSUB0001_pband.B1', field='3C48', solint='inf', refant='ea09', minsnr=3.0, parang = True,
  gaintable=['TSUB0001_pband.K1'], interp=['nearest,nearestflag'])

After solving, we inspect the bandpass calibration amplitudes and phases, and flag any obvious residual outliers. Outliers may be found above or below the average bandpass curves, and tend to arise in the same channel (frequency) ranges for all antennas. Don't spend more than 1 minute per antenna. And don't worry about missing some bad points; these bad channels will likely get flagged in a later stage anyway. Note that the bandpass amplitudes and phases may differ between polarizations of the same antenna. Also note that the bandpass phases across frequency should not have an overall gradient (this was removed by the delay calibration). Also also note that bandpass phases may wrap around from +/-180 to -/+180 degrees. If that is problematic, re-run the phase plotcal() with plotrange = [ 200.,500.,0.,360. ]. Also note that there are no bandpass solutions for antennas ea10, ea13 and ea25. 

plotcal(caltable='TSUB0001_pband.B1', xaxis='freq', yaxis='amp', iteration='antenna', markersize=2.0)

plotcal(caltable='TSUB0001_pband.B1', xaxis='freq', yaxis='phase', iteration='antenna', plotrange=[ 200.,500.,-180.,180. ], markersize = 2.0)

Flagging is made easier by maximizing the plotcal() window. We don't set the scale for plotting the bandpass amplitudes, which is convenient when flagging outliers as it will auto-rescale. While plotting the bandpass phases we do fix the scale, as it provides a better feeling for the relevant magnitude of phase outliers. Flagging of outliers is done using the Mark Region and Flag buttons. We are using this mechanism to efficiently flag data manually, since it is done per antenna rather than per baseline. These 'flags' will not be permanent until we call applycal() later on, but are effective while applying the calibration on the fly, as we will do in subsequent calibration steps.

Gain calibration

With the delay and bandpass calibration in place, we will now look more closely at the time-variable behaviour of the VLA. We use the same scan on our primary calibrator 3C48 over the full effective bandwidth to determine gain calibrations: one complex value per antenna per polarization per integration time. The delay and bandpass tables are applied on the fly. 

gaincal(vis='TSUB0001_pband.ms', caltable='TSUB0001_pband.G1', field='3C48', solint = 'int', refant = reference_antenna, minsnr = 3.0, gaintype = 'G',
calmode = 'ap', gaintable = [ delay_table_v1, bandpass_table_v1 ], interp = ['','','','nearest,nearestflag', 'nearest,nearestflag' ], parang = True )

Similar to the bandpass flagging, we will flag outlier gain solutions in both amplitude and phase. The regular dips in the gain amplitudes every 4 minutes is caused by a very strong RFI signal from the local VLA site radio; this issue is being addressed. For the amplitudes, the final scatter around an average close to one should be a few percent. Note that antenna ea12 has quite a large bump in the gain amplitudes on one polarization. For this and other cases, remember that the gain solutions for both polarizations can be flagged together, as we will loose the equivalent data of the surviving polarization anyway. 

plotcal(caltable='TSUB0001_pband.ms', xaxis='time', yaxis='amp', iteration='antenna', markersize=2.0 )

plotcal(caltable='TSUB0001_pband.ms', xaxis='time', yaxis='phase', iteration='antenna', plotrange=[ None,None,-180.,180.], markersize=2.0)

Note that when plotting gains against time, it is convenient to use auto-scaling of the time axis by putting None in the first two fields of plotrange.

Next, we will prepare a smoothed and interpolated version of the gain calibration table, which will be applied later on to the target field data. This prevents flagging of target field data when one of the edges of a calibrator scan is flagged for one or more antennas. The result can be inspected with plotcal. 

smoothcal(vis='TSUB0001_pband.ms', tablein='TSUB0001_pband.G1', caltable='TSUB0001_pband.Gs1', smoothtype='median', smoothtime = 60.*60.)

plotcal(caltable='TSUB0001_pband.Gs1', xaxis='time', yaxis='amp', iteration='antenna', plotrange=[ None,None,0.8,1.2 ], markersize=2.0)

plotcal(caltable='TSUB0001_pband.Gs1', xaxis='time', yaxis='phase', iteration='antenna', plotrange=[ None,None,-10.,10. ], markersize=2.0)

Instrumental polarization calibration

The VLA at P-band measures the incoming radio waves using dual-dipole feeds mounted near the primary focus of the dishes. While much effort is put into mechanically making these dual-dipoles orthogonal, and to strongly suppress any cross-talk between the electronic signal paths of both dipoles, there is some (frequency-dependent) polarization leakage between them. Furthermore, the P-band feeds are mounted and aligned on each antenna by hand, allowing for 5-10 degrees difference in the orientation of the dipoles between antennas. For both intensity and polarization imaging, it is important to calibrate the polarization products in our data to optimize the (i) orthogonality between X and Y for each antenna and (ii) the alignment between Xs and Ys of different antennas. Note that this does not (yet) include calibration of the absolute polarization angle, for which we need a linearly polarized source with known polarization angle.

First we calibrate the cross delay, which is the difference in delay between the two polarizations of the reference antenna, in our case ea01. Remember that the regular delay calibration arbitrarily set the delays for both polarizations of the reference antenna to zero. In reality, there will be a small but noticeable physical path length difference between these two polarizations, for which we try to correct. This calibration uses the same scan on 3C48, assuming that 3C48 is unpolarized at P-band (which is a good approximation). 

gaincal(vis='TSUB0001_pband.ms', caltable='TSUB0001_pband.Kc1', field='3C48', solint='inf', refant='ea09', minsnr=3.0, gaintype='KCROSS', parang=True,
  gaintable=[ 'TSUB0001_pband.antpos','TSUB0001_pband.rq','TSUB0001_pband.tecim','TSUB0001_pband.K1', 'TSUB0001_pband.B1','TSUB0001_pband.Kc1'], 
  interp=['','','','nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag'] )

plotcal(caltable='TSUB0001_pband.Kc1', xaxis='antenna', yaxis='delay', markersize=2.0)

The plot shows a small ~2 nanosec delay between X and Y of the reference antenna. Note that the table is generated in such a way that this correction gets applied to all antennas.

Next, we determine the leakage of the X signal into Y, and Y into X. For a well-behaved instrument, the magnitude of the leakage (the so-called D-terms) should be much less than one. Since the leakage is frequency dependent, we solve for leakage for each frequency channel, somewhat similar to bandpass calibration. The way to do this in CASA is through the polcal() task with the option poltype = 'Df'. Again, we can inspect (and edit) the leakage calibration table with plotcal() 

polcal(vis='TSUB0001_pband.ms', caltable='TSUB0001_pband.Df1', field='3C48', solint='inf', refant='ea09', minsnr=3.0, poltype = 'Df',
  gaintable=[ 'TSUB0001_pband.antpos','TSUB0001_pband.rq','TSUB0001_pband.tecim','TSUB0001_pband.K1', 'TSUB0001_pband.B1','TSUB0001_pband.Kc1'],
  interp=['','','','nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag' ] )

plotcal(caltable='TSUB0001_pband.Df1', xaxis='freq', yaxis='amp', iteration='antenna', plotrange=[ 200.,500.,0.,1. ], markersize=2.0)

Before starting any flagging, it is best to first cycle past all antennas to get a feel for the magnitude of the leakage terms. For most antennas, both leakage terms are roughly equal in amplitude across the band, with values between 0 and 0.15. Few have higher leakages (up to 0.3). The most likely explanation is a larger misalignment of the P-band feed of these antennas relative to the others. Examples are ea05, ea11, and ea24. These antennas do not need to be flagged, but are merely expected to benefit most from this polarization calibration. Some examples of suspicious data that is probably best to flag: antenna ea03 has a peculiar bump in one of the leakage terms, ea12 has an upturn at the low-frequency end of the band, ea14 has a peak near the low-frequency end of the band.

It is better not to flag based on phases, because the closer the amplitude of the leakage terms get to zero, the more erratic the phase behaviour will be. An alternative is to flag based on the real and imaginary parts, which you may or may not want to do. 

plotcal(caltable='TSUB0001_pband.Df1', xaxis='freq', yaxis='real', iteration='antenna', plotrange=[ 200.,500.,-0.5,0.5 ], markersize=2.0)

plotcal(caltable='TSUB0001_pband.Df1', xaxis='freq', yaxis='imag', iteration='antenna', plotrange=[ 200.,500.,-0.5,0.5 ], markersize=2.0)

Re-iterate bandpass and leakage calibration

The largest instrumental effects on our data are now captured in several relevant calibration tables. But since we have used uncalibrated / partially calibrated data to determine these calibration tables, there are improvements to be obtained when re-determining some calibration tables. In this guide we re-calibrate the bandpass and polarization leakage since these are expected to have changed most. Our re-calibrations will be relative to the existing bandpass and leakage calibrations, so that we can more easily see deviations from the ideal case.

First we update the bandpass calibration by determining a new, relative bandpass table while applying all previous calibrations. In the ideal case that our data and previous calibrations were perfect, the new bandpass table would contain gain amplitudes and phases of one and zero, respectively, for all (unflagged) antennas, polarizations and channels. 

bandpass(vis='TSUB0001_pband.ms', caltable='TSUB0001_pband.B1', field='3C48', solint='inf', refant='ea09', minsnr=3.0, parang = True,
  gaintable=['TSUB0001_pband.antpos','TSUB0001_pband.rq','TSUB0001_pband.tecim','TSUB0001_pband.K1'], interp=['','','','nearest,nearestflag'])

plotcal(caltable='TSUB0001_pband.B1' , xaxis='freq', yaxis='amp', iteration='antenna', plotrange=[ 200.,500.,0.,2. ], markersize=2.0)

plotcal(caltable='TSUB0001_pband.B1', xaxis='freq', yaxis='phase', iteration='antenna', plotrange=[ 200.,500.,-180.,180. ], markersize=2.0)

In the plotcal() calls above, the vertical axis scales have been chosen to indicate the relevance of the ideal case of unity amplitudes and zero phases. For most antennas, the bandpass solutions are close to ideal. Some antennas show some noticeable structure across frequency. A general guideline for flagging (or not) is that we expect the relative bandpass to be close to one in amplitude and vary smoothly in frequency. 

polcal(vis='TSUB0001_pband.ms', caltable='TSUB0001_pband.Df1', field='3C48', solint='inf', refant='ea09', minsnr=3.0, poltype = 'Df',
  gaintable=[ 'TSUB0001_pband.antpos','TSUB0001_pband.rq','TSUB0001_pband.tecim','TSUB0001_pband.K1',
              'TSUB0001_pband.B1','TSUB0001_pband.Kc1','TSUB0001_pband.Df1'],
  interp=['','','','nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag' ] )

plotcal(caltable='TSUB0001_pband.Df1', xaxis='freq', yaxis='amp', iteration='antenna', plotrange=[ 200.,500.,None,None ], markersize=2.0)

Transfer of calibrations to the target field

So far, we have used one scan on the primary calibrator 3C48 to derive various calibration tables. These calibrations will now be applied to all scans of all sources in our measurement set, which includes 3C48 itself and our target field 0313-192. As described above, the task applycal() creates a CORRECTED_DATA column in which the calibrated visibilities get stored. This may take a while. 

applycal(vis='TSUB0001_pband.ms', parang=True, applymode='calflagstrict', flagbackup=True, gaintable=['TSUB0001_pband.antpos','TSUB0001_pband.rq',
             'TSUB0001_pband.tecim','TSUB0001_pband.K1', 'TSUB0001_pband.B1','TSUB0001_pband.Gs1','TSUB0001_pband.Kc1','TSUB0001_pband.Df1'],
         interp = ['','','','nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag', 'nearest,nearestflag','nearest,nearestflag'])

Now we split off the calibrated target field data, meaning that the visibilities of source 0313-192 get copied from the CORRECTED_DATA column to the DATA column of a new measurement set. This is convenient for further processing. 

split(vis='TSUB0001_pband.ms', outputvis='TSUB0001_pband_target.ms', datacolumn='corrected', field='0313-192')

Cross-polarization flagging

All flagging up to this point was done on the calibrator 3C48, and some of these flags got propagated to the target field data. we will now extend the flagging of the target field data by flagging based on the cross-polarization visibilities XY and YX. Most RFI is strongly polarized, and will stand out clearly against the typically very weakly polarized target field. 

flagdata(vis='TSUB0001_pband_target.ms', mode = 'rflag', field = '*', correlation = 'ABS_XY,ABS_YX', datacolumn = 'corrected', winsize = 5,
  action = 'apply', display = 'report', flagbackup = False )


Imaging

At this point, we're ready to make a first image of our target field. Imaging in CASA is done by means of the CLEAN task. The task implements algorithms for wide field imaging such as W-Projection (Cornwell et al. 2008 http://arxiv.org/abs/0807.4161) and Multi Term Multi Frequency Synthesis ( Rau et al. 2011 http://arxiv.org/abs/1106.2745) both of which we will utilize to make out initial image. Even though our source of concern the giant radio galaxy lies in the centre of our field. The observations were carried out in B configuration of P Band resulting in an effective maximum resolution of approximately 5x5 arcsecs. To effectively model the source spectral index in the sky we will utilize the MTMFS algorithm and use the W-Projection algorithm to make a wide field image (The full beam at P Band is about 3 degrees in diameter). So going by these requirements we can now compute the required cell and image size. To Nyqvist sample the beam we need at least 3 pixels across the synthesized beam width, along the major and minor axes of the point spread function. Consulting the resolution guide of NRAO science page we can see that the expected HPBW for B - Configuration, P -Band Observation is 18.5 arcsec. So we can sample it well by using a 5 arcsec cellsize. Having done that if we decide to make a wide field image to account for all the point sources we will set the image size parameter to 4096 pixels. To enable the wide-field algorithm we set the gridmode='widefield' this invokes the W-Projection algorithm, upon which we set the number of W-Projection planes to be 128. We set the imager mode as Multi Frequency Synthesis using the mode='mfs' parameter along with the number of taylor terms to be considered during imaging to be 2. This allows for the source spectral variation to be modeled by a second order polynomial. We launch the interactive clean process. 

clean(vis='TSUB0001_pband_target.ms', imagename='0313-0192_intial_clean', cell=['5.0arcsec','5.0arcsec'], imsize=[4096,4096], mode='mfs',
         nterms=2, gridmode='widefield', wprojplanes=128, stokes='I', niter=10000, spw='3~8,10~15', interactive=True, 
         minpb=0.1, usescrarch=T, weighting='briggs', robust=0.0)

The CLEAN command launches an interactive session after a 100 iterations of clean and produces a wide field map with the source at the center and a lot of bright sources far out in the field. As it is a snapshot image the bright sources have significant side lobes and so tight clean boxing helped. This can done in the interactive viewer interface that pops up. This topic has been covered extensively in the 3C391 imaging tutorial which can be found here (https://casaguides.nrao.edu/index.php?title=EVLA_Continuum_Tutorial_3C391#Initial_Imaging). If we proceed with interactive clean with subsequent steps to keep boxing out the strong sources that pop up in the image we finally see the extended emission from the target radio galaxy begin to emerge. Boxing and cleaning to ensure that the residuals of the boxed cleaning look like nice (~10000 clean iterations). We stop the interactive task and look at the final image it produced.

Note the image has some sources still showing strong side lobes and imaging artifacts around them. We expected this as we have carried out phase calibration so far only on our flux calibrator and have just transferred the solutions over to our target field. Since we used the usescratch=T the MODEL_DATA column in the measurement set now contains the initial image model which we will self calibrate against to produce a better image. Also in clean do notice that I the spw's utilized are the cleanest spectral windows that are totally RFI free.

Self Calibration

We now proceed to compute gain phase solutions for our target field using the gain cal task as the first step in self-calibration.

gaincal(vis='TSUB0001_pband_target.ms', caltable='TSUB0001_pband_target.ScG0', field='3C48', solint='inf', refant='ea09', 
           spw='3~8,10~15',minsnr=3.0, gaintype='G', parang=True, calmode='p')

applycal(vis='TSUB0001_pband_target.ms', gaintable=['TSUB0001_pband.ScG0'], applymode='calflagstrict')

Having applied these gain solutions we will once again image the target measurement set which we now expect to have better gain solutions and consequently a better image. We do this by invoking the CLEAN command. 

clean(vis='TSUB0001_pband_target.ms', imagename='0313-0192_clean_sc0', cell=['5.0arcsec','5.0arcsec'], imsize=[4096,4096], mode='mfs',
         nterms=2, gridmode='widefield', wprojplanes=128, stokes='I', niter=10000, spw='3~8,10~15', interactive=True, 
         minpb=0.1, usescrarch=T, weighting='briggs', robust=0.0)

On boxing and cleaning we already notice that the imaging artifacts has reduced significantly. We also see that the target source appears to contain more structure and a greater amount of flux. Further self calibration iterations involving an amplitude & phase gain calibration, even an on target bandpass and leakage calibration are all possible steps that can be introduced into the self-cal imaging loop. These steps and introducing Multi Scale Clean for the extended source are possible improvements on the current strategy implemented.

Appendix: Some P-band data issues you may want to know about

Unfortunately, there are a few known issues with JVLA P-band observing. Most problems were discovered during the P-band commissioning period and have been fixed for newer data sets, but the archive keeps part of this history alive. Here we will go over some of those issues, and (if possible) provide ways of fixing them.

Polarization labeling

For a long period (until recent), P-band feeds have been labelled as being circular (R and L), while they are linear (X and Y). A contributed task =fixlowband()= is available to recognize and fix this problem.

Swapped polarizations

Throughout the whole history of the new VLA low-band system, even up to this day, there have been mistakes in the cabling of the full signal chain. This results in that some antennas have the X-polarization signal come in as Y and vice versa. The way to notice this in visibility data is that for baselines with one swapped antenna most power will be in the cross-hand correlations (XY and YX) rather than the parallel-hand correlations (XX and YY). Once antennas with this feature are identified (manually), a contributed task swappol() is available to fix this problem.

Double data descriptor entries

The data description table is part of the measurement set, and provides a link between the recorded visibilities, the spectral window information and the polarization information. For some data sets, the data description table contains two entries for each P-band spectral window, one pointing towards a circular (RL) polarization definition, and one pointing towards a linear (XY) polarization definition. A contributed task fixlowband() is available to recognize and fix this problem.

Continuous Radio Frequency Interference (RFI)

The wide bandwidth of the low-band receiver is (unfortunately) guaranteed to contain significant amounts of RFI. There are a few RFI sources that are active all the time, and are visible in all array configurations. The default P-band observing setup of 16 x 16 MHz tries to capture as much of the continuous RFI in as few spectral windows as possible, allowing for a simple RFI mitigation strategy in which these spectral windows (spws 1,2 and 9,10, but possibly more) can be immediately flagged.

Two spectral window setup

In the early commissioning period, the default setup for P-band observing was to use 2 spectral windows of 1024 channels each to cover 256 MHz of bandwidth from 230 - 486 MHz. However, it was noticed that strong, narrow-band RFI events were causing data to be lost for the whole spectral window in which they occured. To make the system more robust against data loss from these events, the frequency range was divided into 16 x 16 MHz, and shifted downwards by 6 MHz to capture ever-present RFI into as few spectral windows as possible. If your data has only two P-band spectral windows, please be aware that a higher data loss due to RFI is possible. There is no way to repair this.

Bandpass ripples

Due to signal reflections in cables and within the VLA dish, sinusoidal amplitude (and phase) modulations are always present in P-band data. This is most easily seen in bandpass calibration plots of amplitude versus frequency. In some cases, mostly due to cable & connector problems, these modulations can be very strong (~up to 50 percent of the average amplitude level). These modulations tend not vary over the duration of an observation and can therefore be removed through bandpass (and polarization) calibration. If they are found to be variable (e.g., by inspecting bandpass solutions for separate calibrator scans), the offending antenna / polarization should be flagged.

High cross-polarization

On each antenna, the P-band feed (dipole) is visually aligned with respect to the primary focus support legs. This is normally done within 5-10 degrees accuracy. On rare occasions the dipole on one antenna has accidentally rotated to much larger angles because the locking bolt on the back of the feed was not completely tightened. The result is a high cross-polarization in the baseline visibilities that include this antenna. If this is the case, it is safest to flag this data, since there is a possibility that the dipole has rotated during the observations.

Useful links

Retrieved from "https://casaguides.nrao.edu/index.php?title=User_talk:Preshanthj&oldid=18501"