#!/bin/csh -f
# handling standard 2GHz uniform data with bandpass ripples.
#
# Responsible: Jun-Hui Zhao (SAO)
# Version 3.0
# 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.
goto continue
#
continue:
echo "Start"
set darea=/pool/tdc4/miriad/050911_04:21:35
echo " "
echo "End-to-end procedure for the reduction of SMA data"
echo "taken from the observations of HII region NGC 7538"
echo "on 2005-9-11 by Pamela Klaassen and Eric Keto etal"
echo " "
echo "ls -l 050911.tar"
ls -l 050911.tar
echo "The file with a size of 1 Gbytes containing"
echo "both side bands. The H30alpha line contains"
echo "in the usb. Also the spectral resolution of"
echo "the original MIR data is 256 channels per "
echo "chunk. The loading is set to resample the "
echo "data to 128 channels per chunk (rsnchan=128)" 
sleep 5 
echo " "
echo "1. READ THE MIR DATA INTO MIRIAD"
echo " "
echo "smalod in=$darea \
      out=050911H30a rxif=0 \
      sideband=1 rsnchan=128 readant=8 nscans=1, "
sleep 1
\rm -r  050911H30a_rx0.usb 
smalod in=$darea out=050911H30a rxif=0 \
       sideband=1 rsnchan=128 readant=8 nscans=1,
sleep 5
echo " "
echo "2. INSPECT THE RAW UV DATA"
echo " "
echo "List the observing scans in the uv data"
echo " "
echo "uvindex vis=050911H30a_rx0.usb"
echo " "
sleep 5 
uvindex vis=050911H30a_rx0.usb
sleep 10
echo " "
echo "Plot the visibilities of the raw data"
echo " "
echo "smauvplt vis=050911H30a_rx0.usb \
        axis=time,ampl device=/xs nxy=3,5 \
         yrange=0,300"
sleep 3
smauvplt vis=050911H30a_rx0.usb \
	axis=time,ampl device=/xs nxy=3,5 \
         yrange=0,300
sleep 10
echo " "
echo "Filter out the 5 sec pointing scans (option)"
echo " "
echo "skip smacheck vis=050911H30a_rx0.usb var=inttime flagval=flag \
        range=10,60 "
echo "The pointing data on 3C454.3 will be used for solving"
echo "for antenna-based bandpass in an earlier time interval"
echo "to check up the bandpass ripples." 
echo " "
echo "Plot the Tsys measurements"
echo " " 
echo "smavarplt vis=050911H30a_rx0.usb device=/xs \
          xaxis=time \
          yaxis=systemp  nxy=2,4"
sleep 15 
smavarplt vis=050911H30a_rx0.usb device=/xs \
          xaxis=time \
          yaxis=systemp  nxy=2,4
sleep 10
echo " "
echo "Filter out the bad tsys"
echo " "
echo "smacheck vis=050911H30a_rx0.usb var=systemp flagval=flag \
        range=0,500"
sleep 5  
smacheck vis=050911H30a_rx0.usb var=systemp flagval=flag \
        range=0,500
sleep 10
echo " "
echo "Check up after Tsys flagging"
echo " "
echo "smavarplt vis=050911H30a_rx0.usb device=/xs \
          xaxis=time \
          yaxis=systemp  yrange=50,150 nxy=2,4"
sleep 5 
smavarplt vis=050911H30a_rx0.usb device=/xs \
	  xaxis=time \
	  yaxis=systemp  yrange=50,150 nxy=2,4
sleep 10
echo " "
echo "3. FRINGE AMPLITUDE CORRECTION"
echo " "
echo "Tsys Corrections for fringe amplitude"
echo " " 
echo "smafix vis=050911H30a_rx0.usb out=050911H30a_rx0.usb.tsys \
	device=/xs xaxis=time yaxis=systemp nxy=2,4 \
	options=tsyscorr "
sleep 10 
\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
sleep 10
echo " "
echo "Inspect vis data and flagging bad data"
echo " "
echo "smablflag vis=050911H30a_rx0.usb.tsys axis=time,ampl device=/xs"
sleep 5 
smablflag vis=050911H30a_rx0.usb.tsys axis=time,ampl device=/xs
echo " "
echo "Inspect spectral data"
echo " "
echo "for 3c454.3"
echo " "
echo "smauvspec vis=050911H30a_rx0.usb.tsys \
        select='source(3c454.3)' \
        interval=1000  \
        axis=freq,both device=/xs nxy=1,5"
