# Antennae band 6 calibration script for CASA 3.3

#South (HA coverage ~ -0.2 - 3.0)
#0. uid___A002_X1fe255_X1f6.ms    RA ~ 1.8 - 3.0
#1. uid___A002_X21487a_X7e.ms     RA ~ -0.2 - 0.9 
#3. uid___A002_X219601_X13c.ms    RA ~ 1.0 - 2.1

#North (HA coverage ~ -1.4 - 4.0)
#4. uid___A002_X1fe255_X51.ms     RA ~ 0.5 - 1.5
#5. uid___A002_X1fe255_X39b.ms    RA ~ 3.2 - 4.0
#6. uid___A002_X20b4e9_X19.ms     RA ~ -0.6 - -0.1
#7. uid___A002_X20bcff_X1c.ms     RA ~ 2.4 - 3.4
#8. uid___A002_X21bf10_X1d8.ms    RA ~ 2.0 - 3.0
#9. uid___A002_X22afc4_X42.ms     RA ~ -1.4 - -0.3
#10. uid___A002_X22afc4_X24f.ms   RA ~ 1.3 - 2.3
#11. uid___A002_X22afc4_X3f4.ms   RA ~ 2.6 - 3.0                 

# Set this option to true if you want to make calibration plots
# - note that this may take a long time
# - default is no plots (you can make them later anyway)
calplots=F

# Set this option to true if you want to run this script interactively to make inspection plots
# - you'll need to hit <enter> at various stages to continue the script
# - default is true (user input is required)
# - set to false if you want no user interaction except clean
interact=T 

basename_all=['uid___A002_X1fe255_X1f6','uid___A002_X21487a_X7e','uid___A002_X219601_X13c','uid___A002_X1fe255_X51','uid___A002_X1fe255_X39b','uid___A002_X20b4e9_X19','uid___A002_X20bcff_X1c','uid___A002_X21bf10_X1d8','uid___A002_X22afc4_X42','uid___A002_X22afc4_X24f','uid___A002_X22afc4_X3f4']

basename_all_south=['uid___A002_X1fe255_X1f6','uid___A002_X21487a_X7e','uid___A002_X219601_X13c']
basename_all_north=['uid___A002_X1fe255_X51','uid___A002_X1fe255_X39b','uid___A002_X20b4e9_X19','uid___A002_X20bcff_X1c','uid___A002_X21bf10_X1d8','uid___A002_X22afc4_X42','uid___A002_X22afc4_X24f','uid___A002_X22afc4_X3f4']

#Create summary file for each dataset
for asdm in basename_all:
  os.system('rm '+asdm+'.listobs.txt')
  listobs(vis=asdm+'.ms', listfile=asdm+'.listobs.txt', verbose=True)

#Print examples of summary for each mosaic
print 'Example south,  uid___A002_X1fe255_X1f6.listobs.txt'
cat uid___A002_X1fe255_X1f6.listobs.txt

print 'Example north,  uid___A002_X1fe255_X51.listobs.txt'
cat uid___A002_X1fe255_X51.listobs.txt

for asdm in basename_all:
  print "Antenna configuration for : "+asdm
  plotants(vis=asdm+'.ms', figfile=asdm+'.plotants.png')
  if(interact):
    dummy_string = raw_input("Hit <Enter> to see the antenna configuration for the next data set.")

# Fix planet/satellite position using fixplanets
for asdm in basename_all:
  fixplanets(vis=asdm+'.ms', field='Titan', fixuvw=True)

# Examining amplitude vs. time for a single dataset 
asdm=basename_all[0]
plotms(vis=asdm+'.ms',
  spw='1',
  xaxis='time', yaxis='amp',
  correlation='XX',
  antenna='*&*',
  avgchannel='10000',coloraxis='field')
dummy_string = raw_input("Hit <Enter> to see the antenna configuration for the next data set.")

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

for asdm in basename_all:
  # To remove all the previous flags
  flagdata(vis=asdm+'.ms',mode='manualflag', unflag=T, flagbackup=F)
  # Flagging shadowed data
  flagdata(vis=asdm+'.ms',mode = 'shadow', diameter=12.0, flagbackup = F)
  # Flagging pointing and atmosphere calibration scans
  flagdata(vis=asdm+'.ms', mode='manualflag', intent='*POINTING*,*ATMOSPHERE*',flagbackup = F)
  # Flagging autocorrelation data
  flagdata(vis=asdm+'.ms',autocorr=True,flagbackup=F)
  # Saving 'a priori' flags
  flagmanager(vis = asdm+'.ms', mode = 'save', versionname = 'Apriori')

