########################################################################
# Data reduction script for HD163296 Band 7:
# - Calibration script: "HD163296_Band7_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

#----------------------------------------------------------------------------------
#----- Some setup steps -----------------------------------------------------------
#----------------------------------------------------------------------------------

version = casalog.version()
print "You are using " + version
if (int(version.split()[4][1:-1]) < 19874):
    print "\033[91m YOUR VERSION OF CASA IS TOO OLD FOR THIS GUIDE."
    print "\033[91m PLEASE UPDATE IT BEFORE PROCEEDING."
else:
    print "Your version of CASA is appropriate for this guide."


#----------------------------------------------------------------------------------
#----- Input calibrated dataset ---------------------------------------------------
#----------------------------------------------------------------------------------

name='HD163296_Band7_concat'

#----------------------------------------------------------------------------------
#----- Image the continuum --------------------------------------------------------
#----------------------------------------------------------------------------------

os.system('rm -rf HD163296_continuum*') 
clean(vis=name+'.ms',
  imagename= 'HD163296_continuum',
  field='0',
  spw  = '0,1:1~1700;2200~3800,2:1~1500;2400~3800,3',
  mode = 'mfs',
  niter = 1000,
  psfmode = 'clark',
  interactive = True,
  threshold = '0.1mJy',  
  npercycle=100,
  cell = '0.05arcsec',
  imsize = 512,
  weighting = 'briggs')


#----------------------------------------------------------------------------------
#----- SELF-CALIBRATION -----------------------------------------------------------
#----------------------------------------------------------------------------------

os.system('rm -rf gaincal_p.sc1') 
gaincal(vis=name+'.ms',
  caltable='gaincal_p.sc1',
  solint='int', 
  refant='DV04',
  spw  = '0,1:1~1700;2200~3800,2:1~1500;2400~3800,3',
  gaintype = 'T',
  minsnr = 2,      
  calmode='p')

if(calplots):
    plotcal(caltable='gaincal_p.sc1',
            xaxis='time',
            yaxis='phase',
            iteration='spw,antenna',
            subplot=421,
        plotrange=[0,0,-180,180])

applycal(vis=name+'.ms',
  gaintable=['gaincal_p.sc1'],
  calwt=F)
 
clean (vis=name+'.ms',
  imagename= 'HD163296_continuum_sc1',
  field='0',
  spw  = '0,1:1~1700;2200~3800,2:1~1500;2400~3800,3',
  mode = 'mfs',
  niter = 1000,
  psfmode = 'clark',
  interactive = True,   
  threshold = '0.1mJy',  
  npercycle=100,
  cell = '0.05arcsec',
  imsize = 512,
  weighting = 'briggs')


#----- Selfcal 2
os.system('rm -rf gaincal_p.sc2') 
gaincal(vis=name+'.ms',
  caltable='gaincal_p.sc2',
  solint='int', 
  refant='DV04',
  spw  = '0,1:1~1700;2200~3800,2:1~1500;2400~3800,3',
  gaintype = 'T',
  minsnr = 2,      
  calmode='p')

if(calplots):
    plotcal(caltable='gaincal_p.sc2',
            xaxis='time',
            yaxis='phase',
            iteration='spw,antenna',
            subplot=421,
        plotrange=[0,0,-180,180])

applycal(vis=name+'.ms',
  gaintable=['gaincal_p.sc2'],
  calwt=F)

clean (vis=name+'.ms',
  imagename= 'HD163296_continuum_sc2',
  field='0',
  spw  = '0,1:1~1700;2200~3800,2:1~1500;2400~3800,3',
  mode = 'mfs',
  niter = 1000,
  psfmode = 'clark',
  interactive = True,     
  threshold = '0.1mJy',  
  npercycle=200,
  cell = '0.05arcsec',
  imsize = 512,
  weighting = 'briggs')


#----- Selfcal 3 (amplitude & phase)
os.system('rm -rf gaincal_ap.sc3') 
gaincal(vis=name+'.ms',
  caltable='gaincal_ap.sc3',
  gaintable='gaincal_p.sc2',      
  solint='inf', 
  refant='DV04',
  spw  = '0,1:1~1700;2200~3800,2:1~1500;2400~3800,3',
  gaintype = 'T',
  minsnr = 2,      
  calmode='ap')

if(calplots):
    plotcal(caltable='gaincal_ap.sc3',
            xaxis='time',
            yaxis='phase',
            iteration='spw,antenna',
            subplot=421,
        plotrange=[0,0,-180,180])

    plotcal(caltable='gaincal_ap.sc3',
            xaxis='time',
            yaxis='amp',
            iteration='spw,antenna',
            subplot=421,
        plotrange=[0,0,-0.6,1.0])

