#! /bin/csh -f
#jhz 2015-10-20 create
#jhz 2016-01-12 modify
#jhz 2016-01-27 edited for CO(2-1)
#jhz 2016-01-29 finalize and document
#Imaging-proces: for calibrated CO(2-1) spectral line data 
#                produced from SWARM and ASIC, SMA hybrid correlator
#                with Miriad4GHz software -
#                <a href="https://www.cfa.harvard.edu/sma/miriad/miriad4GHz/SMAMiriad4GHz.html">(SMA-Miriad 1.5.1)</a>
#
echo "                                                                      "
echo "######################################################################"
echo "# imagingln_CO21.csh: script to image CO(2-1) spectral line extracted "
echo "# from a calibrated subset of SMA data produced from Pre-process 2.   "
echo "# Users are required to set up the line parameters concerning the     "
echo "# target line and target sources (fields). The spectral line imaging  "
echo "# process are recommended to be carried out in Miriad. The spectral   "
echo "# line uv data set should be compatiable with any reliable versions   "
echo "# distributed by CARMA, ATNF, and SAO concerning imaging. And FITS    "
echo "# outputs are produced for the spectral line cube in the end of the   "
echo "# imaging process.                                                    "
echo "# It is option if users choose to make UVFITS output for the specific "
echo "# line uv data and image elsewhere (CASA, Miriad, AIPS and etc.).     "
echo "######################################################################"
#
#Target                (font) - 
#Bandpass calibrator   (bcal) -
#Gain calibrator       (gcal, or gcal1, gcal2) -
#Flux scale calibrator (fcal) -
#
echo "                                                                      "
echo "######################################################################"
echo "# Users need to carefully review the following project-depended setup "
echo "# from lines 37 - 134 prior to executing this script                   "
echo "######################################################################"
echo "                                                                      "
#
#
##set the name of the SMA raw data
#
set file = 160115_03:45:40
#
echo $file   
#
set dt          = 160115
echo ' '
echo 'prefix of data files = ' $dt
echo ' '
#
set hour        = 03:45:40 
set sb          = usb
set nn          = 1 
set rxid        = 0
set rx          = rx${rxid}
set adge        =  7 
set edge        = 500
set nchan_as    = 6130
set maxtsys     = 500
set fname       = "$dt""_""$rx".$sb
set oname       = "$dt""_""$rx"
set fcal        = uranus 
set gcal1       = 0721+713
set gcal2       = 0841+708 
set gcal        = 0721+713  
set bcal        = 3C84
set blcal       = $bcal 
set font        = N2403 
set refant      = 1 
set avetime     = 5
set source1     = ${font}dc1
set source2     = ${font}fd2
set source3     = ${font}fd3
set source4     = ${font}fd4
set source5     = ${font}fd5
#set twogcal    = F    % T two gain calibrators' scheme
set twogcal     = T
set DATPATH     = ../ucla-pp2/ 
                # 
                # The data path contains the calibrated output data in pre-process 2,
                # 160115_rx0.usb.asic.tsys.bp.a                      
                #                     
set fname = ${DATPATH}/${fname}
###########################################
#set spectral line parameters
###########################################
#
#sline is the name of the line transition
#
set sline = CO21
#
#rfreq is the CO(2-1) rest frequency in GHz
#
set rfreq = 230.5379700
#
#input calibrated uv file contains the line
#
set uvin = $fname.asic.tsys.bp.a
#
#target sources
#
set targets = $source1,$source2,$source3,$source4,$source5
#
#output calibrated uv file contains the line
#
set uvout = ${dt}.${sline}.uv
#
#line-free windows 
#
set lfwin = 1,3000,3200,6144
#
#mosaic fields
#
set mosaicA = $source1,$source2,$source3
set mosaicB = $source4,$source5
#
#image cube name
#
set imnameA = ${font}.${sline}.fd123
set imnameB = ${font}.${sline}.fd45
#
#linetyp for imaging line cube
#
set linetype = vel,50,0,5
#
#imaging size and cell size
#
set ims   = 512
set cells = 0.5
#
#clean iter, and fwhm
#
set niter = 2500
set fwhm  = 4.5,4.5
set rms   = 0.025
#
############################################################################
#Imaging-process handles -
############################################################################
#goto INSPECT         # inspect the data calibrated in pre-processing 2
#goto UVSPC-ASIC      # inspect calibrated uv spectra for ASIC data
#goto UVLIN           # separate the line from the vis data
                      # cut the specific line data from the calibrated data set
