# this guide requires the current casapy-stable release (or later):
# At http://casa.nrao.edu/casa_obtaining.shtml the "stable" files contain the revision number 15777

############################################################
########## Initial inspection & a Priori Flagging
############################################################

data = ['uid___A002_X1d7c9e_X314_wvrtsys.ms',
'uid___A002_X1db0ee_X2a8_wvrtsys.ms','uid___A002_X1db0ee_X44d_wvrtsys.ms']

for vis in data:
  listobs(vis=vis,listfile='listobs_'+vis+'.txt') 

for vis in data:
  flagdata(vis=vis,autocorr=True, flagbackup=T)
  flagdata(vis=vis,mode='shadow',diameter=12.0, flagbackup=F)

################################
# Plots ########################

vis=data[0]
plotms(vis=vis,xaxis='time',yaxis='amp',field='',avgchannel='3840',
       antenna='',coloraxis='field',iteraxis='spw',xselfscale=T,ydatacolumn='data')

# spw=1 PM03 YY
# spw=0,1 06:21:05.4; DV06
# spw=0,1   06:02:23.4; DV04
# spw=0,1   06:18:10.4; DV09

vis=data[1]
plotms(vis=vis,xaxis='time',yaxis='amp',field='',avgchannel='3840',
       antenna='',coloraxis='field',iteraxis='spw',xselfscale=T,ydatacolumn='data')

# spw=1; PM03
# spw=0,1; DV09
# spw=0,1; DV04 dropouts: scans 12,14,16,18

vis=data[2]
plotms(vis=vis,xaxis='time',yaxis='amp',field='',avgchannel='3840',
       antenna='',coloraxis='field',iteraxis='spw',xselfscale=T,ydatacolumn='data')

# spw=1; PM03
# spw=0,1; DV09
# spw=0,1; DV04 dropouts: scans 8,14,22

vis = data[0]
plotms(vis=vis,xaxis='frequency',yaxis='amp',field='3C279,Titan,J1037-295=QSO',
       antenna='',avgtime='1e8',avgscan=T,coloraxis='field',iteraxis='spw',
       xselfscale=T,ydatacolumn='data')
# Titan has line in spw0 (no line free channels)
# 3C279 can see the mesospherica absorption in spw=0

vis = data[1]
plotms(vis=vis,xaxis='frequency',yaxis='amp',field='3C279,Titan,J1037-295=QSO',
       antenna='',avgtime='1e8',avgscan=T,coloraxis='field',iteraxis='spw',
       xselfscale=T,ydatacolumn='data')


vis = data[2]
plotms(vis=vis,xaxis='frequency',yaxis='amp',field='3C279,Titan,J1037-295=QSO',
       antenna='',avgtime='1e8',avgscan=T,coloraxis='field',iteraxis='spw',
       xselfscale=T,ydatacolumn='data')

########################
# Initial Flagging #
########################

data = ['uid___A002_X1d7c9e_X314_wvrtsys.ms',
'uid___A002_X1db0ee_X2a8_wvrtsys.ms','uid___A002_X1db0ee_X44d_wvrtsys.ms']

# Specific flagging for all ms files
for vis in data:
  flagdata(vis=vis,spw='1',antenna='PM03',correlation='YY', flagbackup=T) 

# Specific flagging for some ms files

flagdata(vis=data[0],antenna='DV06', timerange='06:21:05.0~06:21:05.9',flagbackup=T) 
flagdata(vis=data[0],antenna='DV04', timerange='06:02:23.0~06:02:23.9',flagbackup=T) 
flagdata(vis=data[0],antenna='DV09', timerange='06:18:10.0~06:18:10.9',flagbackup=T) 

flagdata(vis=data[1],antenna='DV09', flagbackup=T) # low and erratic amplitudes
flagdata(vis=data[1],antenna='DV04', scan='12,14,16,18',flagbackup=T) 

