########################################################################
# ALMA Data reduction script for Comet Lemmon Band 6 data: Imaging
# CometLemmon_Band6_Imaging.py
# Tested using CASA Version 4.2.1 (r29048)
#
# Datasets used:
# uid___A002_X62f759_X2363.ms 
# - Observed on 11 May 2013 with 25 12-m antennas
# --- See Readme file for additional details ---
########################################################################
"""
See accompanying README file for details of the necessary input files
and comments on the data.

Prior to running this script you will need to obtain the Analysis Utilities package.
See http://casaguides.nrao.edu/index.php?title=Analysis_Utilities for details and 
for how to download Analysis Utilities.

"""

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

version = casalog.version()
print "You are using " + version
if (version < '4.2.1'):
    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."


#----------------------------------------------------------------------------------
#----- Make Comet Lemmon continuum image ------------------------------------------
#----------------------------------------------------------------------------------

os.system('rm -rf Lemmon_Cont_natural*' )
clean(vis = 'uid___A002_X62f759_X2363_12m.ms.split.cal',
  imagename = 'Lemmon_Cont_natural',
  field = '2', # Lemmon
  spw = '0:20~1850;2000~3800,1,2,3',
  mode = 'mfs',
  interactive = F,
  imsize = [350, 350],
  cell = '0.04arcsec',
  phasecenter = '00:19:06.3  +17.08.22.7',
  weighting = 'natural',
  robust = 0.5,niter=500,threshold = '0.14mJy')

#----------------------------------------------------------------------------------
#----- Clean the Comet Lemmon line data -------------------------------------------
#----------------------------------------------------------------------------------

os.system('rm -rf Lemmon_HCNLine_natural*' )
clean(vis = 'uid___A002_X62f759_X2363_12m.ms.split.cal',
  imagename = 'Lemmon_HCNLine_natural',
  field = '2', # Lemmon
  spw = '0',
  mode = 'channel',
  start=1800,
  nchan=230,
  interactive = F,
  mask='Lemmon_hcn_natural.mask',
  imsize = [350, 350],
  cell = '0.04arcsec',
  restfreq='265.8864GHz',
  phasecenter = '00:19:06.3  +17.08.22.7',
  weighting = 'natural',
  robust = 0.5,
  threshold = '18mJy',
  niter=1000)

#----------------------------------------------------------------------------------
#----- Make moment maps of the HCN emission -----------------------------------
#----------------------------------------------------------------------------------

myimage='Lemmon_HCNLine_natural.image'
chanstat=imstat(imagename=myimage,chans='4')
rms1= chanstat['rms'][0]
chanstat=imstat(imagename=myimage,chans='226')
rms2= chanstat['rms'][0]
rms=0.5*(rms1+rms2)
print 'rms in a channel = '+str(rms)

os.system('rm -rf Lemmon_HCNLine_natural.image.mom0')

immoments(imagename = 'Lemmon_HCNLine_natural.image',
         moments = [0],
         axis = 'spectral',chans = '102~144',
         includepix = [rms*3,100.],
         outfile = 'Lemmon_HCNLine_natural.image.mom0')

os.system('rm -rf Lemmon_HCNLine_natural.image.mom1')
immoments(imagename = 'Lemmon_HCNLine_natural.image',
         moments = [1],
         axis = 'spectral',chans = '102~144',
         includepix = [rms*5.5,100.],
         outfile = 'Lemmon_HCNLine_natural.image.mom1')


#----------------------------------------------------------------------------------
#----- Export data as fits --------------------------------------------------------
#----------------------------------------------------------------------------------

os.system('rm -rf Lemmon_HCNLine_natural*.fits')
exportfits(imagename='Lemmon_HCNLine_natural.image', fitsimage='Lemmon_HCNLine_natural.image.fits')
exportfits(imagename='Lemmon_HCNLine_natural.image.mom0', fitsimage='Lemmon_HCNLine_natural.image.mom0.fits')
exportfits(imagename='Lemmon_HCNLine_natural.image.mom1', fitsimage='Lemmon_HCNLine_naturals.image.mom1.fits')


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