#!/bin/csh -f
# handling standard 2GHz uniform data
#
# Responsible: Edmund Hodges-Kluck (UMd) & Jun-Hui Zhao (SAO)
# ehk & jhz 2007-12-05 produce the template scripts for handling std SMA data  
##############################################################
# Edmund Hodges-Kluck (UMd) & Jun-Hui Zhao (CfA) 12/5/07
# for the reduction of SMA data taken recently (2007 or later).
# The SMA data, used in this example, were taken in a spectral
# configuration with a hybrid resolution. The observations had 
# multiple pointing for mosaic image of molecular lines. This
# E2E script consists of nine major steps with some detailed 
# comments
# 1. READ THE RAW SMA DATA INTO MIRIAD
# 2. INSPECT THE RAW UV DATA
# 3. FRINGE AMPLITUDE CORRECTION
# 4. BANDPASS CALIBRATION
# 5. ANTENNA GAIN CALIBRATION
# 6. BOOTSTRAP THE FLUX DENSITY SCALE
# 7. CONTINUUM IMAGE
# 8. LINE/CONTINUUM SEPARATION
# 9. CONSTRUCT LINE IMAGE CUBE
##############################################################
### INTRODUCTION #############################################
##############################################################
# This is a full reduction in one line (HCO+ 3-2) looking at
# a data set with data in both HCO+ and HCN 3-2 in the lower
# sideband.  To run the script, sequentially step through and
# uncomment the various 'continue' markers and run through
# until the script exits, then comment the marker out and 
# uncomment the marker for the next step.  Refer to the SMA
# MIRIAD manual for additional information.  Note that we
# find a non-detection on the source, but get nice images of
# Sgr A*.
##############################################################

  set source=S1
  goto continue

### LOAD THE DATA ############################################
# In this section we take the data from its original MIR
# form and load it into MIRIAD.
##############################################################

#continue:

  set darea=/2006/science/mir_data/071004_00:26:54/
  set darea2=/2006/science/mir_data/071004_04:41:13/
  echo "1. READ THE RAW SMA DATA INTO MIRIAD"
  \rm -r 1E1740_1_rx1.lsb 1E1740_2_rx1.lsb 
  smalod in=$darea out=1E1740_1 rxif=1 \
         sideband=0 readant=8 nscans=1
  smalod in=$darea2 out=1E1740_2 rxif=1 \ 
         sideband=0 readant=8 nscans=1

# We had 2 files for the observation, so load in the data
# twice and merge.

  \rm -r 1E1740_rx1.lsb
  uvaver vis=1E1740_1_rx1.lsb,1E1740_2_rx1.lsb out=1E1740_rx1.lsb

exit

### INSPECT RAW DATA ########################################
# In this section we look first at the sources in the
# observation (uvindex) and then the amplitude vs.
# time and the system temperature vs. time
#############################################################

#continue:

  echo "2. INSPECT THE RAW UV DATA"
  echo "List the observing scans in the uv data"

  uvindex vis=1E1740_rx1.lsb
  echo "Press  to continue"
  $<

  smauvplt vis=1E1740_rx1.lsb \
   	axis=time,ampl device=/xs nxy=1,1 \
           yrange=0,300
  echo "Press  to continue"
  $<

  smavarplt vis=1E1740_rx1.lsb device=/xs \
            xaxis=time \
            yaxis=systemp  nxy=2,4
exit

### FRINGE AMPLITUDE CORRECTION #############################
# The changes in the system temperature over the duration
# of the observation essentially track changes in the
# atmosphere.  The system temperature measurements are
# obtained by comparing the receiver output power from
# an ambient temperature load with the sky power, so we
# can calibrate the amplitude gains.
#############################################################

#continue:

  echo "3. FRINGE AMPLITUDE CORRECTION"
  \rm -r 1E1740_rx1.lsb.tsys
  smafix vis=1E1740_rx1.lsb out=1E1740_rx1.lsb.tsys \
   	device=/xs xaxis=time \
         yaxis=systemp nxy=2,4   \
         options=tsyscorr
  echo "Inspect vis data and flag bad data"

# We could filter out various times while we're doing the
# flagging here.  We can also do this later when looking
# at the bandpass calibrator (as is the case in this script)

# Jupiter was just there for pointing and makes it harder to do
# interactive flagging.
  uvflag vis=1E1740_rx1.lsb.tsys select='source(jupiter)' flagval=flag