sleep 5 
smauvspec vis=050911H30a_rx0.usb.tsys \
	select='source(3c454.3)' \
        interval=1000  \
        axis=freq,both device=/xs nxy=1,5
echo " "
echo "for Neptune"
echo " "
echo "smauvspec vis=050911H30a_rx0.usb.tsys \
        select='source(neptune)' \
        interval=1000  \
        axis=freq,both device=/xs nxy=1,5"
sleep 5 
smauvspec vis=050911H30a_rx0.usb.tsys \
        select='source(neptune)' \
        interval=1000  \
        axis=freq,both device=/xs nxy=1,5
echo " "
echo "4. BANDPASS CALIBRATIONS"
echo " "
echo "4.1. Correction for Possible Bandpass Ripples "
echo " "
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) "
echo " "
echo "uvaver vis=050911H30a_rx0.usb.tsys \
	select='source(neptune)' \
        out=sept11bp_1 options=nocal,nopass"
echo "uvaver vis=050911H30a_rx0.usb.tsys \
        select='source(3c454.3),time(6:00,6:30)' \
        out=sept11bp_2 options=nocal,nopass"
echo "uvaver vis=050911H30a_rx0.usb.tsys \
	select='source(3c454.3),time(14:00,14:45)' \
        out=sept11bp_3 options=nocal,nopass"
sleep 15 
\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 " "
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.
#
echo "smamfcal vis=sept11bp_1 edge=5,5" 
echo "        weight=2 refant=3 interval=600 "
sleep 5
smamfcal vis=sept11bp_1 edge=5,5  \
        weight=2 refant=3 interval=600
echo " "
echo "smamfcal vis=sept11bp_2 edge=5,5" 
echo "weight=2 refant=3 interval=600"
sleep 5 
smamfcal vis=sept11bp_2 edge=5,5  \
        weight=2 refant=3 interval=600
echo " "
echo "smamfcal vis=sept11bp_3 edge=5,5" 
echo "	weight=2 refant=3 interval=600"
smamfcal vis=sept11bp_3 edge=5,5  \
        weight=2 refant=3 interval=600
echo " "
sleep 5 
echo "Check up the bandpass solutions"
echo " "
echo "For time interval  1: Neptune (4:45-5:40)"
echo " "
echo "smagpplt vis=sept11bp_1  device=/xs \
        yaxis=amp,phase \
        options=bandpass  polyfit=5 \
        nxy=2,4"
sleep 10 
smagpplt vis=sept11bp_1  device=/xs \
        yaxis=amp,phase \
        options=bandpass  polyfit=5 \
        nxy=2,4
echo " "
echo "For the 2nd time interval: 3c454.3 (6:00-6:30)"
echo "smagpplt vis=sept11bp_2  device=/xs \
        yaxis=amp,phase \
        options=bandpass  polyfit=5 \
        nxy=2,4"
sleep 3
smagpplt vis=sept11bp_2 device=/xs \
        yaxis=amp,phase \
        options=bandpass  polyfit=5 \
        nxy=2,4
echo " "
echo "For the 3rd time interval: 3c454.3 (14:00-14:45)"
echo "smagpplt vis=sept11bp_3  device=/xs \
        yaxis=amp,phase \
        options=bandpass  polyfit=5 \
        nxy=2,4"
sleep 3
continue:
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.
sleep 10 
echo " "
echo "Compare the two bpassies solved from different time"
echo "Plot the bandpass ratio of sept11bp_2 to sept11bp_3"
echo " "
echo "For amplitude:"
echo " "
echo "smagpplt vis=sept11bp_2,sept11bp_3 device=/xs \
        yaxis=amp \
        options=bandpass,ratio polyfit=5 \
        nxy=2,3 yrange=0.5,1.5"
sleep 3
smagpplt vis=sept11bp_2,sept11bp_3 device=/xs \
        yaxis=amp \
        options=bandpass,ratio polyfit=5 \
        nxy=2,3 yrange=0.5,1.5 
sleep 5
echo " "
echo "For phase:"
echo " "
echo "smagpplt vis=sept11bp_2,sept11bp_3 device=/xs \
        yaxis=phase \
        options=bandpass,wrap,ratio  polyfit=5 \
        nxy=2,3  yrange=-50,50 "