# Restoring up 'a priori' flags"+asdm
#for asdm in basename_all:
#  flagmanager(vis = asdm+'.ms', mode = 'restore', versionname = 'Apriori')

######## Examine and Apply Tsys, WVR cal tables ########

#Examine Tsys
if(calplots):
  for asdm in basename_all:
    #Plotting Tsys vs. time
    plotcal(caltable=asdm+'.tsys.cal.fdm', 
      xaxis="time",yaxis="amp",
      spw='1:1200~1200',plotsymbol=".", subplot=421,
      antenna='0~7',
      iteration='antenna', figfile=asdm+'.tsys_vs_time.page1.png',
      fontsize=6.0)    
    plotcal(caltable=asdm+'.tsys.cal.fdm', 
      xaxis="time",yaxis="amp",
      antenna='8~15',
      spw='1:1200~1200',plotsymbol=".", subplot=421,
      iteration='antenna', figfile=asdm+'.tsys_vs_time.page2.png',
      fontsize=6.0)    
    #Plotting Tsys vs. frequency
    plotcal(caltable=asdm+'.tsys.cal.fdm',
      xaxis="freq",yaxis="amp",
      spw='1', plotsymbol=".", subplot=421,
      iteration='antenna', figfile=asdm+'.tsys_vs_freq.page1.png',
      antenna='0~7', fontsize=6.0)    
    plotcal(caltable=asdm+'.tsys.cal.fdm',
      xaxis="freq",yaxis="amp",
      spw='1', plotsymbol=".", subplot=421,
      iteration='antenna', figfile=asdm+'.tsys_vs_freq.page2.png',
      antenna='8~15', fontsize=6.0)    

#Apply Tsys, WVR calibrations.

field_names_north = ['Titan','3c279','J1127-1857','*Antennae*']
field_names_south = ['Titan','3c279','J1127-1857','*Antennae*']

for asdm in basename_all_north:
  for field in field_names_north:
    applycal(vis=asdm+".ms", spw='1', 
      field=field, gainfield=[field,field],
      interp='nearest', 
      gaintable=[asdm+".tsys.cal.fdm",asdm+'.wvr.cal'],
      flagbackup=F)

for asdm in basename_all_south:
  for field in field_names_south:
    applycal(vis=asdm+".ms", spw='1', 
      field=field, gainfield=[field,field],
      interp='nearest', 
      gaintable=[asdm+".tsys.cal.fdm",asdm+'.wvr.cal'],
      flagbackup=F)

# Check the cal weights here with plotms -> ok

# Split out corrected data and create summary file
for asdm in basename_all:
  split(vis=asdm+'.ms', outputvis=asdm+'.wvrtsys.ms', 
    datacolumn='corrected', spw='1', keepflags=F)    
  os.system('rm '+asdm+'.wvrtsys.listobs.txt')
  listobs(vis=asdm+'.wvrtsys.ms', listfile=asdm+'.wvrtsys.listobs.txt', verbose=True)

#   Plot spectral plot, amp and phase, 3c279
if(interact):
  for asdm in basename_all:
    plotms(vis=asdm+'.ms', 
      field='3c279',
      xaxis='frequency', yaxis='amp',
      selectdata=T, spw='1', 
      avgtime='1e8',avgscan=T,
      coloraxis='corr',
      iteraxis='baseline',
      antenna='*&*',
      ydatacolumn='data')
    dummy_string = raw_input("Next plot for "+asdm+" . Hit <Enter> to continue.")
    plotms(vis=asdm+'.ms', 
      field='3c279',
      xaxis='frequency', yaxis='phase',
      selectdata=T, spw='1', 
      avgtime='1e8',avgscan=T,
      coloraxis='corr',
      iteraxis='baseline',
      antenna='*&*',
      ydatacolumn='data')
    dummy_string = raw_input("Next plot for "+asdm+" . Hit <Enter> to continue.")