#goto BLFLAG          # flag bad data points for ASIC chunks
                      #######################################################
                      # below are special for the SMA project 2015B-S041 
#goto LNMAP-FD123     # do mosaic-imaging for fields 1 2 3
#goto CDISP-FD123     # display images of fields 1 2 3
#goto STATS-FD123     # do statistics for the image line cube of fields 1 2 3
#goto LNMAP-FD45      # do mosaic-imaging for fields 4 5
#goto CDISP-FD45      # display images of fields 4 5
#goto STATS-FD45      # do statistics for the image line cube of fields 4 5
                      #######################################################
#goto UVFITS          # export UVFITS from the calibrated specific line uvdata
#goto FITS            # export FITS from the image line cube
#goto TAR             # make tarballs for the miriad uv files and line image cubes
#goto ZAP             # clean up the working files
#
$<
INSPECT:
echo 'Inspect: ASIC'
uvindex  vis=$fname.asic.tsys.bp.a

$<
UVSPC-ASIC:
$<
echo 'Check spectra: ASIC -' $bcal
smauvspec vis=$fname.asic.tsys.bp.a \
        select='source('$bcal')' \
        interval=1000 hann=5 options=nobase,avall \
        axis=freq,both device=/xs nxy=1,1
$<
echo 'Check spectra: ASIC -' $gcal
smauvspec vis=$fname.asic.tsys.bp.a \
        select='source('$gcal')' \
        interval=1000 hann=5 options=nobase,avall \
        axis=freq,both device=/xs nxy=1,1
$<
echo 'Check spectra: ASIC -' $fcal
smauvspec vis=$fname.asic.tsys.bp.a \
        select='source('$fcal')' \
        interval=1000 hann=5 \
        axis=freq,both device=/xs nxy=1,1 options=nobase,avall
$<
echo 'Check spectra: ASIC -' $targets
smauvspec vis=$fname.asic.tsys.bp.a \
        select='source('$targets')' \
        interval=1000 hann=5 \
        axis=freq,both device=/xs nxy=1,1 options=nobase,avall

$<
UVLIN:
echo 'Check spectra of uvdata: ' $uvin
echo 'Reurn for fringe amp&pha vs sky-frequnecy -'
$<
smauvspec vis=$uvin \
        select='source('$targets')' \
        interval=1000 hann=5 options=nobase,avall \
        axis=freq,both device=/xs nxy=1,1

echo 'Reurn for fringe amp&pha vs LSR-velocity -'
$<
smauvspec vis=$uvin \
        select='source('$targets')' \
        interval=1000 hann=5 options=nobase,avall \
        axis=vel,both device=/xs nxy=1,1

echo 'Reurn for fringe amp&pha channel number -'
$<
smauvspec vis=$uvin \
        select='source('$targets')' \
        interval=1000 hann=5 options=nobase,avall \
        axis=chan,both device=/xs nxy=1,1
$<
\rm -r $uvout 
echo 'line-free windows -' $lfwin
echo 'Make sure the line-free windows are correct'
echo 'Reurn for continuum subtraction -'
$<   
uvlin vis=$uvin chans=$lfwin order=1 \
       mode=line select='source('$targets')' \
       out=$uvout options=nowin
$<
echo 'check spectra of continuum-subtracted uvdata: ' $uvout
echo 'Reurn for fringe amp&pha vs sky-frequnecy -'
$<
smauvspec vis=$uvout \
        select='source('$targets')' \
        interval=1000 hann=5 options=nobase,avall \
        axis=freq,both device=/xs nxy=1,1
echo 'Reurn for fringe amp&pha vs LSR-velocity -'
$<
smauvspec vis=$uvout \
        select='source('$targets')' \
        interval=1000 hann=5 options=nobase,avall \
        axis=vel,both device=/xs nxy=1,1
echo 'Reurn for fringe amp&pha channel number -'
$<
smauvspec vis=$uvout \
        select='source('$targets')' \
        interval=1000 hann=5 options=nobase,avall \
        axis=chan,both device=/xs nxy=1,1