sleep 5 
smagpplt vis=sept11bp_2,sept11bp_3 device=/xs \
        yaxis=phase \
        options=bandpass,wrap,ratio  polyfit=5 \
        nxy=2,3  yrange=-50,50
echo " "
echo "Interpolate bandpass solutions from the three time intervals, applying to the data"
echo " "
echo "smatbpass  bpfile=sept11bp nfiles=3 \
        vis=050911H30a_rx0.usb.tsys bptime=30 options=nocal \
        out=050911H30a_rx0.usb.tsys.bp "
sleep 10 
\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 " "
echo "Upto this point, the calibration is done for antenna-based."
echo " "
echo "4.2. Corrections for Baseline-based Errors in Bandpass"
echo " "
echo "Inspect spectral data again with applied antenna-based"
echo "bandpass for any baseline-based problems in bandpass "
echo " "
echo "smauvspec vis=050911H30a_rx0.usb.tsys.bp\
        select='source(3c454.3),time(14:00,14:30)' \
        interval=1 \
        axis=freq,both device=/xs nxy=2,2 yrange=0,40 " 
sleep 15 
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
echo " "
# 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 " "
echo "now we extract source 3c454.3 from the data with corrected" 
echo "for antenna-based bandpass. The source MUST be point source."
echo " "
echo "uvaver vis=050911H30a_rx0.usb.tsys.bp \
        select='source(3c454.3),time(14:00,14:45)' \
        out=sept11bl options=nocal,nopass"
\rm -r sept11bl
uvaver vis=050911H30a_rx0.usb.tsys.bp \
	select='source(3c454.3),time(14:00,14:45)' \
        out=sept11bl options=nocal,nopass
echo " "
echo "blcal vis=sept11bl,050911H30a_rx0.usb.tsys.bp \
        interval=5000 out=050911H30a_rx0.usb.tsys.bp.bl"
echo " "
\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 " "
echo "Inspect spectral data again with corrected bandpass"
echo " "
echo "smauvspec vis=050911H30a_rx0.usb.tsys.bp.bl"
echo "select='source(3c454.3),time(14:00,14:45)'"
echo "interval=1000"
echo "axis=freq,both device=/xs nxy=3,5"
sleep 3
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 " "
echo "flag bad visibilities"
echo " "
echo "smablflag vis=050911H30a_rx0.usb.tsys.bp.bl \
        axis=time,ampl device=/xs nxy=1,1"
sleep 5
smablflag vis=050911H30a_rx0.usb.tsys.bp.bl \
        axis=time,ampl device=/xs 
continue:
echo " "
echo "Plot overall baseline-based spectra"
echo " " 
echo "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"
sleep 5 
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 " "
echo "For the H30a  line in spectral chunks 19,20,21"
echo " "
echo "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"
sleep 3
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
echo " "
echo "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"
sleep 3
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
sleep 10 
echo " " 
echo "For the CO2-1 line in spectral chunks 3,4,5"
echo " "
echo "smauvspec vis=050911H30a_rx0.usb.tsys.bp.bl \
        select='source(n7538),window(3,4,5)' \
        interval=1000  hann=7 yrange=0,6 \
        axis=freq,both device=/xs nxy=3,5"
sleep 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
sleep 10
#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 " "
echo "uvflag edge channels"
echo " "
echo "uvflag vis=050911H30a_rx0.usb.tsys.bp.bl edge=3,3 \
        flagval=flag"
sleep 10 
uvflag vis=050911H30a_rx0.usb.tsys.bp.bl edge=3,3 \
	flagval=flag
echo " "
echo "Check spectrum  on one of the baselines"
echo " "
echo "smauvspec vis=050911H30a_rx0.usb.tsys.bp.bl \
        select='source(n7538),ant(1)(4)' \
        interval=1000  \
        axis=freq,both device=/xs nxy=1,1"
sleep 5 
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 
sleep 5
echo " "
echo "5. ANTENNA GAIN CALIBRATION"
echo " "
echo "Up to this point, we have corrected for bandpass errors."
echo " "
echo "Check and flag visibility "
echo " "
echo "smablflag vis=050911H30a_rx0.usb.tsys.bp.bl \
         select='source(n7538,2202*,0102*)' \
        device=/xs options=nocal axis=time,amp"
