# 2008 NRAO Synthesis Summer School

# VLA Continuum Polarimetry Tutorial - CASA 

#

# gmoellen (2008Jun08)

#

# This is a script that reduces B- and C-configuration VLA continuum

#  polarimetry data at 5 GHz on 3C129 and calibrators, and images the

#  dual-config 3C129 in full polarization.

#

# Observations:  AT166

#                B-array: 1994Jul25

#                C-array: 1994Nov03

#

# Sources:   Science Target: 3C129, a radio galaxy

#            Calibrator:     0420+417, observed altenately with 3C129

#            Calibrator:     0518+165 (3C138), for flux density/poln p.a.

#            Calibrator:     0134+329 (3C48), for flux density/poln p.a.

#

# Obs Modes: Two 50 MHz continuum (single chan) sub-bands at 4585.1 & 4885.1 MHz

#            Full polarization: RR,RL,LR,LL

#



# In the following script, the B array and C array data are separately

#  filled and reduced.  The calibration parameters are fairly standard,

#  and some variations may be suggested.  In general, after each calibration

#  opearation, a data or calibration plotting command is included, again

#  with suggestions on possible variations.

#

# This script is not written to be run automatically; doing so will likely

#  yield a less-than-optimal result.  It is intended that students will

#  cut-and-paste this script one task at a time to progress through the

#  reduction.

# The methods are run in the function-call style:  task(param1=x,param2=y)

#  To see the range of parameters for a function, use 'inp <task>'.  After

#  running a task (in function-call style), type 'tget <task>', and

#  'inp <task>' to review the parameters that were used.





# A few parameters used repeatedly below, are defined here:

#  (Be sure to reset these if you exit CASA and start it up again)

msnameB='at166B.ms';

msnameC='at166C.ms';





# Fill B-config data at C-band (5 GHz)

importvla(archivefiles=['AT166_1', 'AT166_2'],vis=msnameB,bandname='C');



# List a summary of the dataset in the logger

listobs(vis=msnameB);



#--> Note scan sequence, fields, and spectral window information, etc.



# set flux density calibrator total intensity models

#  NB: these sources are resolved, so we use model images provided

#      by the VLA (copy them from the data repository: data/VLA/CalModels/)

#  NB: By default, the model for 0420+417 is a 1 Jy point source

setjy(vis=msnameB,field='0518+165',modimage='3C138_C.im');

setjy(vis=msnameB,field='0134+329',modimage='3C48_C.im');



# Plot data and interactively edit

#  One spw at a time, calibrators only, RR, LL only

plotxy(vis=msnameB,spw='0',field='0420+417,0518+165,0134+329',selectdata=T,correlation='RR,LL');

plotxy(vis=msnameB,spw='1',field='0420+417,0518+165,0134+329',selectdata=T,correlation='RR,LL');



#--> Variations: plot on different axes (e.g., phase vs. time, real vs imag, amp/phase vs. uvdist)

#-->             plot model data for calibrators

#-->             plot 3C129 data



# Solve for gains on calibrators

#  NB: reference phases to VA12; pre-apply parallactic angle correction

gaincal(vis=msnameB,caltable='at166B.gcal',field='0420+417,0518+165,0134+329',refant='VA12',parang=T);



#--> (22/22 good solutions)



# Examine solutions:

plotcal(caltable='at166B.gcal',yaxis='amp');

plotcal(caltable='at166B.gcal',yaxis='phase');

#--> Variations:  plot solutions per antenna



# Scale gain solution from 0420+417 according to f.d. calibrators

fluxscale(vis=msnameB,caltable='at166B.gcal',fluxtable='at166B.fcal',reference='0518+165,0134+329');



# -->  0420: 1.441/1.443 Jy (broadly consistent with VLA Cal man values at L & X)



# Examine solutions:

plotcal(caltable='at166B.fcal',yaxis='amp');



#--> Note that gain amps are ~constant now 



# Solve for instrumental polarization on 0420+417 (and also for source poln)

# NB: IMPORTANT: use gcal---not fcal---here because model is _still_ 1.0 Jy

polcal(vis=msnameB,caltable='at166B.dcal',field='0420+417',refant='VA12',gaintable='at166B.gcal',gainfield='0420+417');



# --> 2/2 good solutions;