applycal(vis=name+'.ms',
  gaintable=['gaincal_p.sc2','gaincal_ap.sc3'],
  calwt=F)

clean (vis=name+'.ms',
  imagename= 'HD163296_continuum_sc3',
  field='0',
  spw  = '0,1:1~1700;2200~3800,2:1~1500;2400~3800,3',
  mode = 'mfs',
  niter = 1000,
  psfmode = 'clark',
  interactive = True,     
  threshold = '0.1mJy',  
  npercycle=200,
  cell = '0.05arcsec',
  imsize = 512,
  weighting = 'briggs')

#----------------------------------------------------------------------------------
#----- SPLIT the sc-data ----------------------------------------------------------
#----------------------------------------------------------------------------------

os.system('rm -rf HD163296_calibrated*') 
split(vis = name+'.ms',
    outputvis = 'HD163296_calibrated.sc',
    datacolumn = 'corrected',    
    keepflags = T)

#----------------------------------------------------------------------------------
#----- Continuum substraction -----------------------------------------------------
#----------------------------------------------------------------------------------

uvcontsub2(vis='HD163296_calibrated.sc',
  field='0',
  fitspw='0,1:300~1800;2200~3800,2:100~1500;2400~3800,3',
  combine='spw',
  solint='int',
  fitorder=1)

#----------------------------------------------------------------------------------
#----- Spectral Line Imaging: CO(3-2) ---------------------------------------------
#----------------------------------------------------------------------------------

os.system('rm -rf HD163296.CO3-2Line*') 
clean(vis = 'HD163296_calibrated.sc.contsub',
  imagename = 'HD163296.CO3-2Line.Clean',
  field='0',
  spw='2',
  mode = 'velocity',
  nchan = 140,  
  start = '-2.5km/s',
  width='0.110km/s',
  niter = 1000,
  psfmode = 'hogbom',
  interactive = False,
  mask = 'HD163296_CO3-2_R2_masked.mask',
  robust = 2,    
  npercycle = 100,
  cell = '0.05arcsec',
  imsize = 512,
  weighting = 'briggs',
  restfreq = '345.795990GHz',
  outframe = 'LSRK')


# HCO+(4-3) map 
clean(vis = 'HD163296_calibrated.sc.contsub',
  imagename = 'HD163296_HCO+',
  field='0',
  spw='1',
  mode = 'velocity',
  nchan = 180,  
  start = '-2.0km/s',
  width='0.110km/s',
  niter = 1000,
  psfmode = 'hogbom',
  interactive = False,
  mask = '',
  npercycle = 100,
  cell = '0.05arcsec',
  imsize = 512,
  weighting = 'briggs',
  restfreq = '356.7341340GHz',
  outframe = 'LSRK')


#----- Create Moment Maps for CO(3-2)

immoments(imagename='HD163296.CO3-2Line.Clean.image', 
  moments=[0],
  chans='', 
  box='', 
  axis='spectral',
  includepix=[0.01, 10000], 
  outfile='HD163296.CO3-2.Clean.image.mom0')

immoments(imagename='HD163296.CO3-2Line.Clean.image', 
  moments=[1],
  chans='', 
  box='', 
  axis='spectral',
  includepix=[0.3, 10000], 
  outfile='HD163296.CO3-2Line.Clean.image.mom1')

os.system('rm -f HD163296.CO3-2Line.Clean.image.mom0.png')
imview (raster={'file': 'HD163296.CO3-2Line.Clean.image.mom0',
                'colorwedge':T,'scaling': -0.5,'range':[-3,30]},
         zoom=1,out='HD163296.CO3-2Line.Clean.image.mom0.png')

os.system('rm -f HD163296.CO3-2Line.Clean.image.mom1.png')
imview (raster={'file': 'HD163296.CO3-2Line.Clean.image.mom1',
                 'colorwedge':T},
         zoom=1,out='HD163296.CO3-2Line.Clean.image.mom1.png')

#----- Export as fits
exportfits(imagename='HD163296_continuum_sc3.image',
     fitsimage='HD163296_continuum_sc3.image.fits')
exportfits(imagename='HD163296.CO3-2Line.Clean.image',
     fitsimage='HD163296.CO3-2Line.Clean.image.fits')
exportfits(imagename='HD163296.CO3-2Line.Clean.image.mom0',
     fitsimage='HD163296.CO3-2Line.Clean.image.mom0.fits')
exportfits(imagename='HD163296.CO3-2Line.Clean.image.mom1',
     fitsimage='HD163296.CO3-2Line.Clean.image.mom1.fits')

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

96.CO3-2Line.Clean.image.mom1.fits')

exportfits(imagename='HD163296.CO3-2Line.Clean.image.mom1',
     fitsimage='HD163296.CO3-2Line.Clean.image.mom1.fits')

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