echo " "
# smablflag vis=050911H30a_rx0.usb.tsys.bp.bl \
#        select='source(n7538,2202*,0102*)' \
#        device=/xs options=nocal axis=time,amp
echo "or"
echo "smauvplt  vis=050911H30a_rx0.usb.tsys.bp.bl"
echo "select='source(n7538,2202*,0102*)'"
echo "device=/xs options=nocal axis=time,amp nxy=3,5"
sleep 15
smauvplt  vis=050911H30a_rx0.usb.tsys.bp.bl \
	select='source(n7538,2202*,0102*)' \
	device=/xs options=nocal axis=time,amp nxy=3,5
echo " "
echo "Split the gain calibrators from the multiple source file"
echo " "
echo "uvsplit vis=050911H30a_rx0.usb.tsys.bp.bl "
echo "select='source(2202+422,0102+584)' options=nowindow"
sleep 10 
\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 " "
echo "rename the calibrator files"
echo " "
echo "\rm -r 2202+422.tsys.bp.bl 0102+584.tsys.bp.bl"
echo "mv 2202+422.230409 2202+422.tsys.bp.bl"
echo "mv 0102+584.230409 0102+584.tsys.bp.bl"
sleep 5 
\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 " "
echo "Solve for antenna gains"
echo " "
echo "For 0102+584"
echo "smamfcal vis=0102+584.tsys.bp.bl  \
select='source(0102+584)'  \
refant=3 interval=2 options=nopassol"
sleep 3
smamfcal vis=0102+584.tsys.bp.bl  \
select='source(0102+584)'  \
refant=3 interval=2 options=nopassol
sleep 5
echo " "
echo "Plot the gain solutions"
echo " "
echo "smagpplt vis=0102+584.tsys.bp.bl device=/xs \
        yaxis=amp,phase log=0102+584gain.log \
        options=gains  polyfit=6 \
        nxy=2,4"
sleep 5 
smagpplt vis=0102+584.tsys.bp.bl device=/xs \
        yaxis=amp,phase log=0102+584gain.log \
        options=gains  polyfit=6 \
        nxy=2,4
sleep 10 
echo " "
echo "For 2202+422"
echo " "
echo "smamfcal vis=2202+422.tsys.bp.bl  \
select='source(2202+422)'  \
refant=3 interval=2 options=nopassol"
sleep 5 
smamfcal vis=2202+422.tsys.bp.bl  \
select='source(2202+422)'  \
refant=3 interval=2 options=nopassol
sleep 5
echo " "
echo "Plot the gain solutions"
echo " "
echo "smagpplt vis=2202+422.tsys.bp.bl device=/xs \
        yaxis=amp,pha  \
        options=gains  polyfit=6 \
        nxy=2,4"
sleep 5 
smagpplt vis=2202+422.tsys.bp.bl device=/xs \
        yaxis=amp,pha  \
        options=gains  polyfit=6 \
        nxy=2,4
sleep 10
echo " "
echo "Merge the two gain tables"
echo " "
echo "smagpplt vis=2202+422.tsys.bp.bl,0102+584.tsys.bp.bl device=/xs \
        yaxis=amp,pha  \
        options=gains,merge polyfit=8 \
        nxy=1,6"
sleep 10 
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
sleep 3
echo " "
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 " "
echo "Gpcopy the merged gain table back"
echo " "
echo "gpcopy vis=0102+584.tsys.bp.bl out=050911H30a_rx0.usb.tsys.bp.bl \
     mode=copy options=nopass"
sleep 5 
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
sleep 5
echo " "
echo "6. BOOTSTRAPE THE FLUX DENSITY SCALE"
echo " "
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."
echo " "
sleep 10
echo "smaflux vis=050911H30a_rx0.usb.tsys.bp.bl"
echo "      select='source(ceres)' mirhome=$MIR"
echo " "
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.
#  
sleep 10
echo " "
echo "Apply the gains to the data"
echo " "
echo "uvaver vis=050911H30a_rx0.usb.tsys.bp.bl out=050911H30a_rx0.usb.tsys.bp.bl.av"
sleep 10  
\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 " "
echo "Check for spectra of NGC7538"
echo " "
echo "for overall"
echo "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 "
sleep 5 
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
sleep 10
echo " "
echo "For H30a line"
echo " "
echo "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"
sleep 5 
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
sleep 10 
echo " "
echo "Uvsplit n7538" 
echo " "
echo "uvsplit vis=050911H30a_rx0.usb.tsys.bp.bl.av select='source(n7538)' " 
echo "options=nowindow"
sleep 5 
\rm -r n7538.230409 n7538.h30a
uvsplit vis=050911H30a_rx0.usb.tsys.bp.bl.av select='source(n7538)' \
	options=nowindow
