#
# This script to explores CASA imaging of spectral line data sets
# using ALMA Science Verification data. So far we have looked at
# imaging continuum data. The main difference with line data is that
# we image each frequency separately. CLEAN includes several
# parameters that let you do this. Also, the continuum that we imaged
# before is still present in each data set and will need to be removed
# before we can study the line.
#

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# SET INPUTS AND OUTPUTS
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

#
# Pick the input data set and set the root name of the output
# data. These are read into variables that are then fed to parameters
# in the calls to CLEAN and continuum subtraction below.
#

# Path to the calibrated data (pick one of these)

data = "../../calib/NGC3256_Band3_CalibratedData/ngc3256_line_target.ms"
local = "ngc3256_line_target.ms"

#data = "../../calib/TWHYA_BAND7_CalibratedData/TWHydra_corrected.ms
#local = "TWHydra_corrected.ms"

#data = "../../calib/TWHYA_BAND6_CalibratedData/TWHydra_corrected.ms"
#local = "TWHydra_corrected.ms"

#data = "../../calib/TWHYA_BAND3_CalibratedData/TWHydra_corrected.ms"
#local = "TWHydra_corrected.ms"

#
# List the basic parameters of the input data set. 
#

listobs(vis=data)

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# COPY THE DATA
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

#
# We're going to be manipulating the data. To keep things clean, use
# split to make a copy that contains ONLY the line that we care
# about. Note that you need to know something about your data before
# you would know the mapping of SPW to line (e.g., you might know that
# SPW 0 contains the CO line because you set up the NGC 3256
# observations to do that). Included below are "SPLIT" calls to
# separate out lines for each data set.
#

# NGC 3256 CO
local = 'ngc3256_co.ms'
os.system('rm -rf '+local) # needed for repeated running
split(vis=data,
      outputvis=local,
      spw='0',
      datacolumn='data')

# NGC 3256 CN
#local = 'ngc3256_cn.ms'
#os.system('rm -rf '+local) # needed for repeated running
#split(vis=data,
#      outputvis=local,
#      spw='1',
#      datacolumn='data')

# TW HYDRA BAND 3
#local = 'tw_hydra_band3_hcoplus.ms'
#os.system('rm -rf '+local) # needed for repeated running
#split(vis=data,
#      outputvis=local,
#      spw='1',
#      datacolumn='data')

# TW HYDRA BAND 6 CO AND DCN3-2
#local = 'tw_hydra_band6_co.ms'
#os.system('rm -rf '+local) # needed for repeated running
#split(vis=data,
#      outputvis=local,
#      spw='0',
#      datacolumn='data')

#local = 'tw_hydra_band6_dcn32.ms'
#os.system('rm -rf '+local) # needed for repeated running
#split(vis=data,
#      outputvis=local,
#      spw='1',
#      datacolumn='data')

# TW HYDRA BAND 7 CO AND HCO+
#local = 'tw_hydra_band7_co.ms'
#os.system('rm -rf '+local) # needed for repeated running
#split(vis=data,
#      outputvis=local,
#      spw='2',
#      datacolumn='data')

#local = 'tw_hydra_band7_hcoplus.ms'
#os.system('rm -rf '+local) # needed for repeated running
#split(vis=data,
#      outputvis=local,
#      spw='0',
#      datacolumn='data')

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# UV CONTINUUM SUBTRACTION
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

#
# Find the line free channels by making an average spectrum that
# combines all times and baselines (and so maximizes
# signal-to-noise). This won't always work but it will work for our
# test data sets. Note the channel ranges where it looks like there is
# no line. You will feed these in to the uv continuum subtraction in
# the next step.  As an example, see the NGC 3256 CN lines here:
# http://casaguides.nrao.edu/index.php?title=File:Amp_vs_chan_spw1.png

plotms(vis = local,
       spw='0',
       xaxis='channel',
       yaxis='amp',
       avgtime='1e9',
       avgscan=True,
       avgbaseline=True,
       ydatacolumn='data')