# --> 0420: 0.027Jy @ 61.2deg / 0.024Jy @ -83.0 Jy



# Examine solutions:

plotcal(caltable='at166B.dcal',xaxis='antenna',yaxis='amp');

plotcal(caltable='at166B.dcal',xaxis='real',yaxis='imag',plotrange=[-0.05,0.05,-0.05,0.05]);



# Set full polarization model for 0518+165 (pol is 11.1% @ -11 deg  [RL = -22])

#  NB: neglecting source structure here)

setjy(vis=msnameB,field='0518+165',spw='0',fluxdensity=[3.688, 0.380, -0.153, 0.0]);

setjy(vis=msnameB,field='0518+165',spw='1',fluxdensity=[3.862, 0.397, -0.161, 0.0]);



# Solve for polarization position angle on 0518+165

polcal(vis='at166B.ms',caltable='at166B.xcal',field='0518+165',refant='VA12',poltype='X',gaintable=['at166B.fcal', 'at166B.dcal'],gainfield=['0518+165','0420+417']);



# --> 2/2 good solutions

# --> 0518:  2.7deg  / 41.8deg  (I need to check these numbers)





# apply all calibration...

#  (NB: different fields are selected from each caltable, depending on selected data fields)

#  (NB: parang=T is set to rotate polarization p.a. frame from antennas to sky



#...to 0420 & 3C129 

applycal(vis=msnameB,field='0420+417,3C129',gaintable=['at166B.fcal','at166B.dcal','at166B.xcal'],gainfield=['0420+417','0420+417','0518+165'],parang=T);

#...to 0518 

applycal(vis=msnameB,field='0518+165',gaintable=['at166B.fcal','at166B.dcal','at166B.xcal'],gainfield=['0518+165','0420+417','0518+165'],parang=T);

#...to 0134 

applycal(vis=msnameB,field='0134+329',gaintable=['at166B.fcal','at166B.dcal','at166B.xcal'],gainfield=['0134+329','0420+417','0518+165'],parang=T);





# Examine (edit?) calibrated data (calibrators)

plotxy(vis=msnameB,datacolumn='corrected',field='0420+417,0518+165,0134+329',selectdata=T,correlation='RR,LL');





#--> Variations: plot on different axes (e.g., phase vs. time, real vs imag)

#-->             overlay model data for calibrators

#-->             plot 3C129 data



# Examine cross-hand data (real vs. imag)

plotxy(vis=msnameB,xaxis='real',yaxis='imag',datacolumn='corrected',selectdata=True,correlation='RL,LR',plotrange=[-0.5,0.5,-0.5,0.5],field = '0518+165');

plotxy(vis=msnameB,xaxis='real',yaxis='imag',datacolumn='corrected',selectdata=True,correlation='RL,LR',plotrange=[-0.5,0.5,-0.5,0.5],field = '0134+329');

plotxy(vis=msnameB,xaxis='real',yaxis='imag',datacolumn='corrected',selectdata=True,correlation='RL,LR',plotrange=[-0.5,0.5,-0.5,0.5],field = '0420+417');



#--> NB: RL and LR signal are complex conjugates of each other (Q+iU & Q-iU)





# do some simple non-interactive imaging of each source

#  NB: using hogbom cleanso pol planes cleaned separately



# 0518:

clean(vis=msnameB,imagename='0518B',field='0518+165',psfmode='hogbom',niter=500,imsize=[512,512],cell=['0.3arcsec','0.3arcsec'],stokes='IQUV',weighting='briggs',robust=0.5,interactive=F);



# 0134:

clean(vis=msnameB,imagename='0134B',field='0134+329',psfmode='hogbom',niter=500,imsize=[512,512],cell=['0.3arcsec','0.3arcsec'],stokes='IQUV',weighting='briggs',robust=0.5,interactive=F);



# 0420:

clean(vis=msnameB,imagename='0420B',field='0420+417',psfmode='hogbom',niter=500,imsize=[512,512],cell=['0.3arcsec','0.3arcsec'],stokes='IQUV',weighting='briggs',robust=0.5,interactive=F);



# 3C129:

clean(vis=msnameB,imagename='3c129B',field='3C129',psfmode='hogbom',niter=2500,imsize=[2048,2048],cell=['0.3arcsec','0.3arcsec'],stokes='IQUV',weighting='briggs',robust=0.5,interactive=F);



