#!/bin/csh -f
# handling standard 2GHz uniform data with bandpass ripples.
#
# Responsible: Jun-Hui Zhao (SAO)
# Version 3.3
# jhz 2005-10-20 Prepared for demo
#                for the reduction of SMA data.
#                for H30a line observations on Sept11_05 (Pamela Klaassen)
# jhz 2005-11-01 made a demo on jove, a linux box at CFA in real time.
#                the demo was broadcast to Hilo site.
# jhz 2005-11-02 corrected for few typos and added a few comments
#                to address the questions received from the audience.
# jhz 2010-2-28  test Miriad_beta for reduction of old SMA data
# jhz 2011-3-08  producing illustration figures and removing
#                the redundant text that was used in the live demo in 
#                2005 fall. Create ch0 data for continuum imaging.
# jhz 2011-4-25  add uvflag to flag the bad points that were left
#                for exercise of using smablflag
# jhz 2011-6-9   change the data path to CF where the test data 
#                located.
# jhz 2011-8-8   added goto TAG
# testing data: 050911_04:21:35 from SMA archive.
#
# demo illustration figures
# fig1.ps = continuum image of NGC7538 IRS1 
# fig2.ps = spectrum from the vis data of NGC7538 IRS1
# fig3.ps = spectrum from H30a line image cube NGC7538 IRS1
# fig4.ps = image of integrated H30a line flux NGC7538 IRS1
set darea=/sma/SMAusers/jzhao/data/050911_04:21:35
goto SMALOD
#goto INSPECT
#goto SMAFIX
#goto BPASS
#goto GCAL
#goto SMAFLUX
#goto CMAP
#goto UVLIN
#goto LMCUBE
#goto MOMENT
SMALOD:
echo "1. LOADING THE DATA" 
\rm -r  050911H30a_rx0.usb 
smalod in=$darea out=050911H30a rxif=0 \
       sideband=1 rsnchan=128 readant=8 nscans=1,
INSPECT:
echo "2. INSPECT THE RAW UV DATA"
echo "List the observing scans in the uv data"
echo "uvindex vis=050911H30a_rx0.usb"
$<
uvindex vis=050911H30a_rx0.usb
echo "Plot the visibilities of the raw data"
$<
smauvplt vis=050911H30a_rx0.usb \
	axis=time,ampl device=/xs nxy=3,5 \
         yrange=0,300
echo "Plot the Tsys measurements"
$<
smavarplt vis=050911H30a_rx0.usb device=/xs \
          xaxis=time \
          yaxis=systemp  nxy=2,4
echo "Filter out the bad tsys"
$<
smacheck vis=050911H30a_rx0.usb var=systemp flagval=flag \
        range=0,500
echo "Check up after Tsys flagging"
$<
smavarplt vis=050911H30a_rx0.usb device=/xs \
	  xaxis=time \
	  yaxis=systemp  yrange=50,150 nxy=2,4
SMAFIX:
$<
echo "3. FRINGE AMPLITUDE CORRECTION"
echo "Tsys Corrections for fringe amplitude"
$<
\rm -r 050911H30a_rx0.usb.tsys
smafix vis=050911H30a_rx0.usb out=050911H30a_rx0.usb.tsys \
	device=/xs xaxis=time \
       	yaxis=systemp nxy=2,4   \
       	options=tsyscorr
echo "Inspect vis data and flagging bad data"
uvflag vis=050911H30a_rx0.usb.tsys select='ant(3)(4),time(12:49:30,12:50:00)' \
	flagval=flag
uvflag vis=050911H30a_rx0.usb.tsys select='ant(4)(5),time(10:26:30,10:27:00)' \
        flagval=flag
uvflag vis=050911H30a_rx0.usb.tsys select='ant(4)(5),time(13:07:00,13:07:30)' \
        flagval=flag
uvflag vis=050911H30a_rx0.usb.tsys select='ant(4)(5),time(13:17:30,13:18:00)' \
        flagval=flag
