AT3 for ERIS 2024

Installation guide and data download

Data reduction will be performed in CASA. You won't need any special modules, so please follow the main CASA installation guide. For the VLBI data, please download the following tar file and unpack it in your AT3 work directory (320 MB): n23c2-eris-vlbi.tar.gz

You can also download the same tar file from the following ftp server, e.g. with wget:

​​​​​wget -r ftp.astron.nl:outgoing/jive/marcote/n23c2*

Imaging of the reduced data will be performed in Difmap, which is a software for VLBI data inspection, flagging, imaging and modeling. We created a virtual machine that contains Difmap and also the data we are going to use. Please follow the instructions here: https://github.com/aardk/difmap_tutorial

You need to have Docker and Vagrant installed for the git to work. The installation for these is described in the github link. Do the installation (git clone) in your AT3 work directory.

For the above virtual machine and Difmap to work, Mac users will also need to have 

• Command Line Tools (part of Xcode): https://mac.install.guide/commandlinetools/
• X Window: https://www.xquartz.org/

Once you installed everything and log into the virtual environment (vagrant ssh), you should get a prompt in the command line starting with vagrant@. This is the virtual machine. The data we will use is in ~/data/ of this machine.

You can start up difmap by typing in the vagrant terminal: difmap 
If succsesful, you should see the prompt change to: 0>
You can quit both Difmap and the virtual environment by typing: exit

---

Data reduction in CASA

We will be working on the continuum dataset of N23C2 (see previous talks by Zsolt and Benito for details). It is part of a network monitoring test experiment for the EVN, which is done regularly to make sure that the telescopes are working well and their data can be correlated without issues. We will look at a ~2-hour section which was testing out phase-referencing, which means the telescopes were nodding between two close-by sources in a C-T-C-T-C-... so we could see how well the phase solutions can be transfered between one (C) to the other (T).

Calibrator (C)	J0854+2006
Target (T)	J0905+2052

 

Go to your AT3 work directory and start up casa from there.

cd ~/ERIS/n23c2-eris-vlbi/

casa

 

Loading and inspecting the data

Let's load in the data. Everything in the boxes are typed directly into the CASA terminal (with the CASA<#>: prompt). We will always call in the tasks with default settings for this tutorial and set every important parameter by hand. You can always look at all the parameters by typing inp() and a help on each task with help(<taskname>).

default(importfitsidi)
fitsidifile='n23c2_1_1.IDI'
vis='n23c2.ms'
constobsid=True
scanreindexgap_s=15
importfitsidi()

Load in the external flag table from the station logs. Make it CASA friendly using casavlbitools, which are specialized scripts written for VLBI data reduction in CASA. For this tutorial, don't worry about it and you can download the file from here (copy it into your AT3 directory): n23c2.flag

./casa-vlbi-master/flag.py n23c2.uvflg n23c2_1_1.IDI > n23c2.flag

Read in the flag table.

default(flagdata)
vis='n23c2.ms'
mode='list'
inpfile='n23c2.flag'
flagdata()

Understand the structure of your observation. What is in the file?

default(listobs)
vis='n23c2.ms'
listfile='n23c2.ms.listobs'
overwrite=True
listobs()

Look at the uncalibrated data (amplitude/phase vs frequency and amplitude/phase vs time). We will plot one thing using the parameters and then inspect the rest in the plotms gui. What do we see? Anything out of the ordinary?

default(plotms)
vis='n23c2.ms'
xaxis='frequency'
yaxis='amp'
field='J0854+2006'
avgtime='3600'      
antenna='EF&*'
gridrows=3   
gridcols=3  
iteraxis='antenna'        
coloraxis='scan'
plotms()

 

A priori amplitude calibration

Generate the system temperature information based on the station logs. This information is already a part of the measurement set for this tutorial.

default(gencal)
vis='n23c2.ms'
caltable='n23c2.tsys'
caltype='tsys'
uniform=False
gencal()

Generate the gain curves (amptlitude vs elevation corrections) for each antenna. Already in the ms file.

default(gencal)
vis="n23c2.ms"
caltable= "n23c2.gcal"
caltype="gc"
gencal()

Look at the Tsys corrections, both along frequency and time.

default(plotms)
vis='n23c2.tsys'    
xaxis='frequency'   
yaxis='tsys'        
gridrows=3   
gridcols=3          
coloraxis='corr'    
iteraxis='antenna'
plotms()

xaxis='time'
plotms()

 

Inspect and flag data

A lot of flagging could be done here manually, but we are just going to do a few things to showcase things. Let's plot first the amplitude vs frequency on a single sensitive baseline.