$<
BLFLAG:
echo 'Baseline-based flagging -' $uvout
smablflag vis=$uvout axis=time,ampl device=/xs

echo 'Reurn for fringe amp&pha vs LSR-velocity '
echo 'with a larger size of hanning smoothing window -'
$<
smauvspec vis=$uvout \
        select='source('$targets')' \
        interval=1000 hann=15 options=nobase,avall \
        axis=vel,both device=/xs nxy=1,1

LNMAP-FD123:
echo 'Reurn for mosaic imaging group A of pointing fields for linetype -' $linetype 
$<
\rm -r $imnameA.map $imnameA.beam $imnameA.icmp
invert vis=$uvout map=$imnameA.map beam=$imnameA.beam \
	imsize=$ims,$ims cell=$cells sup=0 options=systemp,mosaic \
        line=$linetype select='source('$mosaicA')'
echo 'Reurn for examing dirty image line cube -'
$<
cgdisp in=$imnameA.map type=pixel region=arcsec,'boxes(-125,-125,75,100)(2,50)' \
        xybin=1,1 device=/xs nxy=7,7 options=full,beambr,wedge,trlab,3val \
        labtyp=arcsec,arcsec range=0,0,lin,2 cols1=7 csize=0.5,0.5,0.5
echo 'Reurn to clean dirty image line cube -'
$<
clean map=$imnameA.map beam=$imnameA.beam out=$imnameA.icmp gain=0.08 \
        cutoff=$rms niters=$niter \
        region=arcsec,'boxes(-125,-125,75,100)(1,50)'
echo 'Reurn to restore image cube with a circular beam=' $fwhm 
$<
\rm -r $imnameA.icln
restor model=$imnameA.icmp beam=$imnameA.beam map=$imnameA.map  \
        out=$imnameA.icln fwhm=$fwhm
$<
CDISP-FD123:
echo 'Reurn to exam the cleaned  line cube -'
$<
cgdisp in=$imnameA.icln type=pixel \
        region=arcsec,'boxes(-125,-125,75,100)(2,50)' \
        xybin=1,1 device=/xs \
        nxy=7,7 options=beambr,wedge,blacklab,3val \
        labtyp=arcsec,arcsec range=0,0.5,lin,2 cols1=7 csize=0.5,0.5,0.5
echo 'Reurn to exam the clouds 1 -'
$<
cgdisp in=$imnameA.icln type=pixel \
        region=arcsec,'boxes(-95,-5,0,90)(22,30)' \
        xybin=1,1 device=/xs \
        nxy=3,3 options=beambr,wedge,blacklab,3val \
        labtyp=hms,dms range=0,0.5,lin,2 cols1=7 csize=0.5,0.5,0.5
echo 'Reurn to exam the clouds 2 -'
$<
cgdisp in=$imnameA.icln type=pixel \
        region=arcsec,'boxes(-45,-50,50,45)(30,38)' \
        xybin=1,1 device=/xs \
        nxy=3,3 options=beambr,wedge,blacklab,3val \
        labtyp=hms,dms range=0,0.5,lin,2 cols1=7 csize=0.5,0.5,0.5
echo 'Reurn to exam the clouds 3 -'
$<
cgdisp in=$imnameA.icln type=pixel \
        region=arcsec,'boxes(-80,-105,50,5)(35,46)' \
        xybin=1,1 device=/xs \
        nxy=3,4 options=beambr,wedge,blacklab,3val \
        labtyp=hms,dms range=0,0.3,lin,2 cols1=7 csize=0.5,0.5,0.5
$<
STATS-FD123:
echo 'Reurn to do statistics for  the image cube -'
$<
imstat in=$imnameA.icln region='arcsec,box(-5,-5,5,5)'
$<
imstat in=$imnameA.icln region='arcsec,box(-55,-55,-70,-70)'
$<
imstat in=$imnameA.icln region='arcsec,box(-55,55,-70,70)'
exit