if(interact):
  for asdm in basename_all:
    # Plot continuum plot, amp and phase
    plotms(vis=asdm+'.ms', 
      field='0,1,2',
      xaxis='time', yaxis='amp',
      selectdata=T, spw='1:1200~1300', 
      avgchannel='1000',avgscan=F,
      coloraxis='corr',
      iteraxis='baseline',
      antenna='*&*',
      ydatacolumn='data')
    dummy_string = raw_input("Next plot for "+asdm+" . Hit <Enter> to continue.")
    plotms(vis=asdm+'.ms', 
      field='0,1,2',
      xaxis='time', yaxis='phase',
      selectdata=T, spw='1:1200~1300', 
      avgchannel='1000',avgscan=F,
      coloraxis='corr',
      iteraxis='baseline',
      antenna='*&*',
      ydatacolumn='data')
    dummy_string = raw_input("Next plot for "+asdm+" . Hit <Enter> to continue.")

######## Inspect and Flag data ########

for asdm in basename_all:
  # Unflag anything done in wrtsys.ms
  # flagdata(vis = asdm+'.wvrtsys.ms',mode='manualflag', unflag= T, flagbackup = F)
  # Flag the less sensitive edge channels of every data set
  flagdata(vis = asdm+'.wvrtsys.ms',spw = '0:0~7,0:3831~3839', flagbackup = F)
  # Possible contamination of Titan by Saturn
  flagdata(vis=asdm+'.wvrtsys.ms', mode='manualflag', field='Titan', uvrange='0~50', flagbackup = F)
  flagmanager(vis=asdm+'.wvrtsys.ms',mode='save',versionname ='User')

#Individual Flags

asdm='uid___A002_X20b4e9_X19'
flagdata(vis=asdm+'.wvrtsys.ms',antenna='DV05',spw='0',correlation='YY')  
flagdata(vis=asdm+'.wvrtsys.ms',antenna='DV12',spw='0:638~641;1663~1664;2430~2433',correlation='YY')

asdm='uid___A002_X20bcff_X1c'
flagdata(vis=asdm+'.wvrtsys.ms',antenna='DV05',spw='0',correlation='YY')  
flagdata(vis=asdm+'.wvrtsys.ms',antenna='DV12',spw='0:638~641;1663~1664;2430~2433',correlation='YY')

asdm='uid___A002_X21487a_X7e' 
flagdata(vis=asdm+'.wvrtsys.ms',antenna='DV13',spw='0',correlation='YY')
flagdata(vis=asdm+'.wvrtsys.ms',antenna='DV05',spw='0',correlation='YY') 

asdm='uid___A002_X219601_X13c'  
flagdata(vis=asdm+'.wvrtsys.ms',antenna='DV05',spw='0',correlation='YY')  

asdm='uid___A002_X21bf10_X1d8'
flagdata(vis=asdm+'.wvrtsys.ms',antenna='DV05',spw='0',correlation='YY') 

asdm='uid___A002_X20b4e9_X19'
flagdata(vis=asdm+'.wvrtsys.ms',antenna='DV05',spw='0',correlation='YY') 

asdm='uid___A002_X20bcff_X1c'
flagdata(vis=asdm+'.wvrtsys.ms',antenna='DV05',spw='0',correlation='YY') 

asdm='uid___A002_X22afc4_X24f'
flagdata(vis=asdm+'.wvrtsys.ms',antenna='DV05',spw='0',correlation='YY')

asdm='uid___A002_X22afc4_X3f4'
flagdata(vis=asdm+'.wvrtsys.ms',antenna='DV05',spw='0',correlation='YY')

######## Bandpass ########

for asdm in basename_all:
  # Running a short solution interval phase calibration
  os.system('rm -rf '+asdm+'.bpphase.gcal')
  gaincal(vis = asdm+'.wvrtsys.ms',
    selectdata=T,field = '3c279',spw = '0:1100~1300',
    caltable = asdm+'.bpphase.gcal',
    solint = 'int',refant = 'DV11',calmode='p')
  #Running a bandpass calibration
  os.system('rm -rf '+asdm+'.bandpass.bcal')
  bandpass(vis = asdm+'.wvrtsys.ms',
    field = '3c279',
    gaintable = asdm+'.bpphase.gcal',
    caltable = asdm+'.bandpass.bcal',
    bandtype='B',
    solint = 'inf',combine = 'scan,field', solnorm=T,refant = 'DV11',
    minblperant=3,minsnr=2,fillgaps=1)