# Antenna 4 was apparently not on, flag it just in case
  uvflag vis=1E1740_rx1.lsb.tsys select='ant(4)' \
 	flagval=flag

# Interactively flag bad data
  smablflag vis=1E1740_rx1.lsb.tsys axis=time,ampl device=/xs

exit

#continue:
# Look at the bandpass calibrator
  echo "Inspect spectral data"
  echo "for uranus"
  smauvspec vis=1E1740_rx1.lsb.tsys \
 	select='source(uranus)' \
         interval=1000  \
         axis=freq,both device=/xs nxy=1,5 hann=5
# SgrA* too
#continue: 
  smauvspec vis=1E1740_rx1.lsb.tsys \
         select='source(sgra*,Sg*)' \
         interval=1000  \
         axis=freq,both device=/xs nxy=1,5 hann=5
exit

### BANDPASS CALIBRATION ####################################
# Bandpass calibration gives the frequency solution, so
# we use something bright -- in the case of the SMA, a
# planet typically.  In this case, Uranus.
#############################################################

#continue:

  echo "4. BANDPASS CALIBRATION"

# smamfcal makes calibration corrections (antenna gains, delay
# terms and passband shapes) to a multi-frequency observation.
# The gains are calculated periodically using the interval 
# specified.  If phase stability is poor or you are using a
# planet, set weight=2 (see manual).  weight=2 is default for 
# the SMA.

  smamfcal vis=1E1740_rx1.lsb.tsys select='source(uranus)' edge=2.2  \
          weight=2 refant=3 interval=600
  echo "Check up the bandpass solutions"
  smagpplt vis=1E1740_rx1.lsb.tsys device=/xs \
         yaxis=amp,phase \
         options=bandpass  polyfit=5 \
         nxy=2,4
#exit
#continue:

  echo 'apply the bpass solutions to the data'
  \rm -r 1E1740_rx1.lsb.tsys.bp
  uvaver vis=1E1740_rx1.lsb.tsys out=1E1740_rx1.lsb.tsys.bp
  smauvspec vis=1E1740_rx1.lsb.tsys.bp \
         select='source(sgr*,Sgr*)' \
         interval=1000   \
         axis=freq,both device=/xs nxy=2,2

exit
#continue:

# Interactively flag again to get anything still messing up
# the solution.
  smablflag vis=1E1740_rx1.lsb.tsys.bp \
         axis=time,ampl device=/xs
exit
#continue:

  uvflag vis=1E1740_rx1.lsb.tsys.bp edge=15,15 \
         flagval=flag  select='window(1,2,3,4)'

# We have to vary the windows we're looking to clean up.
# For this observation,
#  select='window(20,21,22,23,24)' selects the HCN (3-2)
#  select='window(1,2,3,4)' selects the HCO+ (use this for now,
#         it's brighter!)

  smauvspec vis=1E1740_rx1.lsb.tsys.bp \
         select='source(sgr*,Sgr*),ant(1)(6)' \
         interval=1000  \
         axis=freq,both device=/xs nxy=1,1 
exit

### ANTENNA GAIN CALIBRATION ################################
# The gain calibration sources for this observation were
# 1733-130, SgrAstar.  Here we apply the gain calibration to 
# the data.
#############################################################

#continue:

  echo "5. ANTENNA GAIN CALIBRATION"

# Filter out bad times here -- we decide if/what to filter
# based on the smauvplt below, so may have to run this step
# again with the line below uncommented and changed as
# appropriate.

# NOTE: May need to exclude 1.75h-3h based on the 
# system temperature/phase problems.
#
# uvflag vis=1E1740_rx1.lsb.tsys.bp \ 
#        select='time(1:45,3:00)' flagval=flag
#
# If you are unhappy with the result of the flagging, you can
# return to this step and run:
#
# uvflag vis=1E1740_rx1.lsb.tsys.bp \ 
# 	   select='time(1:45,3:00)' flagval=unflag
#
# However, doing so may unflag other data you didn't want to unflag,
# so you might have to run smablflag again to interactively flag
# the data again.  If you are worried, start over and don't flag
# based on time.  You might also create a new file to do the time
# filtering on.

  smablflag vis=1E1740_rx1.lsb.tsys.bp  device=/xs axis=real,imag \
  	select='source(1733*,Sgr*)' options=nobase

  echo "Press  to continue."
  $<
  smauvplt  vis=1E1740_rx1.lsb.tsys.bp \
 	select='source(1733*,Sgr*)' \
 	device=/xs options=nocal axis=time,amp nxy=1,3

  echo "Split the gain calibrators from the multiple source file"

