########################################################################
# Data reduction script for HD163296 Band 6:
# - Calibration script: "HD163296_Band6_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_Band6.Calib'


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

#spw 1: avoid 3 first channels that have the line
#spw 2: avoid 3 first channels that have the line

clean(vis = name+'.ms',
  imagename = 'HD163296_Band6.Cont',
  field = 'HD163296',
  spw = '0,1:0~107;188~3563;3620~3839,2:0~119;220~3839,3',
  mode = 'mfs',
  niter=1000,
  interactive = T,
  imsize = [432, 432],
  cell = '0.05arcsec',
  weighting = 'briggs',
  phasecenter='HD163296',
  robust = 0.5,
  timerange='')

viewer('HD163296_Band6.Cont.image')

imfit(imagename='HD163296_Band6.Cont.image',
      logfile='HD163296_Band6.Cont.image.imfit.txt')


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

os.system('rm -rf hd163296_selfcal-phase')

gaincal(vis = name+'.ms',
        refant='DV04',
        calmode='p',
        caltable='hd163296_selfcal-phase',
    solint='120s')

if(calplots):
    plotcal(caltable = 'hd163296_selfcal-phase',
            xaxis = 'time',
            yaxis  = 'phase',
            iteration='antenna',
        plotrange=[0,0,-180,180])

applycal(vis = name+'.ms',
         gaintable='hd163296_selfcal-phase')

#----- clean again

clean(vis = name+'.ms',
      imagename = name+'.Cont.sc1',
      field = 'HD163296',
      spw = '0,1:0~107;188~3563;3620~3839,2:0~119;220~3839,3',
      mode = 'mfs',
      niter=1000,
      interactive = T,
      imsize = [432, 432],
      cell = '0.05arcsec',
      weighting = 'briggs',
      robust = 0.5,
    timerange='')

viewer(name+'.Cont.sc1.image')

imfit(imagename = name+'.Cont.sc1.image',
      logfile=name+'.Cont.sc1.image.imfit.txt')


#----- Second Iteration of self-cal now amp and phase

gaincal(vis = name+'.ms',
        refant='DV04',
        calmode='ap',
        caltable='hd163296_selfcal-amp-phase',
        solint='inf',
    gaintable = 'hd163296_selfcal-phase')

if(calplots):
    plotcal(caltable = 'hd163296_selfcal-amp-phase',
        xaxis =  'time',
        yaxis =  'amp',
        iteration='antenna',
        plotrange=[0,0,0.5,1.5],
        iteration='antenna',
        subplot=541,
        figfile='test-plotcal')

applycal(vis = name+'.ms',
         gaintable=['hd163296_selfcal-phase','hd163296_selfcal-amp-phase'])


#----- clean again

clean(vis = name+'.ms',
  imagename = name+'Cont.sc2',
  field = 'HD163296',
  spw = '0,1:0~107;188~3563;3620~3839,2:0~119;220~3839,3',
  mode = 'mfs',
  niter=1000,
  interactive = T,
  imsize = [432, 432],
  cell = '0.05arcsec',
  weighting = 'briggs',
  robust = 0.5,
  timerange='')

viewer(name+'Cont.sc2.image')

imfit(imagename = name+'Cont.sc2.image',
      logfile='hd163296_cont_spw_1234_selfcal2.image.imfit.txt')


#----------------------------------------------------------------------------------
#----- CO(2-1) Line Imaging -----------------------------------------------------------
#----------------------------------------------------------------------------------

linename='HD163296.CO2-1Line'

#----- split

os.system('rm -rf HD163296.CO2-1Line*')

split(vis=name+'.ms',
    outputvis=linename,
    datacolumn = 'corrected',
    field='HD163296',
    spw='2:0~1000')


#----- plot the line

if(interact):
plotms(vis = linename+'.ms',
       field = '0',
       xaxis = 'channel',
       yaxis = 'amp',
       ydatacolumn = 'corrected',
       spw = '0',
       iteraxis = 'spw',
       xselfscale = T,
       yselfscale = T,
       coloraxis = 'corr',
       avgchannel = '',
       avgtime = '30000',
       avgbaseline = T,
       avgscan = T,
       transform=T,
    freqframe='LSRK')


#----- Continuum subtraction 


uvcontsub2(vis = linename,
           fitspw = '0:0~50,0:400~999',
    fitorder = 1,spw = '0')

