A Step-by-Step Procedure of A Step-by-Step Procedure of The Reduction of Spectral Line Data from The ASIC (2GHz) Correlator (Demonstration)

Jun-Hui Zhao (SAO)

2005-11-01



1. INTRODUCTION

A step-by-step procedure of the reduction for the SMA spectral line data in Miriad was presented on one of the CF's Linux computers at CfA in real time based on a recent SMA data set with a size of 1 Gbytes. The SMA data was taken from the observations of the HII region NGC 7538 for the hydrogen recombination line (H30a) at 231.901 GHz. The observation was made on 2005-9-11. In the same side band, the CO (2-1) line is also included. The procedure contains the detailed setups for the Miriad tasks used in the data reduction from loading data to constructing a spectral line image cube. The flow chart below (Section 2) summarizes the nine major steps listed in the first column:
1. Read MIR Data,  
2. Inspect Data,
3. Correct Amplitude,
4. Calibrate Bandpass, 
5. Gain Calibration, 
6. Flux Density Scale, 
7. Continuum Imaging,
8. Continuum Subtraction,
9. Construct Line Cube. 
The second column lists the Miriad tasks used in the process. The third column gives the uv datasets, gain/bandpass tables, image files, and plots which are produced in the course of the data reduction. The realtime demonstration has been scripted following the steps listed in the flow chart. The C-shell script used in the demonstration is shown in Section 3 and can be downloaded. Running the script will repeat the entire data processing shown in the demo.

Courtesy: The SMA data used in this demo can be obtained from the SMA archival. please send RTDC a request for uploading the file 050911_04:21:35 for you.

2. FLOW CHART

Here is the postscript version, which can be downloaded:
fchart1.ps
fchart2.ps
fchart3.ps






3. THE SCRIPT
#!/bin/csh -f
# 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


Updated: 2005-nov-07
Updated: 2008-nov-05 (with L. Zhu's suggestion)
Updated: 2015-jul-24 (for the new Miriad cover page)

Questions: smamiriad@cfa.harvard.edu