print "# Now begins the main calibration"

# DATA INSPECTION AND FLAGGING:
# Note: Atmospheric lines in spw 3, observed in all sources.

print "# Please, do DATA INSPECTION using plotms. Below we report the bad data found and we flag them"

print "# Flag the first 15 channels in all spws."

flagdata(vis = 'uid___A002_X24c120_X18c.ms.split',
 mode='manualflag',
 spw='*:0~15',
 flagbackup = F)

print "# Flag the last 8 channels in all spws."

flagdata(vis = 'uid___A002_X24c120_X18c.ms.split',
 mode='manualflag',
 spw='*:120~128',
 flagbackup = F)

print "# Flag antenna DV12 in polarization XX." 

flagdata(vis = 'uid___A002_X24c120_X18c.ms.split',
 mode='manualflag',
 spw='3',
 antenna='DV12',
 correlation='XX',
 flagbackup = F)


print "# Running fixplanets on fields with 0,0 coordinates."

fixplanets(vis = 'uid___A002_X24c120_X18c.ms.split',
  field = 'Titan',
  fixuvw = T)

print "# Putting a model for the spatially-resolved calibrators."

setjy(vis = 'uid___A002_X24c120_X18c.ms.split',
  field = 'Titan ',
  standard = 'Butler-JPL-Horizons 2010')


print "# Saving flag column to version BeforeBandpassCalibration"

flagmanager(vis = 'uid___A002_X24c120_X18c.ms.split',
  mode = 'save',
  versionname = 'BeforeBandpassCalibration')


print "# BANDPASS CALIBRATION:"


gaincal(vis = 'uid___A002_X24c120_X18c.ms.split',
  caltable = 'uid___A002_X24c120_X18c.ms.split.pre_bandpass',
  field = 'J1517-243',
  refant = 'DV11',
  calmode = 'ap',
  spw = '*:40~70',
  solint = 'int')


print "# Checking gain solutions."

plotcal(caltable = 'uid___A002_X24c120_X18c.ms.split.pre_bandpass',
 xaxis = 'time',
 yaxis = 'amp',
 subplot = 221,
 iteration = 'spw,antenna',
 plotrange = [0,0,0,1])

raw_input ("Please press NEXT button in the CASA plotter pop-up window and then press enter to continue the script.")


print "# Checking gain solutions."

plotcal(caltable = 'uid___A002_X24c120_X18c.ms.split.pre_bandpass',
 xaxis = 'time',
 yaxis = 'phase',
 subplot = 221,
 iteration = 'spw,antenna',
 plotrange = [0,0,-180,180])


raw_input ("Please press NEXT button in the CASA plotter pop-up window and then press enter to continue the script.")

bandpass(vis = 'uid___A002_X24c120_X18c.ms.split',
  caltable = 'uid___A002_X24c120_X18c.ms.split.bandpass',
  field = 'J1517-243',
  solint = 'inf',
  combine='scan',
  refant = 'DV11',
  solnorm = T,
  bandtype = 'B',
  gaintable = 'uid___A002_X24c120_X18c.ms.split.pre_bandpass')


print "# Checking bandpass solutions."

plotcal(caltable = 'uid___A002_X24c120_X18c.ms.split.bandpass',
 xaxis = 'chan',
 yaxis = 'amp',
 subplot = 221,
 plotsymbol = ',',
 iteration = 'spw,antenna',
 plotrange = [0,0,0.85, 1.10])

raw_input ("Please press NEXT button in the CASA plotter pop-up window and then press enter to continue the script.")

print "# Checking bandpass solutions."

plotcal(caltable = 'uid___A002_X24c120_X18c.ms.split.bandpass',
 xaxis = 'chan',
 yaxis = 'phase',
 subplot = 221,
 plotsymbol = ',',
 iteration = 'spw,antenna',
 plotrange = [0,0,-180, 180])