LNMAP-FD45:
echo 'Reurn for mosaic imaging group B of pointing fields for linetype -' $linetype
\rm -r $imnameB.map $imnameB.beam $imnameB.icmp
invert vis=$uvout map=$imnameB.map beam=$imnameB.beam \
        imsize=$ims,$ims cell=$cells sup=0 options=systemp,mosaic \
        line=$linetype select='source('$mosaicB')'
echo 'Reurn for examing dirty image line cube -'
$<
cgdisp in=$imnameB.map type=pixel region=arcsec,'boxes(-125,-125,75,100)(2,50)' \
        xybin=1,1 device=/xs nxy=7,7 options=beambr,wedge,trlab,3val \
        labtyp=arcsec,arcsec range=0,0,lin,2 cols1=7 csize=0.5,0.5,0.5
echo 'Reurn to clean the dirty image cube -'
$<
clean map=$imnameB.map beam=$imnameB.beam out=$imnameB.icmp gain=0.08 \
        cutoff=$rms niters=$niter \
        region=arcsec,'boxes(-50,-75,75,50)(1,50)'
echo 'Reurn to restore image cube with a circular beam=' $fwhm
$<
\rm -r $imnameB.icln
restor model=$imnameB.icmp beam=$imnameB.beam map=$imnameB.map  \
        out=$imnameB.icln fwhm=$fwhm
$<
CDISP-FD45:
echo 'Reurn to exam the cleaned  line cube -'
cgdisp in=$imnameB.icln type=pixel \
        region=arcsec,'boxes(-75,-100,75,75)(2,50)' \
        xybin=1,1 device=/xs \
        nxy=7,7 options=beambr,wedge,blacklab,3val \
        labtyp=arcsec,arcsec range=0,0.5,lin,2 cols1=7 csize=0.5,0.5,0.5
$<
echo 'Reurn to exam the clouds 4 -'
cgdisp in=$imnameB.icln type=pixel \
        region=arcsec,'boxes(-50,-60,60,50)(26,41)' \
        xybin=1,1 device=/xs \
        nxy=4,4 options=beambr,wedge,blacklab,3val \
        labtyp=hms,dms range=0,0.5,lin,2 cols1=7 csize=0.5,0.5,0.5
$<
STATS-FD45:
echo 'Reurn to do statistics for  the image cube -'
$<
imstat in=$imnameB.icln region='arcsec,box(-5,-5,5,5)'
$<
imstat in=$imnameB.icln region='arcsec,box(-40,-40,-50,-50)'
$<
imstat in=$imnameB.icln region='arcsec,box(40,40,50,50)'

echo ''
echo 'Done with imaging process!'
echo ''
echo 'Congratulation if you detected the line that you proposed to observe!'
echo ''
echo 'You may need to go back the editing & calibration process to check'
echo 'your setup if the line has not been detected.'
echo ''
exit

UVFITS:
echo 'export uvfits from ' $uvout
\rm -r $uvout.uvfits
fits in=$uvout op='uvout' out=$uvout.uvfits
echo 'done with uvfits output file =' $uvout.uvfits
$<
FITS:
echo 'export fits from ' $imnameA.icln ' & ' $imnameB.icln
\rm -r $imnameA.icln.fits
fits in=$imnameA.icln op='xyout' out=$imnameA.icln.fits
$<
\rm -r $imnameB.icln.fits
fits in=$imnameB.icln op='xyout' out=$imnameB.icln.fits
echo 'done with fits output files;'
echo 'fits image cubes are ' $imnameA.icln.fits $imnameB.icln.fits
echo ''
echo 'You can continue do analysis with the fits data somewhere else.'
exit
TAR:
echo 'Make gzipped tarballs for the Miriad image cubes -'
\rm $uvout.tar.gz
tar -cvzf $uvout.tar.gz $uvout
\rm $imnameA.tar.gz
tar -cvzf $imnameA.tar.gz $imnameA.icln $imnameA.map $imnameA.beam $imnameA.icmp
\rm $imnameB.tar.gz
tar -cvzf $imnameB.tar.gz $imnameB.icln $imnameB.map $imnameB.beam $imnameB.icmp
echo 'done with tarballs!'
exit
ZAP:
\rm -r $imnameA.map $imnameA.beam $imnameA.icmp 
\rm -r $imnameB.map $imnameB.beam $imnameB.icmp
echo 'done with clean up the working files!'
exit