flagdata(vis=data[2],antenna='DV09', flagbackup=T) # low and erratic amplitudes
flagdata(vis=data[2],antenna='DV04', scan='8,14,22',flagbackup=T) 


# Look at uv-plots of Titan to be sure there aren't any problems with
# it being close to Saturn.

vis=data[0]
plotms(vis=vis,xaxis='uvdist',yaxis='amp',field='Titan',avgchannel='3840',
       antenna='',coloraxis='spw',ydatacolumn='data')

vis=data[1]
plotms(vis=vis,xaxis='uvdist',yaxis='amp',field='Titan',avgchannel='3840',
       antenna='',coloraxis='spw',ydatacolumn='data')

vis=data[2]
plotms(vis=vis,xaxis='uvdist',yaxis='amp',field='Titan',avgchannel='3840',
       antenna='',coloraxis='spw',ydatacolumn='data')

######################################
# Concanenate the data sets together #
######################################

# you may get warning messages saying that "Found invalid antenna ids
# in the POINTING table" -> this is ok
concat(vis=data,concatvis='Band6multi_april.ms')

data='Band6multi_april.ms'
listobs(vis=data, listfile='listobs_concat.txt')
# Note that the SrcId is not unique after concatenating files together
# That means its essential to use field names instead of ids in
# subsequent calibration

###########################
# Fix Titan's Coordinates #
###########################

# Bring each of the Titan observations to a common position

# you may get warning messages saying that "PHASE_DIR is not a variable
# reference column." -> this is ok.

data='Band6multi_april.ms'
fixplanets(vis=data, field='Titan', fixuvw=True, 
       direction='J2000 12h50m00.65515 -02d25m29.0209')  


listobs(vis=data, listfile='listobs_fixplanets.txt')

######################################
# Check calibrator spectral features #
######################################

data='Band6multi_april.ms'

plotms(vis=data,spw='',xaxis='channel',yaxis='amp',field='3C279',
       avgtime='1e8',iteraxis='spw',ydatacolumn='data')
plotms(vis=data,spw='',xaxis='channel',yaxis='amp',field='J1037-295=QSO',
       avgtime='1e8',iteraxis='spw',ydatacolumn='data')
plotms(vis=data,spw='',xaxis='channel',yaxis='amp',field='Titan',
       avgtime='1e8',iteraxis='spw',ydatacolumn='data')
# Nothing to see here, other than CO 2-1 in spw0

flagdata(vis=data,spw='0:1871~1885',flagbackup=T) # Flag atmospheric CO 2-1 line

######################################
# Post calibration flagging #
######################################

# The amplitude vs time plots taught us to flag the following data.

flagdata(vis=data,spw='',field='3C279',antenna='DV04',
         timerange='73:28:05~73:28:10',flagbackup=T)

flagdata(vis=data,spw='',field='Titan',antenna='',
         timerange='06:12:21~06:12:23',flagbackup=T)

################################
# Set up flux calibrator model #
################################

clearcal('Band6multi_april.ms')

data='Band6multi_april.ms'

# have to do them separately to get around issue with setjy and fact
# that Titan now has one name but 3 field ids.

setjy(vis=data,field='3',standard='Butler-JPL-Horizons 2010',scalebychan=T)
setjy(vis=data,field='4',standard='Butler-JPL-Horizons 2010',scalebychan=T)
setjy(vis=data,field='5',standard='Butler-JPL-Horizons 2010',scalebychan=T)

plotms(vis=data,field='Titan',xaxis='uvdist',yaxis='amp',ydatacolumn='model',
       avgchannel='3840',coloraxis='spw')
# Looks pretty constant between the data sets

########################
# Bandpass calibration #
########################
plotms(vis=data,spw='',xaxis='time',yaxis='amp',field='3C279',avgchannel='3840',
       iteraxis='spw',ydatacolumn='data')

plotms(vis=data,spw='0',xaxis='time',yaxis='phase',
       antenna='2',coloraxis='corr',field='3C279',avgchannel='3840',iteraxis='baseline')