# Examine images in viewer...



#  Questions:  Do calibrator polarizations come out right?

#              Is calibrator structure as expected?

#              Is 3C129 image any good?  (cf C-config imaging below)



# Variations:  clean interactively (clean boxes)

#              try different weighting 





# - - - - - - - - - - - - - - -



# NOW, reduce C-config data in much the same way...



# Fill C-config data at C-band (5 GHz)

importvla(archivefiles='AT166_3',vis=msnameC,bandname='C');



# List a summary of the dataset in the logger

listobs(vis=msnameC);



# set flux density calibrator total intensity models

#  NB: these sources are resolved, so we use model images provided

#      by the VLA

#  NB: By default, the model for 0420+417 is a 1 Jy point source

setjy(vis=msnameC,field='0518+165',modimage='3C138_C.im');

setjy(vis=msnameC,field='0134+329',modimage='3C48_C.im');



# Plot data and interactively edit

#  One spw at a time, calibrators only, RR, LL only

plotxy(vis=msnameC,spw='0',field='0420+417,0518+165,0134+329',selectdata=T,correlation='RR,LL');

plotxy(vis=msnameC,spw='1',field='0420+417,0518+165,0134+329',selectdata=T,correlation='RR,LL');



#--> Variations: plot on different axes (e.g., phase vs. time, real vs imag, amp/phase vs. uvdist)

#-->             plot model data for calibrators (datacolumn='model')

#-->             plot 3C129 data

#-->             plot individual baselines (iteration='baseline')



# Solve for gains on calibrators

#  NB: reference phases to VA12; pre-apply parallactic angle correction

gaincal(vis=msnameC,caltable='at166C.gcal',field='0420+417,0518+165,0134+329',refant='VA12',parang=T,solint=500);



#--> (12/14 good solutions)  NB: no solutions for 0134+329 (too little data?)



# Examine solutions:

plotcal(caltable='at166C.gcal',yaxis='amp');

plotcal(caltable='at166C.gcal',yaxis='phase');

#--> Variations:  plot solutions per antenna, spw, etc.



# Scale gain solution from 0420+417 according to f.d. calibrators

fluxscale(vis=msnameC,caltable='at166C.gcal',fluxtable='at166C.fcal',reference='0518+165');

# -->  0420:  1.282/1.292 Jy (broadly consistent with VLA Cal man values at L & X)



# Examine solutions:

plotcal(caltable='at166C.fcal',yaxis='amp');



#--> gain amplitudes now ~constant



# Solve for instrumental polarization on 0420+417 (and also for source poln)

# NB: IMPORTANT: use gcal---not fcal---here because model is _still_ 1.0 Jy

polcal(vis=msnameC,caltable='at166C.dcal',field='0420+417',refant='VA12',gaintable='at166C.gcal',gainfield='0420+417');

# --> 2/2 good solutions;

# --> 0420: 0.035Jy @ -41.4deg / 0.033Jy @ 18.1deg



# Examine solutions:

plotcal(caltable='at166C.dcal',xaxis='antenna',yaxis='amp');

plotcal(caltable='at166C.dcal',xaxis='real',yaxis='imag',plotrange=[-0.05,0.05,-0.05,0.05]);



# Set full polarization model for 0518+165 (pol is 11.1% @ -11 deg  [RL = -22])

#  NB: neglecting source structure here)

setjy(vis=msnameC,field='0518+165',spw='0',fluxdensity=[3.688, 0.380, -0.153, 0.0]);

setjy(vis=msnameC,field='0518+165',spw='1',fluxdensity=[3.862, 0.397, -0.161, 0.0]);



# Solve for polarization position angle on 0518+165

polcal(vis='at166C.ms',caltable='at166C.xcal',field='0518+165',refant='VA12',poltype='X',gaintable=['at166C.fcal', 'at166C.dcal'],gainfield=['0518+165','0420+417']);

# --> 2/2 good solutions

# --> 0518:  77.0deg  / -42.0deg 





# apply all calibration...

#  (NB: different fields are selected from each caltable, depending on selected data fields)

#  (NB: parang=T is set to rotate polarization p.a. frame from antennas to sky



#...to 0420 & 3C129