uvflag vis=050911H30a_rx0.usb.tsys select='ant(4)(5),time(13:19:00,13:19:10)' \
        flagval=flag
uvflag vis=050911H30a_rx0.usb.tsys select='ant(4)(5),time(13:25:45,13:26:00)' \
        flagval=flag
uvflag vis=050911H30a_rx0.usb.tsys select='ant(4)(5),time(13:47:45,13:48:10)' \
        flagval=flag
smablflag vis=050911H30a_rx0.usb.tsys axis=time,ampl device=/xs
$<
echo "Inspect spectral data"
echo "for 3c454.3"
smauvspec vis=050911H30a_rx0.usb.tsys \
	select='source(3c454.3)' \
        interval=1000  \
        axis=freq,both device=/xs nxy=1,5
echo "for Neptune"
$<
smauvspec vis=050911H30a_rx0.usb.tsys \
        select='source(neptune)' \
        interval=1000  \
        axis=freq,both device=/xs nxy=1,5
BPASS:
echo "4. BANDPASS CALIBRATIONS"
echo "4.1. Correction for Possible Bandpass Ripples "
echo "Take the data of neptune scan from the multiple source file"
echo "into a single source file sept11bp_1 (in the time range 4:45-5:40);"
echo "take the two scans data of 3c454.3 from the multiple source"
echo "file into a single source files sept11bp_2(for the first scan"
echo "in the time range 6:00-6:30 and sept11bp_3 (for the second scan"
echo "in the time range 14:00-14:45) "
$<
\rm -r sept11bp_1 sept11bp_2 sept11bp_3
uvaver vis=050911H30a_rx0.usb.tsys select='source(neptune)' \
	out=sept11bp_1 options=nocal,nopass
uvaver vis=050911H30a_rx0.usb.tsys select='source(3c454.3),time(6:00,6:30)' \
       out=sept11bp_2 options=nocal,nopass
uvaver vis=050911H30a_rx0.usb.tsys select='source(3c454.3),time(14:00,14:45)' \
       out=sept11bp_3 options=nocal,nopass
echo "Solve for the antenna-based bandpass in the three time intervals"
#
# watch out for the effect from a broad, shallow absorption on Neptune at the CO transition.
#
$<
smamfcal vis=sept11bp_1 edge=5,5  \
        weight=2 refant=3 interval=600
$<
smamfcal vis=sept11bp_2 edge=5,5  \
        weight=2 refant=3 interval=600
$<
smamfcal vis=sept11bp_3 edge=5,5  \
        weight=2 refant=3 interval=600
echo "Check up the bandpass solutions"
echo "For time interval  1: Neptune (4:45-5:40)"
$<
smagpplt vis=sept11bp_1  device=/xs \
        yaxis=amp,phase \
        options=bandpass  polyfit=5 \
        nxy=2,4
echo "For the 2nd time interval: 3c454.3 (6:00-6:30)"
$<
smagpplt vis=sept11bp_2 device=/xs \
        yaxis=amp,phase \
        options=bandpass  polyfit=5 \
        nxy=2,4
echo "For the 3rd time interval: 3c454.3 (14:00-14:45)"
$<
smagpplt vis=sept11bp_3 device=/xs \
        yaxis=amp,phase \
        options=bandpass  polyfit=5 \
        nxy=2,4
# Note that the polynomial fitting curve will be not applied
# to bandpass solutions. For using the polynomial fits to
# replace the original solutions in the case of lower S/N,
# please look at the inline doc for smagpplt.
echo "Compare the two bpassies solved from different time"
echo "Plot the bandpass ratio of sept11bp_2 to sept11bp_3"
echo "For amplitude:"
$<
smagpplt vis=sept11bp_2,sept11bp_3 device=/xs \
        yaxis=amp \
        options=bandpass,ratio polyfit=5 \
        nxy=2,3 yrange=0.5,1.5 
echo "For phase:"
$<
smagpplt vis=sept11bp_2,sept11bp_3 device=/xs \
        yaxis=phase \
        options=bandpass,wrap,ratio  polyfit=5 \
        nxy=2,3  yrange=-50,50