default(plotms)
vis='n23c2.ms'
xaxis='channel'
yaxis='amp'
field='J0854+2006'   
avgtime='3600'
antenna='EF&JB'
iteraxis='spw'
gridrows=1   
gridcols=4          
plotms()

Let's flag the edges of the bands.

default(flagdata)
vis='n23c2.ms'
mode='manual'
spw='*:0~3;60~63'
flagdata()

While we are at it, let's flag also the auto correlations in the data as we don't need that for imaging.

default(flagdata)
vis='n23c2.ms'
mode='manual'
autocorr=True
flagdata()

Let's look at the data again.

default(plotms)
vis='n23c2.ms'
xaxis='channel'
yaxis='amp'
field='J0854+2006'   
avgtime='3600'
antenna='EF&JB'
iteraxis='spw'
gridrows=1   
gridcols=4          
plotms()

 

Fringe fitting

Fringe fitting is where we solve for the residual phase errors across frequency and time on each baseline, so we can align the telescopes to the same wavefront and add them up coherently (see talks from Zsolt and Benito). We will do this in two steps. Let's look at the data first and identify a good timerange for the first delay calculation (I preselected a good timerange already here.)

default(plotms)
vis='n23c2.ms'
xaxis='frequency'
yaxis='phase'
ydatacolumn='data'
antenna='EF&*'
coloraxis='corr'
timerange='14:32:00~14:33:00'
averagedata=True
avgtime='120'
iteraxis='antenna'
gridrows=3   
gridcols=3          
plotms()

Let's use this good chunk of data to calculate the delays (no rates) and model mainly the instrumental errors of the stations.

default(fringefit)
vis='n23c2.ms'
caltable='n23c2.sbd'
timerange='14:32:00~14:33:00'
solint='inf'
zerorates=True
refant='EF'
minsnr=10
gaintable=['n23c2.gcal','n23c2.tsys']
interp=['nearest','nearest,nearest']
parang=True
fringefit()

Apply the sbd solution table to the data (updates the corrected column in the ms file). 

default(applycal)
vis='n23c2.ms'
gaintable=['n23c2.gcal','n23c2.tsys','n23c2.sbd']
interp=['nearest','nearest,nearest','nearest']
parang=True
applycal()

Let's look at the data again before and after the first fringe fitting (removing a lot of phase errors vs frequency).

tget(plotms)
ydatacolumn='corrected'
plotms()

Next: fringe fitting on the phase reference calibrator (in this case it is the same source, but usually this is not the case) to solve for phase residuals across frequency and time.

default(fringefit)
vis='n23c2.ms'
caltable='n23c2.mbd'
field='J0854+2006'
solint='120s'
zerorates=False
refant='EF'
combine='spw'
minsnr=7
gaintable=['n23c2.gcal','n23c2.tsys','n23c2.sbd']
interp=['nearest','nearest,nearest','nearest']
parang=True
fringefit()

Let's look at the mbd solution table (phase and rate, too). Also have a look how is it different to the sbd table.

default(plotms)
vis='n23c2.mbd'
xaxis='time'
yaxis='delay'
gridrows=4
gridcols=4
coloraxis='corr'
iteraxis='antenna'
plotms()

Now apply it to the target (aka phase referencing the target to the calibrator). This is called phase referencing, which is a differential phase calibration method (see the VLBI talks and also the tutorial explanation in person).

default(applycal)
vis='n23c2.ms'
field='J0854+2006,J0905+2052'
gaintable=['n23c2.gcal','n23c2.tsys','n23c2.sbd','n23c2.mbd']
interp=['nearest', 'nearest,nearest','nearest','linear']
spwmap=[[],[],[],[0,0,0,0]]
parang=True
applycal()

Finally, let's look at the calibrated phases across frequency. Are they flat and zero?

default(plotms)
vis='n23c2.ms'
xaxis='frequency'
yaxis='phase'
ydatacolumn='corrected'
antenna='EF&*'
coloraxis='corr'
timerange='14:32:00~14:33:00'
averagedata=True
avgtime='120'
iteraxis='antenna'
gridrows=3   
gridcols=3          
plotms()

And also the phases vs time. Again, are they flat and near zero?

default(plotms)
vis='n23c2.ms'
xaxis='time'
yaxis='phase'
ydatacolumn='corrected'
averagedata=True
avgchannel='64'
avgtime='120'
gridrows=3
gridcols=3
coloraxis='field'
iteraxis='antenna'
plotms()

 