applycal(vis=msnameC,field='0420+417,3C129',gaintable=['at166C.fcal','at166C.dcal','at166C.xcal'],gainfield=['0420+417','0420+417','0518+165'],parang=T);

#...to 0518

applycal(vis=msnameC,field='0518+165',gaintable=['at166C.fcal','at166C.dcal','at166C.xcal'],gainfield=['0518+165','0420+417','0518+165'],parang=T);





# Examine (edit?) calibrated data (calibrators)

plotxy(vis=msnameC,datacolumn='corrected',field='0420+417,0518+165',selectdata=T,correlation='RR,LL');





#--> Variations: plot on different axes (e.g., phase vs. time, real vs imag)

#-->             overlay model data for calibrators

#-->             plot 3C129 data



# Examine cross-hand data (real vs. imag)

plotxy(vis=msnameC,xaxis='real',yaxis='imag',datacolumn='corrected',selectdata=True,correlation='RL,LR',plotrange=[-0.5,0.5,-0.5,0.5],field = '0518+165');

plotxy(vis=msnameC,xaxis='real',yaxis='imag',datacolumn='corrected',selectdata=True,correlation='RL,LR',plotrange=[-0.5,0.5,-0.5,0.5],field = '0420+417');



#--> NB: RL and LR signal are complex conjugates of each other (Q+iU & Q-iU)



# do some simple non-interactive imaging of each source

#  NB: using hogbom cleanso pol planes cleaned separately



# 0518:

clean(vis=msnameC,imagename='0518C',field='0518+165',psfmode='hogbom',niter=500,imsize=[512,512],cell=['0.3arcsec','0.3arcsec'],stokes='IQUV',weighting='briggs',robust=0.5,interactive=F);



# 0420:

clean(vis=msnameC,imagename='0420C',field='0420+417',psfmode='hogbom',niter=500,imsize=[512,512],cell=['0.3arcsec','0.3arcsec'],stokes='IQUV',weighting='briggs',robust=0.5,interactive=F);



# 3C129:

clean(vis=msnameC,imagename='3c129C',field='3C129',psfmode='hogbom',niter=2500,imsize=[2048,2048],cell=['0.3arcsec','0.3arcsec'],stokes='IQUV',weighting='briggs',robust=0.5,interactive=F);



# Examine images in viewer...



#  Questions:  Do calibrator polarizations come out right?

#              Is calibrator structure as expected?

#              Is 3C129 image good?  (cf B-config imaging above)



# Variations:  clean interactively (clean boxes)

#              try different weighting 

#              in viewer, blink between B and C 3C129 images (in all polarizations)





# - - - - - - - - - - - - - - -



# The next steps extract the 3C129 data from the above datasets, and combine

#  them to permit a dual-config imaging run



# split out 3C129

split(vis=msnameB,outputvis='3C129B.ms',field='3C129');

split(vis=msnameC,outputvis='3C129C.ms',field='3C129');



# make one MS so we can image the combined config

concat(vis=['3C129B.ms','3C129C.ms'],concatvis='3C129BC.ms');

listobs(vis='3C129BC.ms');



# simple clean

clean(vis='3C129BC.ms',imagename='3c129BC',field='3C129',psfmode='hogbom',niter=5000,imsize=[2048,2048],cell=['0.3arcsec','0.3arcsec'],stokes='IQUV',weighting='briggs',robust=0.5,interactive=F);



# Examine image in viewer...



#  Questions:  How has the 3C129 image changed?  Is it better?



# Variations:  clean interactively (clean boxes)

#              try different weighting 

#              split/concat/image dual-config calibrator data

#              self-cal?



# Extract Q an U images

immath(stokes='Q', outfile='3c129BC.Q', mode='evalexpr', expr='\'3c129BC.image\'');

immath(stokes='U', outfile='3c129BC.U', mode='evalexpr', expr='\'3c129BC.image\'');



# Form poln intensity and pos ang

immath(stokes='', outfile='3c129BC.P', mode='poli', imagename=['3c129BC.Q','3c129BC.U'], sigma='0.0mJy/beam');

immath(stokes='', outfile='3c129BC.X', mode='pola', imagename=['3c129BC.Q','3c129BC.U'], sigma='0.0mJy/beam');



# In viewer 'Load Data' window, use the following LEL expression to load the p.a. image:

#   '3c129BC.X'['3c129BC.P'>0.0001]