echo "Interpolate bandpass solutions from the three time intervals, applying to the data"
$<
\rm -r 050911H30a_rx0.usb.tsys.bp
smatbpass  bpfile=sept11bp nfiles=3 \
        vis=050911H30a_rx0.usb.tsys bptime=30 options=nocal \
        out=050911H30a_rx0.usb.tsys.bp
echo "Upto this point, the calibration is done for antenna-based."
echo "4.2. Corrections for Baseline-based Errors in Bandpass"
echo "Inspect spectral data again with applied antenna-based"
echo "bandpass for any baseline-based problems in bandpass "
$<
smauvspec vis=050911H30a_rx0.usb.tsys.bp \
        select='source(3c454.3),time(14:00,14:30)' \
        interval=1000   \
        axis=freq,both device=/xs nxy=2,2 yrange=0,40
# we do see amplitude offset on some baseline;
# also phase drift on baseline related to antenna 7
# for a number of chunks are large.
echo "Correction for baseline-based bandpass errors" 
echo "now we extract source 3c454.3 from the data with corrected" 
echo "for antenna-based bandpass. The source MUST be point source."
$<
\rm -r sept11bl
uvaver vis=050911H30a_rx0.usb.tsys.bp \
	select='source(3c454.3),time(14:00,14:45)' \
        out=sept11bl options=nocal,nopass
$<
\rm -r 050911H30a_rx0.usb.tsys.bp.bl
blcal vis=sept11bl,050911H30a_rx0.usb.tsys.bp \
	interval=5000 out=050911H30a_rx0.usb.tsys.bp.bl
echo "Inspect spectral data again with corrected bandpass"
$<
smauvspec vis=050911H30a_rx0.usb.tsys.bp.bl \
        select='source(3c454.3),time(14:00,14:45)'  \
        interval=1000  \
        axis=freq,both device=/xs nxy=3,5 
echo "flag bad visibilities"
$<
smablflag vis=050911H30a_rx0.usb.tsys.bp.bl \
        axis=time,ampl device=/xs 
echo "Plot overall baseline-based spectra"
$<
smauvspec vis=050911H30a_rx0.usb.tsys.bp.bl \
        select='source(n7538)' \
        interval=1000 hann=3 yrange=0,12 \
        axis=freq,both device=/xs nxy=1,1
echo "For the H30a  line in spectral chunks 19,20,21"
$<
smauvspec vis=050911H30a_rx0.usb.tsys.bp.bl \
        select='source(n7538),window(19,20,21)' \
        interval=1000  hann=5 yrange=0,12 \
        axis=freq,both device=/xs nxy=3,5
$<
smauvspec vis=050911H30a_rx0.usb.tsys.bp.bl \
        select='source(n7538),window(19,20,21),ant(1)(4)' \
        interval=1000  hann=5 yrange=0,12 \
        axis=freq,both device=/xs nxy=1,1
echo "For the CO2-1 line in spectral chunks 3,4,5"
$<
smauvspec vis=050911H30a_rx0.usb.tsys.bp.bl \
        select='source(n7538),window(3,4,5)' \
        interval=1000  hann=7 yrange=0,10 \
        axis=freq,both device=/xs nxy=3,5
#Note upto this point, we have done edge channel flag;
#The spectra on the overlapped channels between two 
#adjacent spectral chunks match very well over the line feature.
#The calibration for bandpass is successful.
echo "uvflag edge channels"
$<
uvflag vis=050911H30a_rx0.usb.tsys.bp.bl edge=3,3 \
	flagval=flag
echo "Check spectrum  on one of the baselines"
smauvspec vis=050911H30a_rx0.usb.tsys.bp.bl \
        select='source(n7538),window(19,20,21),ant(1)(4)' \
        interval=1000  \
        axis=freq,both device=/xs nxy=1,1 