gaincal(vis=data,caltable='bpphase.gcal',field='3C279',spw='*:900~1100',refant='DV07',
        calmode='p',solint='int',minsnr=2.0,minblperant=4)

bandpass(vis=data,caltable='bandpass.bcal',field='3C279',spw='',combine='',refant='DV07',
        solint='inf',solnorm=T,minblperant=4,fillgaps=17,gaintable=['bpphase.gcal'])


plotcal(caltable='bpphase.gcal',xaxis='time',yaxis='phase',spw='',antenna='', 
        iteration='antenna',subplot=421,plotrange=[0,0,-180,180],poln='')

plotcal(caltable='bandpass.bcal',xaxis='freq',yaxis='phase',antenna='', 
        iteration='antenna',subplot=421,spw='0')

####################
# Gain calibration #
####################

data='Band6multi_april.ms'

# phase only solution on integration timescale
gaincal(vis=data,caltable='intphase.gcal',
        field='Titan,3C279,J1037-295=QSO',spw='*:40~3800',refant='DV07',
        calmode='p',solint='int',minsnr=2.0,minblperant=4,
        gaintable=['bandpass.bcal'])

# phase only solution on scan timescale
gaincal(vis=data,caltable='scanphase.gcal',
        field='Titan,3C279,J1037-295=QSO',spw='*:40~3800',refant='DV07',
        calmode='p',solint='inf',minsnr=2.0,minblperant=4,
        gaintable=['bandpass.bcal'])

# amp & phase solution on scan timescale
gaincal(vis=data,caltable='amp.gcal',
        field='Titan,3C279,J1037-295=QSO',spw='*:40~3800',refant='DV07',
        calmode='ap',solint='inf',minsnr=2.0,minblperant=4,
        gaintable=['bandpass.bcal','intphase.gcal'])

plotcal(caltable='intphase.gcal',xaxis='time',yaxis='phase',antenna='', 
        spw='',field='',iteration='antenna',subplot=521,
        plotrange=[0,0,-180,180],poln='X',figfile='intphase_X.png')

# Use this plot to see what to flag
plotcal(caltable='intphase.gcal',xaxis='time',yaxis='phase', 
        spw='',antenna='DV04',iteration='field',subplot=111,
        plotrange=[0,0,-180,180])

plotcal(caltable='scanphase.gcal',xaxis='time',yaxis='phase',antenna='', 
        spw='',field='',iteration='antenna',subplot=521,
        plotrange=[0,0,-180,180],poln='X',figfile='scanphase_X.png')

plotcal(caltable='scanphase.gcal',xaxis='time',yaxis='phase',antenna='', 
        spw='',field='',iteration='antenna',subplot=521,
        plotrange=[0,0,-180,180],poln='Y',figfile='scanphase_Y.png')

plotcal(caltable='amp.gcal',xaxis='time',yaxis='phase',antenna='', 
        spw='',field='',plotrange=[0,0,-1,1],
        iteration='antenna',subplot=521,figfile='amp_phase.png')

plotcal(caltable='amp.gcal',xaxis='time',yaxis='amp',antenna='', 
        iteration='antenna',subplot=521,spw='',poln='X',
        plotrange=[0,0,0.0,1.0],figfile='amp_X.png')
plotcal(caltable='amp.gcal',xaxis='time',yaxis='amp',antenna='', 
        iteration='field',subplot=521,spw='',poln='Y',
        plotrange=[0,0,0.0,1.0])


######################
# Set the flux scale #
######################

# don't use spw 0 for flux calibration because of CO 2-1 emisssion
fluxscale(vis=data,caltable='amp.gcal',
          fluxtable='flux.cal',reference='Titan',refspwmap=[1,1])

plotcal(caltable='flux.cal',xaxis='time',yaxis='amp',antenna='', 
        iteration='antenna',subplot=521,spw='',poln='X',
        plotrange=[0,0,0.0,0.3],figfile='flux_X.png')