os.system('cp -r HD163296.CO2-1Line.contsub HD163296.CO2-1Line.spw2.contsub')

listobs(linename+'.contsub')


#----- Cvel 

os.system('rm -rf HD163296.CO2-1Line.contsub.vslr')

cvel(vis = linename+'.contsub',
     outputvis = linename+'.contsub.vlsr',
     field='0',
     spw='0',
     selectdata  =  F,
     mode  =  'channel',
     nchan = 250,
     start = 50,
     width = 1,
     interpolation  = 'linear',
     phasecenter =  0,
     restfreq  = '230.538GHz',
     outframe  =  'lsrk',
     veltype  = 'radio',
    hanning = F,
    async = F)


#----- Clean

os.system('rm -rf HD163296_CO_2_1*')

clean(vis = linename+'.contsub.vlsr',
      imagename = 'HD163296_CO_2_1',
      field = 'HD163296',
      phasecenter = 0,
      cell = '0.05arcsec',
      imsize = [432,432],
      spw = '0',
      mode = 'channel',
      start = 50,
      width = 1,
      nchan = 250,
      restfreq='230.538GHz',
      weighting = 'briggs',
      niter = 1500,
      interactive = T,
      imagermode = 'csclean',
      cyclefactor = 5)

viewer('HD163296_CO_2_1.image')


#----------------------------------------------------------------------------------
#----- 13CO Line Imaging ----------------------------------------------------------
#----------------------------------------------------------------------------------

#----- split 13CO line

linename='HD163296.13COLine'

os.system('rm -rf HD163296.13COLine*')

split(vis=name+'.ms',
      outputvis=linename,
    datacolumn = 'corrected',
    field='HD163296',
    spw='1:0~1000')

if(interact):
    plotms(vis =linename+'.ms',
           field = '0',
           xaxis = 'channel',
           yaxis = 'amp',
           ydatacolumn = 'corrected',
           spw = '0',
           iteraxis = 'spw',
           xselfscale = T; yselfscale = T,
           coloraxis = 'corr',
           avgchannel = '',
           avgtime = '30000',
           avgbaseline = T,
           avgscan = T,
           transform=T,
           freqframe='LSRK')


#----- Continuum subtraction

uvcontsub2(vis = linename,
           fitspw = '0:0~50,0:400~999',
    fitorder = 1,
    spw = '0')


listobs(linename+'.contsub')

if(interact):
    plotms(vis = linename+'.contsub',
           field = '0',
           xaxis = 'channel',
           yaxis = 'amp',
           ydatacolumn = 'data',
           spw = '0',
           iteraxis = 'spw',
           xselfscale = T,
           yselfscale = T,
           coloraxis = 'corr',
           avgchannel = '',
           avgtime = '30000',
           avgbaseline = T,
           avgscan = T,
           uvrange='0~100m',
           transform=T,
           freqframe='LSRK')


#----- Cvel

os.system('rm -rf HD163296.13COLine.contsub.vslr')

cvel(vis = linename+'.contsub',
     outputvis = linename+'.contsub.vlsr',
     field='0',
     spw='0',
     selectdata  =  F,
     mode  =  'channel',
     nchan = 110,
     start = 100,
     width = 1,
     interpolation  = 'linear',
     phasecenter =  0,
     restfreq  = '220.39868GHz',
     outframe  =  'lsrk',
     veltype  = 'radio',
    hanning = F,
    async = F)

listobs(linename+'.contsub.vlsr')


#----- clean

os.system('rm -rf HD163296_13CO*')

clean(vis = linename+'.contsub.vlsr',
      imagename = 'HD163296_13CO',
      field = 'HD163296',
      phasecenter = 0,
      cell = '0.05arcsec',
      imsize = [432,432],
      spw = '0',
      mode = 'channel',
      start = 1,
      width = 1,
      nchan = 110,
      restfreq='220.39868GHz',
      weighting = 'briggs',
      niter = 1500,
      interactive = T,
      imagermode = 'csclean',
      cyclefactor = 5)

viewer('HD163296_13CO.image')


#----------------------------------------------------------------------------------
#----- C18O Imaging ---------------------------------------------------------------
#----------------------------------------------------------------------------------

linename='HD163296.C18OLine'

os.system('rm -rf HD163296.C18OLine*')