Bandpass calibration

Finally, let's correct for the bandpass shapes. We already did a little before by clipping the edges, but let's do a bit more.

default(bandpass)
vis='n23c2.ms'
caltable='n23c2.bpass'
field='J0854+2006'
gaintable=['n23c2.gcal','n23c2.tsys','n23c2.sbd','n23c2.mbd']
interp=['nearest','nearest,nearest','nearest','linear']
solnorm=True
solint='inf'
refant='EF'
bandtype='B'
spwmap=[[],[],[],[0,0,0,0]]
parang=True
bandpass()

Look at the bandpass solution (it also has a small phase component in this case).

default(plotms)
vis='n23c2.bpass'
xaxis='frequency'
yaxis='amp'
coloraxis='corr'
iteraxis='antenna'
plotms()

Let's apply it to the data.

default(applycal)
vis='n23c2.ms'
field='J0854+2006,J0905+2052'
gaintable=['n23c2.gcal','n23c2.tsys','n23c2.sbd','n23c2.mbd','n23c2.bpass']
interp=['nearest','nearest,nearest','nearest','linear','nearest,nearest']
spwmap=[[],[],[],[0,0,0,0],[]]
parang=True
applycal()

Phases won't change really from before, but let's look at the final corrected amplitudes. First across frequency.

default(plotms)
vis='n23c2.ms'
xaxis='frequency'
yaxis='amp'
ydatacolumn='corrected'
antenna='EF&*'
coloraxis='corr'
timerange='14:32:00~14:33:00'
averagedata=True
avgtime='120'
iteraxis='antenna'
gridrows=3   
gridcols=3          
plotms()

Then across time.

default(plotms)
vis='n23c2.ms'
xaxis='time'
yaxis='amp'
ydatacolumn='corrected'
averagedata=True
avgantenna=True
avgtime='120'
gridrows=3
gridcols=3
coloraxis='field'
iteraxis='antenna'
plotms()

 

Write out calibrated data.

Now that the initial calibration is done, let's apply all the amplitude and phase corrections to the data, so we can average it together in frequency. First make the calibrated measurement set of the calibrator source.

default(mstransform)
vis='n23c2.ms'
outputvis='n23c2_J0854+2006.ms'
field='J0854+2006'
datacolumn='corrected'
keepflags=True
chanaverage=True
chanbin=64
mstransform()

The do the same for the target (making an ms file with a single channel).

default(mstransform)
vis='n23c2.ms'
outputvis='n23c2_J0905+2052.ms'
field='J0905+2052'
datacolumn='corrected'
keepflags=True
chanaverage=True
chanbin=64
mstransform()

Convert the calibrator ms file into UV FITS format, so we could export it from the CASA environment and image it in Difmap.

default(exportuvfits)
vis='n23c2_J0854+2006.ms'
fitsfile='n23c2_J0854+2006.uvfits'
combinespw=True
padwithflags=True
datacolumn='corrected'
exportuvfits()

Finally do the same with the target dataset. We are done with CASA.

default(exportuvfits)
vis='n23c2_J0905+2052.ms'
fitsfile='n23c2_J0905+2052.uvfits'
combinespw=True
padwithflags=True
datacolumn='corrected'
exportuvfits()

---

Imaging in Difmap

Self-calibrate and image in Difmap using the calibrated data in UV FITS format. This is a very interactive thing. I'll write some things here for guidance, but most will be shown directly in the tutorial.

Starting up difmap

Go to your AT3 working directory with the data files and type

difmap

In difmap, set up some macros. You could do this in an external file too, but this is easier now. Again, what you see in the boxes go directly in the 0> prompt.

! radplot settings
rflags="m3"

! vplot settings
vflags="fbm3"

! change model color from red to white
pgscr 2,1,1,1

! snr, calculates & prints the signal-to-noise ratio & peak location in the residual image
#+snr print "signal-to-noise:"; print peak(flux)/imstat(rms); print "peak position:"; print peak(x),peak(y)

! clear model and calibration 
#+clr unshift; delwin; clrmod true; uncal true,true

Main calibration steps

Load in the file.

observe n23c2_J0854+2006.uvfits

Inspect that data.

header, tplot, uvplot, radplot, projplot, vplot

Clean and selfcal (phase-only)

mapsize, maplot, clean, selfcal, radpl, vpl, uvw

Amplitude selfcal

gscale, mapl, clean, selfcal true,true

Final map (making a nice map)

cmul, imstat, mapl cln

Model fitting (getting out some parameters from the data)

clrmod, modelfit, save