raw_input ("Please press NEXT button in the CASA plotter pop-up window and then press enter to continue the script.")

print "# GAIN CALIBRATION:"

print "# First let's look for outliers."

gaincal(vis = 'uid___A002_X24c120_X18c.ms.split',
 caltable = 'uid___A002_X24c120_X18c.ms.split.gaincal0',
 gaintable = 'uid___A002_X24c120_X18c.ms.split.bandpass',
 field = 'Titan,J1604-446',
 refant = 'DV11',
 solint = '10s')

print "# Plotting amplitude as a function of time to look for outliers."

plotcal(caltable = 'uid___A002_X24c120_X18c.ms.split.gaincal0',
 xaxis = 'time',
 yaxis = 'amp',
 subplot = 221,
 iteration = 'spw,antenna')

raw_input ("Please press NEXT button in the CASA plotter pop-up window and then press enter to continue the script.")

# looks O.K. no outliers

print "# Getting gain calibration for Titan and the phase calibrator."

gaincal(vis = 'uid___A002_X24c120_X18c.ms.split',
 caltable = 'uid___A002_X24c120_X18c.ms.split.gaincal',
 gaintable = 'uid___A002_X24c120_X18c.ms.split.bandpass',
 field = 'Titan,J1604-446',
 refant = 'DV11',
 solint = 'inf')

print "# Checking the gaincal solutions."

plotcal(caltable = 'uid___A002_X24c120_X18c.ms.split.gaincal',
 xaxis = 'time',
 yaxis = 'amp',
 subplot = 221,
 iteration = 'spw,antenna')

raw_input ("Please press NEXT button in the CASA plotter pop-up window and then press enter to continue the script.")


plotcal(caltable = 'uid___A002_X24c120_X18c.ms.split.gaincal',
 xaxis = 'time',
 yaxis = 'phase',
 subplot = 221,
 iteration = 'spw,antenna')


print "# Fluxscale using Titan."

fluxscale(vis = 'uid___A002_X24c120_X18c.ms.split',
 caltable = 'uid___A002_X24c120_X18c.ms.split.gaincal',
 fluxtable = 'uid___A002_X24c120_X18c.ms.split.fluxscale',
 reference = 'Titan',
 transfer = '')

print "# Application of the bandpass and gain calibration tables with the correct flux scaling."


applycal(vis = 'uid___A002_X24c120_X18c.ms.split',
 flagbackup=F,
 field='3,4',   # Phase calibrator and science target
 gainfield=['2','3'], # Bandpass calibrator and phase calibrator
 gaintable = ['uid___A002_X24c120_X18c.ms.split.bandpass','uid___A002_X24c120_X18c.ms.split.fluxscale'])


print "#  Let's produce a calibrator image to check the calibration."


for i in ['0','1','2','3']:
  clean(vis = 'uid___A002_X24c120_X18c.ms.split',
    cell = '0.2arcsec',
    imsize = 256,
    field = 'J1604-446',
    niter = 100,
    psfmode='hogbom',
    spw = i,
    imagename = 'J1604-446_check_cal_sp_'+i)

 
for i in ['0','1','2','3']:
  viewer(infile='J1604-446_check_cal_sp_'+i+'.image',
    displaytype='raster')


print "# Splitting the science data."

split(vis = 'uid___A002_X24c120_X18c.ms.split',
 outputvis = 'GRB110715A_calibrated.ms',
 datacolumn = 'corrected',
 field = '4',
 timebin = '')

print "# Done calibration of science target. Now let's produce an image."

clean(vis = 'GRB110715A_calibrated.ms',
 imagename = 'GRB110715A',
 imsize=512,
 cell = '0.05arcsec',
 mode='mfs',
 mask='GRB110715A.mask',
 spw='0,1,2,3',
 weighting = 'briggs',
 psfmode='hogbom',
 robust=2,
 niter = 500,
 interactive = F)


viewer(infile='GRB110715A.image',
  displaytype='raster')