#######################################################################
###### tuCASA – ASIAA ALMA Imaging Workshop 2013 (January 22–24) ######
######     hands-on session I & II: Data Inspection & Imaging    ######
#######################################################################

# To reduce the data size, the data of CH=1500-2000 at spw=17 were extracted.
# This frequency range covers the HCO+(4-3) line.

# Link to the casaguides for TW Hydra:
# http://casaguides.nrao.edu/index.php?title=TWHydraBand7_Calibration



#########################
#### Data Inspection ####
#########################

# Data are already calibrated.

# observation sequence (see plot plot_observation_sequence.png):
# field 1 (pink):   Titan     - flux calibrator
# field 0 (black):  3C279     - bandpass calibrator
# field 3 (green):  J1147-382 - secondary phase calibrator
# field 4 (brown):  J1037-295 - primary phase calibrator
# field 2 (orange): TWHya     - science target

# The TW Hydra data are already split out and are used
# for the following data inspection and imaging.
# The file name of the measurement set is:

msname='TWHydra_calibrated'

# antenna positions
plotants(vis=msname+'.ms',figfile='plot_antenna_positions.png')

# list data summary
listobs(vis=msname+'.ms',verbose=T,listfile=msname+'.listobs')


# see structure (in this case decreasing amplitude with uv-distance)
# => the source is resolved
plotms(vis=msname+'.ms',spw='',xaxis='uvdist',yaxis='amp',
       field='',avgchannel='501', antenna='*&*')

# check where the line is
# => the line is between channels 180 and 220
plotms(vis=msname+'.ms',spw='',xaxis='channel',yaxis='amp',
       field='',avgtime='1e9',avgscan=T, antenna='*&*')

# split out the line and average over 2 channels
os.system('rm -rf TWHydra_line.ms')
split(vis=msname+'.ms',outputvis='TWHydra_line.ms',field='0',
      spw='0:197~220',width=2,datacolumn='data')



#########################
####     Imaging     ####
#########################

# 1: clean and map continuum
# 2: subtract continuum
# 3: clean and map line

msname='TWHydra_calibrated'


# 1: CLEAN AND MAP THE CONTINUUM
# ------------------------------

# clean the continuum interactively:
# - draw a circle around the central emission
# - click inside the circle to select this as cleaning box (it turns from green to white)
# - click the green circle arrow to start the cleaning
#   => only a few iterations are enough
# - stop with X button

os.system('rm -rf TWHydra_continuum.*')
clean(vis=msname+'.ms',
      imagename='TWHydra_continuum',
      imagermode='csclean',
      spw='0:0~180;220~500',
      field='0',
      mode='mfs',
      imsize=100,
      cell=['0.3arcsec'],
      weighting='briggs',
      robust=0.5, 
      mask=[],
      interactive=T,
      threshold='1mJy',
      niter=10000)

# output files:
# TWHydra_continuum.flux            # the effective primary beam response
# TWHydra_continuum.image           # the final restored image
# TWHydra_continuum.model           # the model image
# TWHydra_continuum.psf             # the synthesized (dirty) beam 
# TWHydra_continuum.residual        # the residual image

# view image and get the rms noise level
imview('TWHydra_continuum.image')
imstat(imagename='TWHydra_continuum.image',box='0,0,40,40')


# 2: SUBSTRACT CONTINUUM USING UVCONTSUB
# --------------------------------------

# fit first order baseline and subtract continuum
# the output file name is the input file name plus '.contsub' added
uvcontsub(vis=msname+'.ms', fitorder=1, fitspw='0:0~160;240~500')

# check the continuum subtracted uv data
plotms(vis=msname+'.ms.contsub',spw='',xaxis='channel',yaxis='amp',
       field='',avgtime='1e9',avgscan=T, antenna='*&*')


# 3: CLEAN AND MAP LINE
# ---------------------

# clean the line interactively
# - use the slider to find the channels with emission
# - select 'All Channels' before drawing a clean box
# - only a few iterations are enough

os.system('rm -rf TWHydra_line.*')
clean(vis=msname+'.ms.contsub',
      imagename='TWHydra_line',
      imagermode='csclean',
      field='0',
      spw='',
      imsize=100,
      cell=['0.3arcsec'],
      mode='velocity',
      start='-4km/s',
      nchan=118,
      width='0.12km/s',
      restfreq='356.7342GHz',
      outframe='LSRK',
      weighting='briggs',
      robust=0.5,
      mask=[],
      interactive=T,
      threshold='1mJy',
      niter=100000)

# view image and get the rms noise level
# => the rms noise level should be at around 45 mJy/beam
imview('TWHydra_line.image')
imstat(imagename='TWHydra_line.image',chans='10~40')


# NOTES:
# ---------------------
# - if you want to see the dirty image only, set 'niter=0'
# - if you want to do the clean non-interactively,
#   set 'interactive=T' and set a reasonable number for
#   'threshold' or 'niter'
# - to get help, type 'help clean' (type 'q' to exit the help)
# - to see all available parameters for clean type:
#   > default clean
#   > inp
#   after setting all parameters, one can run clean by typing:
#   > go

#######################################################################