#continue:
# The number after these files will differ based on the dataset.
# Therefore, this step may need to be run twice.

  \rm -r 1733-130.267745 SgrAstar.267745 
  uvsplit vis=1E1740_rx1.lsb.tsys.bp \
 	select='source(1733-130,SgrAstar)' options=nowindow

  echo "rename the calibrator files"
  \rm -r 1733-130.tsys.bp sgra.tsys.bp
  mv 1733-130.267745 1733-130.tsys.bp 
  mv SgrAStar.267745 sgra.tsys.bp 

exit
#continue:

  echo " "
  echo "Solve for antenna gains"
  echo " "
# Can we use selfcal?  Let's hope so.  They're both strong point
# sources in this case.
# Note that here we do phase only so that we can solve for
# phase only.
  selfcal vis=1733-130.tsys.bp interval=1 options=phas refant=3
  selfcal vis=sgra.tsys.bp interval=1 options=phas refant=3
  echo "Plot the gain solutions"
  smagpplt vis=1733-130.tsys.bp device=/xs \
        yaxis=phase log=0102+584gain.log \
        options=gains  polyfit=6 \
        nxy=2,4
  echo "Press  to continue."
  $<
  smagpplt vis=sgra.tsys.bp device=/xs \
        yaxis=phas \
        options=gains  polyfit=6 \
        nxy=2,4

exit
#continue:

  echo "Merge the two gain tables"
  smagpplt vis=1733-130.tsys.bp,sgra.tsys.bp device=/xs \
        yaxis=amp,pha  \
        options=gains,merge polyfit=8 \
        nxy=2,4
#the merged gain table is place in the second file, namely,
#sgra.tsys.bp
  smagpplt vis=sgra.tsys.bp device=/xs \
        yaxis=phas \
        options=gains  polyfit=6 \
        nxy=2,4
# Note that the default interpolation interval equals to
# the solution interval. One may make a longer interval
# in order to apply (interpolate) the selfcal solutions to
# the target sources.
# Modify the SgrA* header file so that the interval value
# is changed to 0.2day, then copy the gain table back to the
# main file with all sources.

  puthd in=sgra.tsys.bp/interval value=0.2
  gpcopy vis=sgra.tsys.bp out=1E1740_rx1.lsb.tsys.bp \
     mode=copy options=nopass

# See what the results are.

  smauvplt vis=1E1740_rx1.lsb.tsys.bp axis=time,amp device=/xs nxy=1,4
  smauvspec vis=1E1740_rx1.lsb.tsys.bp axis=freq,both \
         interval=1000000 options=nobase,avall \
         select='source(sgra*)' device=/xs nxy=1,1 hann=15

  echo 'apply the gain to the data'
  \rm -r 1E1740_rx1.lsb.tsys.bp.gn
  uvaver vis=1E1740_rx1.lsb.tsys.bp out=1E1740_rx1.lsb.tsys.bp.gn

exit

### BOOTSTRAP THE FLUX DENSITY SCALE ########################
# Find the scaling factor by looking at the bandpass
# calibrator, then apply the gain to the data
#############################################################

#continue:
 
  echo "6. BOOTSTRAP THE FLUX DENSITY SCALE"

# This will calculate the scaling factor based on your
# bandpass calibrator (the SMA function looks specifically
# for a planet) and then applies it to the data.  To
# calculate but not apply, use the option "noapply".

  smaflux vis=1E1740_rx1.lsb.tsys.bp.gn \
	select='source(uranus)' mirhome=$MIR  

# Specific to this dataset:
# Found solar system objects: uranus
# Frequency(GHz)   Tb(Kelvin)
#   267.797           93.61
# Scaling the data by 0.838

  echo "Apply the scaling factor gains to the data"
  \rm -r 1E1740_rx1.lsb.tsys.bp.gn.av
  uvaver vis=1E1740_rx1.lsb.tsys.bp.gn out=1E1740_rx1.lsb.tsys.bp.gn.av

# Check for the spectrum of Sgr A*

  smauvspec vis=1E1740_rx1.lsb.tsys.bp.gn.av \
         select='source(sgra*),window(1,2,3,4)' \
         interval=1000 hann=5  \
         axis=freq,both device=/xs nxy=1,1 options=nobase,avall
exit

#############################################################
# At this point, the data has had all the calibrations applied 
# to it. We now look at various sources within the data set.
#############################################################