echo " "
echo "Check up the rest frequency"
echo " "
echo "prthd in=n7538.230409"
sleep 5 
prthd in=n7538.230409
echo " "
echo "Rename the output file"
echo " "
mv n7538.230409 n7538.h30a 
sleep 3
echo "skip the following step for identification of line transitions "
echo "smauvspec vis=n7538.h30a select='ant(3)(4)' \
        interval=1000 catpath='$MIRCAT/jplcat' options=jplcat \
        axis=freq,both device=/xs nxy=1,1"
sleep 10
echo " "
echo "for CO2-1, using the following steps to reset the rest frequency"
echo "mv n7538.230409 n7538.co21"
echo "puthd in=n7538.co21/restfreq value=230.538 type=double"
echo " "
echo " "
echo "7. CONTINUUM IMAGE"
echo " "
\rm -r n7583con.map n7583con.beam n7583con.icmp
sleep 10 
echo " "
echo "Make a dirty image"
echo "using the channel data between CO line and H30a line,"
echo "i.e., from spectral chunks 6-18"
echo " "
echo "invert vis=n7538.h30a map=n7583con.map beam=n7583con.beam \
        imsize=256,256 cell=0.2 sup=0 \
        line=channel,1,641,1664"
sleep 5 
invert vis=n7538.h30a map=n7583con.map beam=n7583con.beam \
        imsize=256,256 cell=0.2 sup=0 \
        line=channel,1,641,1664 
echo " "
echo "Display the dirty image"
echo " "
echo "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"
sleep 5 
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 " "
echo "Clean  the dirty image"
echo " "
sleep 3
echo "clean map=n7583con.map beam=n7583con.beam out=n7583con.icmp gain=0.08 \
        cutoff=0 niters=3500 \
        region='boxes(110,110,146,146)'"
clean map=n7583con.map beam=n7583con.beam out=n7583con.icmp gain=0.08 \
        cutoff=0 niters=3500 \
        region='boxes(110,110,146,146)'
sleep 10
echo " "
echo "Restore cleaned image"
echo " "
echo "restor model=n7583con.icmp beam=n7583con.beam map=n7583con.map  \
        out=n7583con.icln"
sleep 5 
\rm -r n7583con.icln
restor model=n7583con.icmp beam=n7583con.beam map=n7583con.map  \
        out=n7583con.icln
sleep 5
echo " "
echo "Display the restored clean image"
echo " "
echo "cgdisp in=n7583con.icln type=pixel region=arcsec,'boxes(-15,-15,15,15)' \
        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"
sleep 5 
cgdisp in=n7583con.icln type=pixel region=arcsec,'boxes(-15,-15,15,15)' \
        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
#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.
echo " "
echo "8. LINE/CONTINUUM SEPARATION"
echo " "
echo "Plot H30a spectrum as function of velocity"
echo " "
echo "smauvspec vis=n7538.h30a "
echo "select='source(n7538),window(19,20,21)' "
echo "interval=10000  options=nobase,avall "
echo "axis=vel,both device=/xs nxy=1,1 yrange=0,10"
sleep 10 
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
sleep 10
echo " "
echo "Convert sky-frequency to velocity -220 to 100 at"
echo "resolution of 1 km/s"
echo " "
echo "uvaver vis=n7538.h30a line=velocity,320,-220,1 out=n7538.h30a.vel"
sleep 10 
\rm -r n7538.h30a.vel
uvaver vis=n7538.h30a line=velocity,320,-220,1 out=n7538.h30a.vel
echo " "
echo "Plot the velocity spectrum"
echo " "
echo "smauvspec vis=n7538.h30a.vel \
        interval=10000  options=nobase,avall \
        axis=chan,both device=/xs nxy=1,1 yrange=0,10"
sleep 10 
smauvspec vis=n7538.h30a.vel \
        interval=10000  options=nobase,avall \
        axis=chan,both device=/xs nxy=1,1 yrange=0,10