GCAL:
echo "5. ANTENNA GAIN CALIBRATION"
echo "Up to this point, we have corrected for bandpass errors."
echo "Check and flag visibility "
$<
smauvplt  vis=050911H30a_rx0.usb.tsys.bp.bl \
	select='source(n7538,2202*,0102*)' \
	device=/xs options=nocal axis=time,amp nxy=3,5
echo "Split the gain calibrators from the multiple source file"
$<
\rm -r 2202+422.230409 0102+584.230409
uvsplit vis=050911H30a_rx0.usb.tsys.bp.bl \
	select='source(2202+422,0102+584)' options=nowindow
echo "rename the calibrator files"
$<
\rm -r 2202+422.tsys.bp.bl 0102+584.tsys.bp.bl
mv 2202+422.230409 2202+422.tsys.bp.bl 
mv 0102+584.230409 0102+584.tsys.bp.bl
echo "Solve for antenna gains"
echo "For 0102+584"
smamfcal vis=0102+584.tsys.bp.bl  \
select='source(0102+584)'  \
refant=3 interval=2 options=nopassol
echo "Plot the gain solutions"
$<
smagpplt vis=0102+584.tsys.bp.bl device=/xs \
        yaxis=amp,phase log=0102+584gain.log \
        options=gains  polyfit=6 \
        nxy=2,4
echo "For 2202+422"
$<
smamfcal vis=2202+422.tsys.bp.bl  \
select='source(2202+422)'  \
refant=3 interval=2 options=nopassol
echo "Plot the gain solutions"
$<
smagpplt vis=2202+422.tsys.bp.bl device=/xs \
        yaxis=amp,pha  \
        options=gains  polyfit=6 \
        nxy=2,4
echo "Merge the two gain tables"
$<
smagpplt vis=2202+422.tsys.bp.bl,0102+584.tsys.bp.bl device=/xs \
        yaxis=amp,pha  \
        options=gains,merge polyfit=8 \
        nxy=2,3
echo "We fitted the baselines for antenna position."
echo "The position offsets between the antenna position"
echo "used on-line and the fittted position using bee"
echo "in miriad are small (for antenna 4):"
echo "     Delta x = 0.0007 nanosec ( 0.21 mm)"
echo "     Delta y =-0.0007 nanosec (-0.21 mm)" 
echo "     Delta z =-0.0004 nanosec (-0.12 mm)"
# the gain offsets between two calibrators are not due to
# antenna position errors. The offsets are present not only 
# in phase but also in amplitude. They are likely due to pointing
# errors.
echo "Gpcopy the merged gain table back"
$<
gpcopy vis=0102+584.tsys.bp.bl out=050911H30a_rx0.usb.tsys.bp.bl \
     mode=copy options=nopass
#apply the gains to visibilities
\rm -r tmp.uv
uvaver vis=050911H30a_rx0.usb.tsys.bp.bl out=tmp.uv
\rm -r 050911H30a_rx0.usb.tsys.bp.bl
mv tmp.uv 050911H30a_rx0.usb.tsys.bp.bl
SMAFLUX:
echo "6. BOOTSTRAPE THE FLUX DENSITY SCALE"
echo "Up to this point, we have corrected for amplitue (with Tsys),"
echo "bandpass, antenna gains and baseline-based errors. "
echo "We now bootstrape the flux density scale from  Ceres"
echo "to the visibility data set."
$<
smaflux vis=050911H30a_rx0.usb.tsys.bp.bl \
	select='source(ceres)' mirhome=$MIR
#
#The bootstraping relys on the Planck whole disk brightness 
#temperature from the SMA planetary model.
#  
echo "Apply the gains to the data"
$<
\rm -r 050911H30a_rx0.usb.tsys.bp.bl.av
uvaver vis=050911H30a_rx0.usb.tsys.bp.bl out=050911H30a_rx0.usb.tsys.bp.bl.av
echo "Check for spectra of NGC7538"
echo "for overall"
$<
smauvspec vis=050911H30a_rx0.usb.tsys.bp.bl.av \
        select='source(n7538)' \
        interval=1000 hann=3 yrange=0,12 \
        axis=freq,both device=/xs nxy=3,5