# Change this to look at different sources, but remember to change
# manually in other places too as necessary

#continue:

  set source=SgrAStar

### CONTINUUM IMAGE #########################################
# We make an image (via inversion) of the source in the 
# continuum bands -- in this case we put lines in the 1,2,3,4 
# and 20,21,22,23,24 channels and had the rest as wide channels 
# with no lines.  Thus, we look at the source in all of these 
# wide channels.
#############################################################

  echo "7. CONTINUUM IMAGE"
  \rm -r $source.con.map $source.con.beam $source.con.icmp
  invert vis=1E1740_rx1.lsb.tsys.bp.gn.av    \
        map=$source.con.map beam=$source.con.beam \
        imsize=512 cell=0.6 sup=0 \
        select='window(5,6,7,8,9,10,11,12,13,14,15,16,17,18,19),source(SgrA*)' \
         line=chan,1,1,480 
  clean map=$source.con.map beam=$source.con.beam \
        out=$source.con.icmp gain=0.08 \
        cutoff=0 niters=3500 \
        region=arcsec,'box(-25,-25,25,25)'
  \rm -r $source.con.icln
  restor model=$source.con.icmp beam=$source.con.beam \
        map= $source.con.map  \
        out=$source.con.icln
  smacgdisp in=$source.con.icln type=pixel \
        region=arcsec,'boxes(-25,-25,25,25)' \
        xybin=1,1 device=/xs nxy=1,1 \ 
        options=full,beambr,wedge,trlab,3val \
        labtyp=arcsec,arcsec range=0,0,lin,2 cols1=7 csize=1,1,1
  imstat in=$source.con.icln region='box(20,20,80,80)'

exit

### LINE/CONTINUUM SEPARATION ################################
# In this section we subtract the continuum from the lines to 
# enhance contrasts in what is left. 
##############################################################

continue:

  set source=S1

  echo "8. LINE/CONTINUUM SEPARATION"

# In this case, our rest frequency was set to 267.55762 (HCO+ 3-2)

  \rm -r $source.hco+.uv

# You must set source here -- $source won't work

  uvaver vis=1E1740_rx1.lsb.tsys.bp.gn.av \
   select='source(S1_CR),window(1,2,3,4)'\
   out=$source.hco+.uv 

# Plot of amplitude vs. LSR velocity and phase vs. velocity
  smauvspec vis=$source.hco+.uv \
        interval=10000  options=nobase,avall \
        axis=vel,both device=/xs nxy=1,1 hann=15
  echo "Convert sky-frequency to velocity "
  echo "resolution of 1 km/s"
  $<

# The line=velocity argument selects the range in terms of
# velocities in the data; in this case we take essentially
# the full range.  The values used are determined by
# looking at the smauvspec plot from above, so you may run
# through this section a few times.

  \rm -r $source.hco+vel.uv
  uvaver vis=$source.hco+.uv line=velocity,390,-292,1 \
    out=$source.hco+vel.uv

# Plot of Amplitude vs. channels and phase vs. channels
  smauvspec vis=$source.hco+vel.uv \
         interval=10000  options=nobase,avall \
         axis=chan,both device=/xs nxy=1,1 hann=5

  echo "Take out the continuum using linear fitting to" 
  echo "the line free channels"
  \rm -r $source.hco+vel.uv.lin 

# uvlin performs the continuum subtraction by fitting a polynomial
# to the real and imaginary parts of the line-free channels (here
# specified)

  uvlin vis=$source.hco+vel.uv chans=1,100,290,390 \
      out=$source.hco+vel.uv.lin\
	order=1 mode=line

# chans=1,100,290,390 selects channels 1-100 and 290-390 for the
# continuum channels we use for subtraction -- from those we use
# linear interpolation to do the subtraction
# SgrA*: chans=5,80,150,220 (because there are different
# molecular features towards this source at various velocities than
# our target)
# 1E1740.7-2942: chans=1,100,290,390

  echo "Press  to continue."
  $<
 
  echo "Plot the continuum-free spectrum"
  smauvspec vis=$source.hco+vel.uv.lin options=nobase,avall  \
         interval=1000 \
         axis=vel,real device=/xs nxy=1,1 hann=15
  echo "Press  to continue."
  $<

### CONSTRUCT LINE IMAGE CUBE ################################
# Invert the visibilities in the line of the source and 
# deconvolve. Apply the model and look at the images.
##############################################################

  echo "9. CONSTRUCT LINE IMAGE CUBE" 
  \rm -r $source.hco+.map $source.hco+.beam $source.hco+.icmp