sleep 10
echo " "
echo "Take out the continuum using linear fitting to" 
echo "the line free channels"
echo " "
echo "uvlin vis=n7538.h30a.vel chans=1,30,250,280 out=n7538.h30a.vel.lin \
        order=1 mode=line"
sleep 10 
\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
sleep 3
echo " "
echo "Plot the continuum-free spectrum"
echo " "
echo "smauvspec vis=n7538.h30a.vel.lin options=nobase,avall  \
        interval=1000 \
        axis=chan,both device=/xs nxy=1,1"
sleep 10 
smauvspec vis=n7538.h30a.vel.lin options=nobase,avall  \
        interval=1000 \
        axis=vel,both device=/xs nxy=1,1
sleep 10
echo " "
echo "9. CONSTRUCT LINE IMAGE CUBE" 
echo " "
sleep 5
echo "Make a dirty image cube"
echo " "
echo "invert vis=n7538.h30a.vel.lin map=n7538.h30alin.map beam=n7538.h30alin.beam \
        imsize=256,256 cell=0.2 sup=0 \
        line=channel,25,100,5"
sleep 15 
\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,25,100,5
echo " "
echo "Clean the dirty image cube"
echo " "
echo "clean map=n7538.h30alin.map beam=n7538.h30alin.beam out=n7538.h30alin.icmp gain=0.08 \
        cutoff=0 niters=2500 \
        region='boxes(110,110,146,146)(1,25)'"
sleep 5 
clean map=n7538.h30alin.map beam=n7538.h30alin.beam out=n7538.h30alin.icmp gain=0.08 \
        cutoff=0 niters=2500 \
	region='boxes(110,110,146,146)(1,25)'
echo " "
echo "Restore the cleaned image cube"
echo " "
echo "restor model=n7538.h30alin.icmp beam=n7538.h30alin.beam map=n7538.h30alin.map  \
        out=n7538.h30alin.icln"
sleep 10 
\rm -r n7538.h30alin.icln
restor model=n7538.h30alin.icmp beam=n7538.h30alin.beam map=n7538.h30alin.map  \
        out=n7538.h30alin.icln
echo " "
echo "Display the restired clean image cube"
echo " "
echo "cgdisp in=n7538.h30alin.icln type=pixel region=arcsec,'boxes(-5,-5,5,5)(1,25)' \
        xybin=1,1 device=/xs nxy=5,5 options=full,beambr,wedge,trlab,3val \
        labtyp=arcsec,arcsec range=0,2.5,lin,2 cols1=7 csize=0.5,0.5,0.5"
sleep 5 
cgdisp in=n7538.h30alin.icln type=c region=arcsec,'boxes(-5,-5,5,5)(1,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.1 levs1=-5,3,4,8,16,32,64,128,256 
sleep 10 
echo " "
echo "statistics on the channel images"
echo " "
echo "On source region"
echo " "
echo "imstat in=n7538.h30alin.icln region=quarter \
        axes=ra,dec"
sleep 10 
imstat in=n7538.h30alin.icln region=quarter \
        axes=ra,dec
echo " "
echo "Off source region"
echo " "
echo "imstat in=n7538.h30alin.icln region='boxes(30,30,60,60)(1,25)' \
        axes=ra,dec"
sleep 10  
imstat in=n7538.h30alin.icln region='boxes(30,30,60,60)(1,25)' \
        axes=ra,dec
echo "Display the restired clean image cube"
echo " "
echo "cgdisp in=n7538.h30alin.icln type=pixel region=arcsec,'boxes(-5,-5,5,5)(1,25)' \
        xybin=1,1 device=/xs nxy=5,5 options=full,beambr,wedge,trlab,3val \
        labtyp=arcsec,arcsec range=0,2.5,lin,2 cols1=7 csize=0.5,0.5,0.5"
sleep 5
cgdisp in=n7538.h30alin.icln type=pixel region=arcsec,'boxes(-5,-5,5,5)(1,25)' \
        xybin=1,1 device=/xs nxy=5,5 options=full,beambr,wedge,trlab,3val \
        labtyp=arcsec,arcsec range=0,4.0,lin,2 cols1=7 csize=0.5,0.5,0.5
echo " "
echo "The end"
# One can continue on for making a high velocity resolution line cube
# For position-velocity and moment analysis, please reference ATNF's
# users guide http://www.atnf.csiro.au/computing/software/miriad/userhtml.html
echo " "
exit