echo "For H30a line"
$<
smauvspec vis=050911H30a_rx0.usb.tsys.bp.bl.av \
        select='source(n7538),window(19,20,21)' \
        interval=1000 hann=3 yrange=0,12 \
        axis=freq,both device=/xs nxy=3,5
echo "Uvsplit n7538" 
$<
\rm -r n7538.230409 n7538.h30a
uvsplit vis=050911H30a_rx0.usb.tsys.bp.bl.av select='source(n7538)' \
	options=nowindow
echo "Check up the rest frequency"
$<
prthd in=n7538.230409
echo "Rename the output file"
mv n7538.230409 n7538.h30a
echo "For continuum"
\rm -r n7538.usb.ch0
uvlin vis=n7538.h30a chans=1,3072 out=n7538.usb.ch0 \
        order=1 mode=chan0 options=nowin
echo "skip the following step for identification of line transitions "
echo "for CO2-1, using the following steps to reset the rest frequency"
echo "mv n7538.230409 n7538.co21"
CMAP:
echo "7. CONTINUUM IMAGE"
$<
\rm -r n7583con.map n7583con.beam n7583con.icmp
echo "Make a dirty image"
echo "using the channel data between CO line and H30a line,"
echo "i.e., from spectral chunks 6-18"
invert vis=n7538.usb.ch0 map=n7583con.map beam=n7583con.beam \
        imsize=256,256 cell=0.2 sup=0 
echo "Display the dirty image"
$<
cgdisp in=n7583con.map type=pixel region=arcsec,'boxes(-10,-10,10,10)' \
        xybin=1,1 device=/xs nxy=1,1 options=full,beambr,wedge,trlab,3val \
        labtyp=arcsec,arcsec range=0,3,lin,2 cols1=7 csize=0.5,0.5,0.5
echo "Clean  the dirty image"
$<
clean map=n7583con.map beam=n7583con.beam out=n7583con.icmp gain=0.08 \
        cutoff=0.001 niters=5000 \
        region=arcsec,'boxes(-1,-1,1,1)'
$<
echo "Restore cleaned image"
\rm -r n7583con.icln
restor model=n7583con.icmp beam=n7583con.beam map=n7583con.map  \
        out=n7583con.icln
$<
echo "Display the restored clean image"
cgdisp in=n7583con.icln type=pixel region=arcsec,'boxes(-10,-10,10,10)' \
        xybin=1,1 device=/xs nxy=1,1 options=beambr,wedge,blacklab \
        labtyp=arcsec,arcsec range=0.1,3,lin,2 cols1=7 csize=1,1
$<
cgdisp in=n7583con.icln type=pixel region=arcsec,'boxes(-10,-10,10,10)' \
        xybin=1,1 device=fig1.ps/cps nxy=1,1 options=beambr,wedge,blacklab \
        labtyp=arcsec,arcsec range=0.1,3,lin,2 cols1=7 csize=1,1
#by quick looking at the continuum image
#the gain calibration went very well.
#For better image, need to take out the line contamination 
#from the unknown line features between h30a and CO2-1.
UVLIN:
echo "8. LINE/CONTINUUM SEPARATION"
echo "plot spectrum of the usb data"
smauvspec vis=n7538.h30a \
        select='source(n7538)' hann=5 \
        interval=10000  options=nobase,avall \
        axis=freq,both device=/xs nxy=1,1 yrange=0,5 
$<
smauvspec vis=n7538.h30a \
        select='source(n7538)' hann=5 \
        interval=10000  options=nobase,avall \
        axis=freq,both device=fig2.ps/cps nxy=1,1 yrange=0,5
echo "Plot H30a spectrum as function of velocity"
$<
smauvspec vis=n7538.h30a \
        select='source(n7538),window(19,20,21)' \
        interval=10000  options=nobase,avall \
        axis=vel,both device=/xs nxy=1,1 yrange=0,10