#Plotting bandpass solutions
if(calplots):
  for asdm in basename_all:
    plotcal(caltable = asdm+'.bpphase.gcal',
      xaxis = 'time', yaxis = 'phase',
      iteration = 'antenna', plotrange=[0,0,-180,180],
      showgui=False, subplot=421, figfile=asdm+'.bpphase.page1.png',
      antenna='0~7')
    plotcal(caltable = asdm+'.bpphase.gcal',
      xaxis = 'time', yaxis = 'phase',
      iteration = 'antenna', plotrange=[0,0,-180,180],
      showgui=False, subplot=421, figfile=asdm+'.bpphase.page2.png',
      antenna='8~15')
    plotcal(caltable = asdm+'.bandpass.bcal', 
      xaxis = 'freq',yaxis = 'amp',
      plotrange = [0,0,0.8,1.2],
      antenna='0~7', iteration='antenna',
      showgui=False, subplot=421, figfile=asdm+'.bcal_amp.page1.png')
    plotcal(caltable = asdm+'.bandpass.bcal', 
      xaxis = 'freq',yaxis = 'amp',
      plotrange = [0,0,0.8,1.2],
      antenna='8~15', iteration='antenna',
      showgui=False, subplot=421, figfile=asdm+'.bcal_amp.page2.png')
    plotcal(caltable = asdm+'.bandpass.bcal', 
      xaxis = 'freq',yaxis = 'phase', iteration='antenna',
      antenna='0~7', subplot=421, figfile=asdm+'.bcal_phase.page1.png',
      plotrange = [0,0,-180,180], showgui=False)
    plotcal(caltable = asdm+'.bandpass.bcal', 
      xaxis = 'freq',yaxis = 'phase', iteration='antenna',
      antenna='8~15', subplot=421, figfile=asdm+'.bcal_phase.page2.png',
      plotrange = [0,0,-180,180], showgui=False)

#Reading model for Titan
for asdm in basename_all:
  setjy(vis = asdm+'.wvrtsys.ms',field = 'Titan',
    standard = 'Butler-JPL-Horizons 2010')

###Creating calibration solutions

#Carrying out short timescale (integration, ~ 10 sec) phase solution
for asdm in basename_all:
  os.system('rm -rf '+asdm+'.intphase.gcal')
  gaincal(vis=asdm+'.wvrtsys.ms',
    gaintable=asdm+'.bandpass.bcal', 
    caltable=asdm+'.intphase.gcal',
    calmode='p',
    field='Titan,3c279,J1127*',
    spw='0:40~3800',
    refant='DV11', solint='int',minsnr=2.0,minblperant=4)

#Carrying out longer timescale (by scan) phase solution
for asdm in basename_all:
  os.system('rm -rf '+asdm+'.scanphase.gcal')
  gaincal(vis=asdm+'.wvrtsys.ms',
    gaintable=asdm+'.bandpass.bcal', 
    caltable=asdm+'.scanphase.gcal',
    calmode='p',
    field='Titan,3c279,J1127*',
    spw='0:40~3800',
    refant='DV11', solint='inf',minsnr=2.0,minblperant=4)

#Solving for longer (scan) interval amplitude solution
for asdm in basename_all:
  os.system('rm -rf '+asdm+'.amp.cal')
  gaincal(vis = asdm+'.wvrtsys.ms',
    gaintable =[asdm+'.bandpass.bcal',asdm+'.intphase.gcal'],
    caltable = asdm+'.amp.cal',
    calmode='a',
    field = 'Titan,3c279,J1127*',
    spw='0:40~3800',
    refant = 'DV11',solint = 'inf')

#Scaling amplitude calibration to match Titan
for asdm in basename_all:
  os.system('rm -rf '+asdm+'.flux.cal')
  fluxscale(vis = asdm+'.wvrtsys.ms',
    caltable = asdm+'.amp.cal',
    fluxtable = asdm+'.flux.cal',
    reference = 'Titan',
    transfer = '3c279,J1127*')

