#ALMA moment mapping script created by George J. Bendo for the ALMA Data
#Processing Workshop at UCLan (23-24 January 2019).  This is in part based
#on scripts created by Anita Richards for the UK ALMA Workshop in Durham 
#(14-16 November 2018).
#
#This is a custom-written script created for use with calibrated
#visibility data of RXCJ0439+05.  The script cound be executed using
#the command execfile('RXCJ0439_selfcalibration.py') at the command
#line, or the individual lines could be copied and pasted into a CASA
#session.



#Run listobs.
os.system('rm -rf RXCJ0439.split.listobs')
listobs(vis='RXCJ0439.split',listfile='RXCJ0439.split.listobs')

#Remove any existing calibration tables.
clearcal(vis='RXCJ0439.split')



#Create an initial image.
os.system('rm -rf RXCJ0439_selfcal0.cont.*')
tclean(vis = 'RXCJ0439.split',
  imagename = 'RXCJ0439_selfcal0.cont',
  field = '0',
  stokes = 'I',
  spw = '0,1,2,3:0~44;70~119',
  outframe = 'LSRK',
  specmode = 'cont',
  nterms = 2,
  imsize = [300, 300],
  cell = '0.05arcsec',
  deconvolver = 'mtmfs',
  niter = 500,
  weighting = 'briggs',
  robust = 0.5,
  gridder = 'standard',
  interactive = True,
  )



#Insert the model from tclean into the measurement set.
ft(vis='RXCJ0439.split',
  model=['RXCJ0439_selfcal0.cont.model.tt0','RXCJ0439_selfcal0.cont.model.tt1'],
  nterms=2,
  field='0',
  usescratch=True)

#Derive a phase calibration table based on this model.
os.system('rm -rf RXCJ0439.split.selfcal.phase1')
gaincal(vis='RXCJ0439.split',
  caltable='RXCJ0439.split.selfcal.phase1',
  field='0',
  gaintype='T',
  refant='DA47',
  calmode='p',
  solint='inf')

#Plot the phase calibration table.
plotcal(caltable='RXCJ0439.split.selfcal.phase1',
  xaxis='time',
  yaxis='phase',
  timerange='',
  iteration='antenna',
  subplot=421,
  plotrange=[0,0,-180,180])

#Apply the calibration.
applycal(vis='RXCJ0439.split',
  field='0',
  gaintable=['RXCJ0439.split.selfcal.phase1'],
  gainfield=['0'],
  calwt=False,
  flagbackup=False,
  applymode='calflag')



#Image the data after the phase self-calibration.
os.system('rm -rf RXCJ0439_selfcal1.cont.*')
tclean(vis = 'RXCJ0439.split',
  imagename = 'RXCJ0439_selfcal1.cont',
  field = '0',
  stokes = 'I',
  spw = '0,1,2,3:0~44;70~119',
  outframe = 'LSRK',
  specmode = 'cont',
  nterms = 2,
  imsize = [300, 300],
  cell = '0.05arcsec',
  deconvolver = 'mtmfs',
  niter = 500,
  weighting = 'briggs',
  robust = 0.5,
  gridder = 'standard',
  interactive = True,
  )



#Insert the model from tclean into the measurement set.
ft(vis='RXCJ0439.split',
  model=['RXCJ0439_selfcal1.cont.model.tt0','RXCJ0439_selfcal1.cont.model.tt1'],
  field='0',
  nterms=2,
  usescratch=True)

#Derive an amplitude calibration table based on this model.
os.system('rm -rf RXCJ0439.split.selfcal.amp1')
gaincal(vis='RXCJ0439.split',
  caltable='RXCJ0439.split.selfcal.amp1',
  field='0',
  gaintype='T',
  refant='DA47',
  calmode='a',
  solint='inf',
  gaintable='RXCJ0439.split.selfcal.phase1')

#Plot the amplitude calibration table.
plotcal(caltable='RXCJ0439.split.selfcal.amp1',
  xaxis='time',
  yaxis='amp',
  timerange='',
  iteration='antenna',
  subplot=421)

#Apply the calibration.  Include the phase calibration as well.
applycal(vis='RXCJ0439.split',
  field='0',
  gaintable=['RXCJ0439.split.selfcal.phase1','RXCJ0439.split.selfcal.amp1'],
  gainfield=['0','0'],
  calwt=False,
  flagbackup=False,
  applymode='calflag')



#Image the data after the phase and amplitude self-calibration.
os.system('rm -rf RXCJ0439_selfcal2.cont.*')
tclean(vis = 'RXCJ0439.split',
  imagename = 'RXCJ0439_selfcal2.cont',
  field = '0',
  stokes = 'I',
  spw = '0,1,2,3:0~44;70~119',
  outframe = 'LSRK',
  specmode = 'cont',
  nterms = 2,
  imsize = [300, 300],
  cell = '0.05arcsec',
  deconvolver = 'mtmfs',
  niter = 500,
  weighting = 'briggs',
  robust = 0.5,
  gridder = 'standard',
  interactive = True,
  )