########################################################################
# Data reduction script for Orion Spectral Sweep in Band 6:
# - Imaging script: "Orion_SpectralSweep_Imaging.py" -
# Tested in CASA Version 3.3.0 (r16856)
########################################################################

basename='Orion_all'

#split out and clean phasecal data as a test, use spw 4 only (the one containing CO)
os.system('rm -rf phasecal_spw4')
split(vis=basename+'.ms',outputvis='phasecal_spw4',
      field='0',spw='4',width=1,datacolumn='data')

os.system('rm -rf phasecal.*')
clean(vis='phasecal_spw4',imagename='phasecal',
       field='',spw='',mode='mfs',calready=T,
       imagermode='csclean',imsize=100,cell=['0.4arcsec'],
       weighting='natural',psfmode='hogbom',mask='',
       interactive=T,
       threshold='1.5mJy',niter=10000)

#split out and image the continuum at ~231 GHz
os.system('rm -rf Orion.cont231')
split(vis=basename+'.ms',outputvis='Orion.cont231',
       field='1',spw='4:1850~1880',width=1,datacolumn='data')

os.system('rm -rf Orion.cont.*')
clean(vis='Orion.cont231',imagename='Orion.cont',
       field='',spw='',mode='mfs',calready=T,
       imagermode='csclean',imsize=100,cell=['0.4arcsec'],
       weighting='natural',psfmode='hogbom',mask='',
       interactive=T,threshold='1.5mJy',niter=10000)

# split out the isolated and strong CH3OH (Methanol) J=8(-1,8)-7(0,7) line at 229.76 GHz.
# split out line data, from 229.74-229.77 GHz

os.system('rm -rf Orion.methanol')
split(vis=basename+'.ms',outputvis='Orion.methanol',
      field='1',spw='19:1150~1210',
      width=1,datacolumn='data')

# get rms from one of the line-free channels: rms=0.023 Jy.

os.system('rm -rf Orion.methanol.dirty')
clean(vis='Orion.methanol',imagename='Orion.methanol.dirty',
       field='',spw='',mode='channel',
       outframe='LSRK',restfreq='229.759GHz',calready=T,
       imagermode='csclean',imsize=100,cell=['0.4arcsec'],
       weighting='natural',psfmode='hogbom',mask='',
       interactive=T,
       threshold='1.5mJy',niter=0)

#In what follows try to clean to ~3sigma, so use a threshold of 0.07 Jy.
# subtract the continuum,split out the central ~40 channels of the methanol line,
# cleaning channel-by-channel with making a new clean box for every channel.
# Such a channel-by-channel mask is provided for convenience.

uvcontsub2(vis='Orion.methanol',fitorder=1,fitspw='0:0~5,0:52~60')

os.system('rm -rf Orion.methanol.line.contsub')
split(vis='Orion.methanol.contsub',outputvis="Orion.methanol.line.contsub",
       spw='0:10~50',datacolumn='data')

# Clean using the provided channel-by-channel mask
os.system('rm -rf Orion.methanol.cbc.contsub.*')
clean(vis='Orion.methanol.line.contsub',imagename='Orion.methanol.cbc.contsub',
       field='',spw='',mode='channel',
       outframe='LSRK',restfreq='229.759GHz',calready=T,
       imagermode='csclean',imsize=100,cell=['0.4arcsec'],
       weighting='natural',psfmode='hogbom',mask='Orion.methanol.cbc.mask',
       interactive=F,threshold='70mJy',niter=100000)

# make a moment map from this; use only the central channels
# moment 0: total intensity
os.system('rm -rf Orion.methanol.mom0')
immoments(imagename='Orion.methanol.cbc.contsub.image',
          moments=[0],outfile='Orion.methanol.mom0',chans='5~34')

# moment 1: velocity field
# set flux threshold.
os.system('rm -rf Orion.methanol.mom1')
immoments(imagename='Orion.methanol.cbc.contsub.image',
          moments=[1],outfile='Orion.methanol.mom1',
          chans='5~34',includepix=[1.6,100])

#----- Export as fits
os.system('rm -rf Orion.cont.image.fits')
exportfits(imagename='Orion.cont.image',
     fitsimage='Orion.cont.image.fits')

os.system('rm -rf Orion.methanol.cbc.image.fits')
exportfits(imagename='Orion.methanol.cbc.contsub.image',
     fitsimage='Orion.methanol.cbc.contsub.image.fits')

os.system('rm -rf Orion.methanol.mom0.fits')
exportfits(imagename='Orion.methanol.mom0',
     fitsimage='Orion.methanol.mom0.fits')

os.system('rm -rf Orion.methanol.mom1.fits')
exportfits(imagename='Orion.methanol.mom1',
     fitsimage='Orion.methanol.mom1.fits')

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