plotcal(caltable='flux.cal',xaxis='time',yaxis='amp',antenna='', 
        iteration='antenna',subplot=521,spw='',poln='Y',
        plotrange=[0,0,0.0,0.3],figfile='flux_Y.png')
 

#Spw=0 will be referenced to spw=1
#Flux density for J1037-295=QSO in SpW=0 (ref SpW=1) is: 1.28088 +/- 0.00361902 (SNR = 353.93, nAnt= 7)
#Flux density for J1037-295=QSO in SpW=1 is: 1.36127 +/- 0.0036501 (SNR = 372.942, nAnt= 7)
#Flux density for 3C279 in SpW=0 (ref SpW=1) is: 13.4248 +/- 0.060322 (SNR = 222.552, nAnt= 8)
#Flux density for 3C279 in SpW=1 is: 14.2002 +/- 0.0682971 (SNR = 207.918, nAnt= 8)

#################################
# Apply calibration and inspect #
#################################

data='Band6multi_april.ms'

flagmanager(vis=data,mode='save',versionname='before_applycal')

applycal(vis=data,field='3C279',
        gaintable=['bandpass.bcal','intphase.gcal','flux.cal'],
        interp=['nearest','nearest','nearest'],
        gainfield=['3C279','3C279','3C279'],flagbackup=T)

applycal(vis=data,field='J1037-295=QSO',
        gaintable=['bandpass.bcal','intphase.gcal','flux.cal'],
        interp=['nearest','nearest','nearest'],
        gainfield=['3C279','J1037-295=QSO','J1037-295=QSO'],flagbackup=T)

applycal(vis=data,field='Titan',
        gaintable=['bandpass.bcal','intphase.gcal','flux.cal'],
        interp=['nearest','nearest','nearest'],
        gainfield=['3C279','Titan','Titan'],flagbackup=T)

applycal(vis=data,field='TW Hya',
        gaintable=['bandpass.bcal','scanphase.gcal','flux.cal'],
        interp=['nearest','nearest','nearest'],
        gainfield=['3C279','J1037-295=QSO','J1037-295=QSO'],flagbackup=T)

### Inspection section!

## Look at each field, amplitude vs time
plotms(vis=data,spw='',xaxis='time',yaxis='amp',avgchannel='3840',
       coloraxis='field',iteraxis='spw',ydatacolumn='corrected')

plotms(vis=data,spw='',xaxis='time',yaxis='phase',field='3C279,J1037-295=QSO',
       avgchannel='3840',coloraxis='field',iteraxis='spw',ydatacolumn='corrected',
       xselfscale=T)

## Look at each spectral window, frequency vs time
plotms(vis=data,spw='',xaxis='frequency',yaxis='amp',field='',avgtime='1e8',
       coloraxis='baseline',iteraxis='field',ydatacolumn='corrected')

plotms(vis=data,spw='',xaxis='uvdist',yaxis='amp',field='Titan',avgchannel='3840',
       coloraxis='baseline',iteraxis='spw',ydatacolumn='corrected')

## Inspect: does velocity match what's expected for TW Hya?  2.88km/s in LSRK
plotms(vis=data,spw='0',xaxis='velocity',yaxis='amp',field='TW Hya',
       avgtime='1e8',avgscan=T,transform=T,freqframe='LSRK',restfreq='230.5380GHz',
       plotrange=[-10,15,0,0],coloraxis='spw')


# Split off the data
# you may get warning messages saying that "Invalid SOURCE ID
# in SOURCE table row x" -> this is ok
split(vis=data, outputvis='TWHydra_corrected.ms', datacolumn='corrected', field='TW Hya')
split(vis=data, outputvis='Titan_corrected.ms', datacolumn='corrected', field='Titan')
split(vis=data, outputvis='3C279_corrected.ms', datacolumn='corrected', field='3C279')
split(vis=data, outputvis='J1037-295_corrected.ms', datacolumn='corrected', field='J1037-295=QSO')

# Imaging of calibrators is left as an exercise for the reader.