#----- split C18O line

split(vis=name+'.ms',
      outputvis=linename,
      datacolumn = 'corrected',
      field='HD163296',
    spw='1:2800~3800')


#----- plot the line

if(interact):
    plotms(vis=linename+'.ms',
       field = '0',
       xaxis = 'channel',
       yaxis = 'amp',
       ydatacolumn = 'corrected',
       spw = '0',
       iteraxis = 'spw',
       xselfscale = T,
       yselfscale = T,
       coloraxis = 'corr',
       avgchannel = '',
       avgtime = '30000',
       avgbaseline = T,
       avgscan = T,
       transform=T,
       freqframe='LSRK')


#----- Continuum subtraction

uvcontsub2(vis = linename,
           fitspw = '0:0~650,0:920~980',
           fitorder = 1,
           spw = '0')


listobs(linename+'.contsub')

if(interact):
    plotms(vis =linename+'.contsub',
       field = '0',
       xaxis = 'channel',
       yaxis = 'amp',
       ydatacolumn = 'data',
       spw = '0',
       iteraxis = 'spw',
       xselfscale = T; yselfscale = T,
       coloraxis = 'corr',
       avgchannel = '',
       avgtime = '30000',
       avgbaseline = T,
       avgscan = T,
       uvrange='0~100m',
       transform=T,
       freqframe='LSRK')


#----- cvel

os.system('rm -rf HD163296.C18OLine.contsub.vslr')

cvel(vis = linename+'.contsub',
     outputvis = linename+'.contsub.vlsr',
     field='0',
     spw='0',
     selectdata  =  F,
     mode  =  'channel',
     nchan = 150,
     start = 700,
     width = 1,
     interpolation  = 'linear',
     phasecenter =  0,
     restfreq  = '219.56036GHz',
     outframe  =  'lsrk',
     veltype  = 'radio',
    hanning = F,
    async = F)


listobs(linename+'.contsub.vlsr')


#----- clean

os.system('rm -rf HD163296_C18O*')

clean(vis = linename+'.contsub.vlsr',
      imagename = 'HD163296_C18O',
      field = 'HD163296',
      phasecenter = 0,
      cell = '0.05arcsec',
      imsize = [432,432],
      spw = '0',
      mode = 'channel',
      start = 1,
      width = 1,
      nchan = 150,
      restfreq='219.56036GHz',
      weighting = 'briggs',
      niter = 1500,
      interactive = T,
      imagermode = 'csclean',
      cyclefactor = 5)

viewer('HD163296_C18O.image')

#----- Export as fits
exportfits(imagename='HD163296_13CO.image',
           fitsimage='HD163296_13CO.image.fits')
exportfits(imagename='HD163296_C18O.image',
           fitsimage='HD163296_C18O.image.fits')
exportfits(imagename='HD163296_CO_2_1.image',
           fitsimage='HD163296_CO_2_1.image.fits')
exportfits(imagename='HD163296_Band6.CalibCont.sc2.image',
           fitsimage='HD163296_Band6.CalibCont.sc2.image.fits')
           

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

name='HD163296_CO_2_1.image'

immoments(imagename=name, 
  moments=[0],
  chans='', 
  box='', 
  axis='spectral',
  includepix=[0.01, 10000], 
  outfile=name+'.mom0')

immoments(imagename=name, 
  moments=[1],
  chans='', 
  box='', 
  axis='spectral',
  includepix=[0.3, 10000], 
  outfile=name+'.mom1')


#----- Export as fits

exportfits(imagename=name+'.mom0',
     fitsimage=name+'.mom0.fits')
exportfits(imagename=name+'.mom1',
     fitsimage=name+'.mom1.fits')


#----- Create Moment Maps for 13CO

name='HD163296_13CO.image'

immoments(imagename=name, 
  moments=[0],
  chans='', 
  box='', 
  axis='spectral',
  includepix=[0.01, 10000], 
  outfile=name+'.mom0')

immoments(imagename=name, 
  moments=[1],
  chans='', 
  box='', 
  axis='spectral',
  includepix=[0.3, 10000], 
  outfile=name+'.mom1')


#----- Export as fits

exportfits(imagename=name+'.mom0',
     fitsimage=name+'.mom0.fits')
exportfits(imagename=name+'.mom1',
     fitsimage=name+'.mom1.fits')

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