# A note regarding inversion: see the manual page.  The
# sidelobe suppression area (sup) is given in arcseconds, and
# gives the area around a source where inversion tries to
# suppress the sidelobes.  The default shape is square, the
# default behavior is suppression across the field.
# The visibility weighting robustness parameter (robust=x)
# assigns weight to visibilities based on uv coverage.
# Natural weighting reduces the noise but degrades the
# beamshape. 
# We can also change the fwhm parameter (fwhm=theta_maj,theta_min, pa); 
# the SNR in the output is optimized when the parameter is set to the
# FWHM of features in the image.  Careful: In deep cleaning, you 
# can deconvolve noise this way too and make what looks like a convincing
# detection. One may set cutoff to stop the cleaning.
# The cutoff can be chosen based on the theoretical noise
# calculated from invert. Try cutoff= 2-3 sigma.

  invert vis=$source.hco+vel.uv.lin map=$source.hco+.map \
        beam=$source.hco+.beam \
        imsize=512,512 cell=0.6 sup=0 \
        line=vel,10,-150,2

# For S1, line=vel,10,-160,5, for S2 line=vel,10,-120,5 (we know
# where the features are from a previous CARMA observation)

  \rm -r $source.hco+.icmp
  clean map=$source.hco+.map beam=$source.hco+.beam \
        out=$source.hco+.icmp gain=0.08 \
        cutoff=.4 niters=5000 \
	region=arcsec,'boxes(-25,-25,25,25)'

  \rm -r $source.hco+.icln
  restor model=$source.hco+.icmp beam=$source.hco+.beam \
        map=$source.hco+.map  \
        out=$source.hco+.icln

  cgdisp in=$source.hco+.icln type=c \
        region=arcsec,'boxes(-25,-25,25,25)' \
        xybin=1,1 device=/xs nxy=5,5 options=full,beambr,trlab,3val \
        labtyp=arcsec,arcsec cols1=7 csize=0.5,0.5,0.5 \
	slev=a,0.2 levs1=-5,3,4,8,16,32,64,128,256 

  echo "statistics on the channel images"
  imstat in=$source.hco+.icln region=quarter \
         axes=ra,dec
  imstat in=$source.hco+.icln region='boxes(30,30,60,60)' \
         axes=ra,dec
  echo "Press  to continue."
  $<
  cgdisp in=$source.hco+.icln type=pixel \
        region=arcsec,'boxes(-25,-25,25,25)' \
        xybin=1,1 device=/xs nxy=3,3 options=beambr,wedge,trlab,3val \
        labtyp=arcsec,arcsec range=0,0,lin,2 cols1=7 csize=0.5,0.5,0.5
  echo "Press  to continue."
  $<

# \rm -r $source.hco+.mom
# moment in=$source.hco+.icln region='image(3,6)' out=$source.hco+.mom axis=3 \
#  mom=0
# cgdisp in=$source.hco+.mom device=/xs options=3val labtyp=arcsec

  \rm -r $source.hco+.m0
  moment in=$source.hco+.icln region='quarter(1,20)' \
        out=$source.hco+.m0 axis=3 \
        mom=0 clip=-1.2,1.2
  smacgdisp in=$source.hco+.m0 type=pixel \
        region=arcsec,'boxes(-30,-30,30,30)' \
        xybin=1,1 device=/xs \
        nxy=1,1 options=beambr,wedge,trlab,3val \
        labtyp=arcsec,arcsec range=0,0,lin,2 cols1=7 csize=1,1,1
  $<
  \rm -r $source.hco+.m1
  moment in=$source.hco+.icln region='quarter(1,20)' \
        out=$source.hco+.m1 axis=3 \
        mom=1 clip=-1.2,1.2
  smacgdisp in=$source.hco+.m1 type=pixel \
        region=arcsec,'boxes(-30,-30,30,30)' \
        xybin=1,1 device=/xs \
        nxy=1,1 options=beambr,wedge,trlab,3val \
        labtyp=arcsec,arcsec range=-200,-100,lin,2 cols1=7 csize=1,1,1

# check the spectum from the image cube 
# imspect in=$source.hco+.icln region=arcsec,'box(-3,-3,3,3)' device=/xs  

  echo " "
  echo "The end"
  echo " "
exit
#jhz: updated the name of the script