echo "Convert sky-frequency to velocity -220 to 100 at"
echo "resolution of 1 km/s"
$<
\rm -r n7538.h30a.vel
uvaver vis=n7538.h30a line=velocity,320,-220,1 out=n7538.h30a.vel
echo "Plot the velocity spectrum"
$<
smauvspec vis=n7538.h30a.vel \
        interval=10000  options=nobase,avall \
        axis=chan,both device=/xs nxy=1,1 yrange=0,10
echo "Take out the continuum using linear fitting to" 
echo "the line free channels"
$<
\rm -r n7538.h30a.vel.lin
uvlin vis=n7538.h30a.vel chans=1,30,250,280 out=n7538.h30a.vel.lin \
	order=1 mode=line
echo "Plot the continuum-free spectrum"
$<
smauvspec vis=n7538.h30a.vel.lin options=nobase,avall  \
        interval=1000 \
        axis=vel,both device=/xs nxy=1,1
LMCUBE:
echo "9. CONSTRUCT LINE IMAGE CUBE" 
echo "Make a dirty image cube"
$<
\rm -r n7538.h30alin.map n7538.h30alin.beam n7538.h30alin.icmp
invert vis=n7538.h30a.vel.lin map=n7538.h30alin.map beam=n7538.h30alin.beam \
        imsize=256,256 cell=0.2 sup=0 \
        line=channel,175,75,1
echo "Clean the dirty image cube"
$<
clean map=n7538.h30alin.map beam=n7538.h30alin.beam \
	out=n7538.h30alin.icmp gain=0.08 \
        cutoff=0 niters=2500 \
	region=arcsec 'boxes(-1,-1,1,1)(1,175)'
echo "Restore the cleaned image cube"
$<
\rm -r n7538.h30alin.icln
restor model=n7538.h30alin.icmp beam=n7538.h30alin.beam map=n7538.h30alin.map  \
        out=n7538.h30alin.icln
echo "Display the restired clean image cube"
$<
cgdisp in=n7538.h30alin.icln type=c region=arcsec,'boxes(-5,-5,5,5)(1,175)' \
        xybin=1,1 device=/xs nxy=10,10 options=full,beambr,trlab,3val \
        labtyp=arcsec,arcsec cols1=7 csize=0.5,0.5,0.5 \
	slev=a,0.1 levs1=-5,3,4,8,16,32,64,128,256 
echo "statistics on the channel images"
echo "On source region"
$<
imstat in=n7538.h30alin.icln region=quarter \
        axes=ra,dec
echo "Off source region"
$<
imstat in=n7538.h30alin.icln region='boxes(30,30,60,60)(1,175)' \
        axes=ra,dec
echo "make spectrum from the line image cube."
imspect in=n7538.h30alin.icln region=arcsec,'boxes(-0.5,-0.5,0.5,0.5)(1,175)' \
        device=/xs yrange=-0.5,3 xaxis=channel
imspect in=n7538.h30alin.icln region=arcsec,'boxes(-0.5,-0.5,0.5,0.5)(1,175)' \
	device=fig3.ps/ps yrange=-0.5,3
$<
MOMENT:
\rm -r n7538.flux
moment in=n7538.h30alin.icln region='quarter(10,165)' out=n7538.flux \
mom=0 clip=-10,0.2
$<
cgdisp in=n7538.flux,n7538.flux type=c,p \
        region=arcsec,'boxes(-5,-5,5,5)' \
        xybin=1,1 device=/xs nxy=1,1 options=beambr,wedge,blacklab \
        labtyp=arcsec,arcsec range=1,300,lin,2 cols1=7 csize=1.2,1.2 \
        slev=p,1 levs1=10,15,25,35,45,55,65,75,85,95
$<
cgdisp in=n7538.flux,n7538.flux type=c,p \
        region=arcsec,'boxes(-5,-5,5,5)' \
        xybin=1,1 device=fig4.ps/cps nxy=1,1 options=beambr,wedge,blacklab \
        labtyp=arcsec,arcsec range=1,300,lin,2 cols1=7 csize=1.2,1.2 \
        slev=p,1 levs1=10,15,25,35,45,55,65,75,85,95	
echo "The end"
exit