# Plotting solutions
if(calplots):
  for asdm in basename_all:
    plotcal(caltable = asdm+'.scanphase.gcal',
      xaxis = 'time', yaxis = 'phase',
      iteration = 'antenna', plotrange=[0,0,-180,180],
      showgui=True, subplot=421, figfile=asdm+'.scanphase.page1.png',
      antenna='0~7', fontsize=6.0)
    plotcal(caltable = asdm+'.scanphase.gcal',
      xaxis = 'time', yaxis = 'phase',
      iteration = 'antenna', plotrange=[0,0,-180,180],
      showgui=False, subplot=421, figfile=asdm+'.scanphase.page2.png',
      antenna='8~15', fontsize=6.0)
    plotcal(caltable = asdm+'.flux.cal', 
      xaxis = 'time',yaxis = 'amp',
      plotrange = [0,0,0,0],
      antenna='0~7', iteration='antenna',
      showgui=False, subplot=421, figfile=asdm+'.flux.page1.png',
      fontsize=6.0)
    plotcal(caltable = asdm+'.flux.cal', 
      xaxis = 'time',yaxis = 'amp',
      plotrange = [0,0,0,0],
      antenna='8~15', iteration='antenna', fontsize=6.0,
      showgui=False, subplot=421, figfile=asdm+'.flux.page2.png')

#Applying calibrations to the different sources:
for asdm in basename_all_north:
  applycal(vis=asdm+'.wvrtsys.ms',field='3c279',
    gaintable=[asdm+'.bandpass.bcal',asdm+'.intphase.gcal',asdm+'.flux.cal'],
    interp=['nearest','nearest','nearest'],
    gainfield=['3c279','3c279','3c279'],flagbackup=F,calwt=F)
  applycal(vis=asdm+'.wvrtsys.ms',field='Titan',
    gaintable=[asdm+'.bandpass.bcal',asdm+'.intphase.gcal',asdm+'.flux.cal'],
    interp=['nearest','nearest','nearest'],
    gainfield=['3c279','Titan','Titan'],flagbackup=F,calwt=F)
  applycal(vis=asdm+'.wvrtsys.ms',field='J1127*',
    gaintable=[asdm+'.bandpass.bcal',asdm+'.intphase.gcal',asdm+'.flux.cal'],
    interp=['nearest','nearest','nearest'],
    gainfield=['3c279','J1127*','J1127*'],flagbackup=F,calwt=F)
  applycal(vis=asdm+'.wvrtsys.ms',field='*Antennae*',
    interp=['nearest','linear','linear'],
    gaintable=[asdm+'.bandpass.bcal',asdm+'.scanphase.gcal',asdm+'.flux.cal'],
    gainfield=['3c279','J1127*','J1127*'],flagbackup=F,calwt=F)

for asdm in basename_all_south:
  applycal(vis=asdm+'.wvrtsys.ms',field='3c279',
    gaintable=[asdm+'.bandpass.bcal',asdm+'.intphase.gcal',asdm+'.flux.cal'],
    interp=['nearest','nearest','nearest'],
    gainfield=['3c279','3c279','3c279'],flagbackup=F,calwt=F)
  applycal(vis=asdm+'.wvrtsys.ms',field='Titan',
    gaintable=[asdm+'.bandpass.bcal',asdm+'.intphase.gcal',asdm+'.flux.cal'],
    interp=['nearest','nearest','nearest'],
    gainfield=['3c279','Titan','Titan'],flagbackup=F,calwt=F)
  applycal(vis=asdm+'.wvrtsys.ms',field='J1127*',
    gaintable=[asdm+'.bandpass.bcal',asdm+'.intphase.gcal',asdm+'.flux.cal'],
    interp=['nearest','nearest','nearest'],
    gainfield=['3c279','J1127*','J1127*'],flagbackup=F,calwt=F)
  applycal(vis=asdm+'.wvrtsys.ms',field='*Antennae*',
    interp=['nearest','linear','linear'],
    gaintable=[asdm+'.bandpass.bcal',asdm+'.scanphase.gcal',asdm+'.flux.cal'],
    gainfield=['3c279','J1127*','J1127*'],flagbackup=F,calwt=F)