# SUGGESTED CALLS TO UV CONTINUUM SUBTRACTION

# Try to figure the line-free channels out on your own from the plots
# above but these are the suggested calls to uvcontsub2 from the CASA
# guides and associated scripts. Notice what happens: these calls
# produce a NEW measurement set with the extension ".contsub". We will
# image that measurement set.
# Be sure to use the correct call for the line you specified in the
# "split" call above. 

# ... for NGC 3256 CO
uvcontsub(vis = local,
          fitspw='0:20~53;71~120',
          solint ='int',
          fitorder = 1,
          combine='')

# ... for NGC 3256 CN
#uvcontsub(vis = local,
#          fitspw='0:70~120',
#          solint ='int',
#          fitorder = 1,
#          combine='')

# ... for TW Hydra Band 3
#uvcontsub(vis=local,
#          fitspw='0:20~400,0:2300~3800',
#          solint='int',
#          fitorder=1,
#          combine='')

# ... for TW Hydra Band 6 CO
#uvcontsub(vis=local,
#          fitspw='0:20~2000,0:2400~3800',
#          solint='int',
#          fitorder=1,
#          combine='')

# ... for TW Hydra Band 6 DCN3-2
#uvcontsub(vis=local,
#          fitspw='0:40~800;1200~3800',
#          solint='int',
#          fitorder=1,
#          combine='')

# ... for TW Hydra Band 7 CO
#uvcontsub(vis=local,
#          fitspw='0:20~2000,0:2400~3800',
#          solint='int',
#          fitorder=1,
#          combine='')

# ... for TW Hydra Band 7 HCOplus
#uvcontsub(vis=local,
#          fitspw='0:20~1500,0:1900~3800',
#          solint='int',
#          fitorder=1,
#          combine='')

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# IMAGE THE CONTINUUM SUBTRACTED DATA
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

# Set the parameters that CLEAN will use to grid and deconvolve the
# u-v data into your final image. The main difference between this and
# the continuum case is that we will specify the third (spectral) axis
# of the data. We will see several ways to do this.

default('clean')

# image the continuum-subtracted data

vis=local+'.contsub'

# specify the output (NB we cheat on the name here using a minor
# python trick to strip the trailing '.ms')

imagename=str.rstrip(local, '.ms')
print imagename

# we want to do this interactively
interactive=True

# ... STOPPING CONDITIONS
thresh = '1mJy'
niter = 10000

# ... CELL SIZE (WANT ~3-5 CELLS PER SYNTHESIZED BEAM)
# Good for Bands 6 and 7
# cell = '0.3arcsec'
# cell = '0.5arcsec'
# Good for Band 3
cell = '1.0arcsec'

# ... IMAGE SIZE (TRY TO IMAGE OUT TO ~ PRIMARY BEAM BY DEFAULT)
# Good for Bands 6 and 7
# imsize = 100
# Good for Band 3
imsize = 300

imagermode = 'csclean'
psfmode='hogbom'

# ... U-V WEIGHTING
# weighting = 'natural'
weighting = 'briggs'
robust = 0.5

# ----------------------
# THE LINE-SPECIFIC PART
# ----------------------

# You can specify your spectral range and gridding in a few ways. The
# CASA guides adopt different approaches, see them for more details.

# For each of these cases you will want to start by specifying the
# rest frequency (again, remember to match the line you specified in
# the "split" call):

# NGC 3256 CO 1-0
restreq = '115.271201800GHz'

# NGC 3256 CN N=1-0 J=3/2-1/2
#restfreq = '113.48812GHz'

# NGC 3256 CN N=1-0 J=1/2-1/2
#restfreq = '113.17049GHz'

# TW Hydra HCO+ 1-0
#restfreq = '89.188526GHz'

# TW Hydra CO 2-1
#restfreq = '230.538GHz'

# TW Hydra DCN
#restfreq = '217.2385GHz'

# TW Hydra CO 3-2
#restfreq='345.79599GHz'

# TW Hydra HCO+ 4-3
#restfreq='356.7342GHz'

