########################################################################
# Data reduction script for SgrA* Band 3:
# - Imaging script: "SgrA_Band3_Imaging.py" -
# Tested in CASA Version 3.4.0 (r19988)
########################################################################
"""
See accompanying README file for details of the necessary input files
and comments on the data

"""
#----------------------------------------------------------------------------------
#----- Optional Steps -------------------------------------------------------------
#----------------------------------------------------------------------------------

#----- Set this option to true if you want to make diagnostic plots along the way.
#----- Note that this will slow down the reduction significantly.
#----- 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


#----- Make clean continuum image (all spw!!!!!)
# We use SgraB3.ms before splitting which has the 4 spw
os.system('rm -rf SgrA.image.c*')
clean(vis='SgrA.ms',
      cell='0.80arcsec',
      imsize=128,
      field='SgrA',
      interactive=T,
      threshold='30mJy',
      selectdata=T,
      mode='mfs',
      niter=10000,
      npercycle=50,
      imagename='SgrA.image.c',
      spw='0:151~1100;1600~2500;3000~3699,1:151~3000,2:151~3699,3:151~1100;2500~3699')

exportfits(imagename = 'SgrA.image.c.image',
	fitsimage = 'SgrA.image.c.fits')
exportfits(imagename = 'SgrA.image.c.psf',
	fitsimage = 'SgrA.image.c.psf.fits')

#------ IMAGING OF H39a at 106.7373. Using the spw3 splitted data only.
#------ Subtract Continuum in spw 3
uvcontsub2(vis = 'SgrA.mssplitSpw3',
	field='SgrA',
	fitspw='0:300~1700;2400~3500',
	solint = '120s',
	fitorder=1,
	spw='0:1700~2400')

#------ Spectral line plotting: The line is not evident when averaging all baselines so we plot baseline 0&4.
if(interact):
   plotms(vis = 'SgrA.mssplitSpw3.contsub',
	xaxis = 'chan',yaxis = 'amp',
	avgtime = '1e8',
	avgbaseline = T,
	avgscan = T,
	field='SgrA',
	overwrite=T,
	antenna='0&4',
	avgchannel='5',
	transform=T,
	freqframe='LSRK',
	restfreq='106.7373GHz',
	plotfile = 'SgrA.mssplitSpw3.contsub.H39a.line.png')

#------ Cleaning the line cube
os.system('rm -rf SgrA.image.line*')

clean(vis='SgrA.mssplitSpw3.contsub',
      imagename='SgrA.image.line',
      cell='0.80arcsec',
      imsize=128,
      field='SgrA',
      outframe='LSRK',
      restfreq='106.7373GHz',
      interactive=T,
      threshold='1.0mJy',
      selectdata=T,
      spw='0',
      mode='channel',
      interpolation='nearest',
      width=5,
      niter=1000,
      npercycle=100,
      psfmode = 'hogbom',
      robust = 2,
      weighting = 'briggs',
      usescratch=False)


exportfits(imagename = 'SgrA.image.line.image',
	fitsimage = 'SgrA.image.line.image.fits')
exportfits(imagename = 'SgrA.image.line.psf',
	fitsimage = 'SgrA.image.line.psf.fits')

#----- Creating moment maps
immoments(imagename = 'SgrA.image.line.image',
	  outfile='SgrA.image.line',
	  moments=[0,1],
	  axis='spectral',
	excludepix=[-0.021,0.021])

exportfits(imagename = 'SgrA.image.line.integrated',
	fitsimage = 'SgrA.image.line.integrated.fits')
exportfits(imagename = 'SgrA.image.line.weighted_coord',
	fitsimage = 'SgrA.image.line.weighted_coord.fits')

#----------------------------------------------------------------------------------
#----- End of imaging script.
#----------------------------------------------------------------------------------