# Inspecting calibrated data (both quasars)
if(interact):
  for asdm in basename_all:
    plotms(vis = asdm+'.wvrtsys.ms', xaxis='uvdist', yaxis='amp',
      ydatacolumn='corrected', field='3c279',
      averagedata=True, avgchannel='3840', avgtime='',
      avgscan=F, avgbaseline=F, coloraxis='corr')
    plotms(vis = asdm+'.wvrtsys.ms', xaxis='time', yaxis='amp',
      ydatacolumn='corrected', field='3c279',
      averagedata=True, avgchannel='3840', avgtime='',
      avgscan=F, avgbaseline=F, coloraxis='corr')
    plotms(vis = asdm+'.wvrtsys.ms', xaxis='uvdist', yaxis='phase',
      ydatacolumn='corrected', field='3c279',
      avgchannel='3840', avgscan=F, avgbaseline=F, coloraxis='corr')
    plotms(vis = asdm+'.wvrtsys.ms', xaxis='time', yaxis='phase',
      ydatacolumn='corrected', field='3c279',
      avgchannel='3840', avgscan=F, avgbaseline=F, coloraxis='corr')

# Plot spectral plot, amp and phase
if(interact):
  for asdm in basename_all:
    plotms(vis= asdm+'.wvrtsys.ms', 
      field='3c279',
      xaxis='frequency', yaxis='amp',
      selectdata=T, spw='0', 
      avgtime='1e8',avgscan=T,
      coloraxis='corr',
      iteraxis='baseline',
      antenna='*&*',
      ydatacolumn='corrected')
    plotms(vis= asdm+'.wvrtsys.ms', 
      field='3c279',
      xaxis='frequency', yaxis='phase',
      selectdata=T, spw='0', 
      avgtime='1e8',avgscan=T,
      coloraxis='corr',
      iteraxis='baseline',
      antenna='*&*',
      ydatacolumn='corrected')

# Plot continuum plot, amp and phase, corrected
if(interact):
  for asdm in basename_all:
    plotms(vis= asdm+'.wvrtsys.ms', 
      field='0,1,2',
      xaxis='time', yaxis='amp',
      selectdata=T, spw='0:1200~1300', 
      avgchannel='1000',avgscan=F,
      coloraxis='corr',
      iteraxis='baseline',
      antenna='*&*',
      ydatacolumn='corrected')
    plotms(vis= asdm+'.wvrtsys.ms', 
      field='0,1,2',
      xaxis='time', yaxis='phase',
      selectdata=T, spw='0:1200~1300', 
      avgchannel='1000',avgscan=F,
      coloraxis='corr',
      iteraxis='baseline',
      antenna='*&*',
      ydatacolumn='corrected')

###### Split and concatenate #######

# Split, including spectral smoothing. width = 23 corresponds to ~15 km/s.
for asdm in basename_all_north:
  os.system('rm -rf '+asdm+'.cal.ms')
  split(vis = asdm+'.wvrtsys.ms',outputvis = asdm+'.cal.ms',
    field = '*Antennae*',spw='0',width=23, keepflags=False)
  os.system('rm '+asdm+'.cal.listobs.txt')
  listobs(asdm+'.cal.ms',listfile=asdm+'.cal.listobs.txt')

for asdm in basename_all_south:
  os.system('rm -rf '+asdm+'.cal.ms')
  split(vis = asdm+'.wvrtsys.ms',outputvis = asdm+'.cal.ms',
    field = '*Antennae*',spw='0',width=23, keepflags=False)
  os.system('rm '+asdm+'.cal.listobs.txt')
  listobs(asdm+'.cal.ms',listfile=asdm+'.cal.listobs.txt')

# Zero pointing table
for asdm in basename_all:
  tb.open(asdm+'.cal.ms/POINTING',nomodify=False)
  a = tb.rownumbers()
  tb.removerows(a)
  tb.close()

#Concat into two datasets, one for southern mosaic and another for the northern mosaic.

cal_south_vis = [vis+'.cal.ms' for vis in basename_all_south]
cal_north_vis = [vis+'.cal.ms' for vis in basename_all_north]

os.system('rm -rf Antennae_South.cal.ms')
concat(vis=cal_south_vis, concatvis='Antennae_South.cal.ms', timesort=T)

os.system('rm -rf Antennae_North.cal.ms')
concat(vis=cal_north_vis, concatvis='Antennae_North.cal.ms', timesort=T)