#
# APPROACH 1: mode 'channel'. Pick a channel range via spw
# selection. You can refine this with a specified start channel,
# width, and number of channels to image.  As an example, you can look
# at the NGC 3256 CASA guide:
# http://casaguides.nrao.edu/index.php?title=NGC3256_Band3_-_Imaging
# with an image of the interactive clean viewer for the CO line here:
# http://casaguides.nrao.edu/index.php?title=File:Interactive_clean_channel.png

# ... Here for NGC 3256 CO
mode='channel'
start=''
spw='0:38~87'
width=''
nchan=50

# ... Here for NGC 3256 CN "low" line
#mode='channel'
#start=''
#spw = '0:29~54'
#width=''
#nchan=26

# ... Here for NGC 3256 CN "high" line
#mode='channel'
#start=''
#spw = '0:50~76'
#width=''
#nchan=27

#
# APPROACH 2: mode 'velocity'. Pick a starting velocity, step size,
# and number of channel maps to image. See examples in the TW Hydra
# CASA guide. CASA handles the appropriate regridding.  As an example,
# a screen shot of the clean viewer for the TW Hydra Band 7 data is
# here: http://casaguides.nrao.edu/index.php?title=File:Viewer_55.png

# ... for TW Hydra Band 7 (either HCO+ or CO 3-2)
#mode='velocity'
#start='-4km/s'
#width='0.12km/s'
#nchan=118
#outframe='LSRK'

# ... for TW Hydra Band 6
#mode='velocity'
#start='-4km/s'
#width='0.2km/s'
#nchan=70
#outframe='LSRK'

# ... for TW Hydra Band 3
#mode='velocity'
#start='-1km/s'
#width='0.25km/s'
#nchan=36
#outframe='LSRK'

# You also have the option to specify mode="Frequency"

# ----------------------
# FIRE OFF CLEAN
# ----------------------

# erase previous images. Skip this step if you want to continue
# cleaning or skip the mask step specifically if you want to save your
# clean mask for a new CLEAN call. Note that the mask *must* match
# your image cube in shape, so if you change the line parameter you
# must remake the mask.
os.system('rm -rf '+imagename+'.image')
os.system('rm -rf '+imagename+'.model')
os.system('rm -rf '+imagename+'.psf')
os.system('rm -rf '+imagename+'.mask')
os.system('rm -rf '+imagename+'.residual')
os.system('rm -rf '+imagename+'.flux')

# review inputs
inp

# execute CLEAN
go

# Remember we are still performing an interactive CLEAN. Be sure you make
# a CLEAN box (or boxes) around the emission in the cube before continuing
# to clean. 

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# VIEW
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

#
# Display the resulting image. Notice that this is a cube. Use the
# movie buttons and the spectral window browser to look along the
# third dimenstion.
#

viewer(imagename+'.image')

# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%
# MAKE MOMENT MAPS
# =%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%=%

#
# You can collapse the cube into an image using the "immoments"
# task. This allows a number of operations along the axis of your
# choice (traditionally the spectral one). The "includepix" option
# specifies the intensity range of pixels to include. Optionally, a
# more sophisticated mask can be supplied as an argument (as either a
# box, region file, or image).
#

# Make a moment-0 (integrated intensity) image ...
pixrange = [0.25,1e3] # for NGC3256
# pixrange = [0.3,1e3] # for TW Hydra Band 7
os.system('rm -rf '+imagename+'.mom0')
immoments(imagename=imagename+'.image',
          moments=0,
          includepix=pixrange,
          outfile=imagename+'.mom0')
viewer(imagename+'.mom0')

# and a moment-1 (intensity-weighted velocity) image. 
pixrange = [0.25,1e3] # for NGC3256
# pixrange = [0.3,1e3] # for TW Hydra Band 7
os.system('rm -rf '+imagename+'.mom1')
immoments(imagename=imagename+'.image',
          moments=1,
          includepix=pixrange,
          outfile=imagename+'.mom1')
viewer(imagename+'.mom1')