#############################################################
#
# This imaging script starts with fully calibrated SMA data (see "Sample MIR Calibration Script" in "Searching and Requesting Data: Overview").
# It uses the MIRIAD software package to produce the following data products:
#
# 1) A spectral line image cube for each sideband, FITS format.  
#    If the experiment has multiple spectral resolutions, there will be multiple image cubes per sideband.
# 2) An integrated-intensity (pseudo-continuum), single channel image, FITS format.
# 3) Representive image plots, which are output in postscript format but made available to the user in jpg format.
#
#############################################################


#! /bin/csh


#------------------------------------------------------------
#Run correct version of miriad:
source /opt/miriad/miriad_cvs/MIRRC.linux


#------------------------------------------------------------
#Get Rest frequency and velocity resolution(s): 
uvlist vis='070913_125512_B5IRS1_lsb' options=spectra,full
uvlist vis='070913_125512_B5IRS1_usb' options=spectra,full


#------------------------------------------------------------
#Modify header.  
puthd in='070913_125512_B5IRS1_lsb'/restfreq value=230.53799
puthd in='070913_125512_B5IRS1_usb'/restfreq value=230.53799


#------------------------------------------------------------
#Flag edge channels and check.  The value of "edge" is taken to be 1/16 of the total number of channels per band.
#If there are multiple resolutions, there will be multiple values of "edge."

uvflag vis='070913_125512_B5IRS1_lsb' flagval=f edge=8
smauvspec vis='070913_125512_B5IRS1_lsb' device=/xs interval=1E21 axis=vel,amp options=avall
smauvspec vis='070913_125512_B5IRS1_lsb' device=/xs interval=1E21 axis=vel,amp options=avall, select='window(1,2,3)'

uvflag vis='070913_125512_B5IRS1_usb' flagval=f edge=8
smauvspec vis='070913_125512_B5IRS1_usb' device=/xs interval=1E21 axis=vel,amp options=avall
smauvspec vis='070913_125512_B5IRS1_usb' device=/xs interval=1E21 axis=vel,amp options=avall, select='window(1,2,3)'


#Image (line = velocity, [#chan], [start (km/s)], [width (km/s)], [step (km/s)]):

#LINE: Image, clean down to 3 sigma (theoretical), restore, and get statistics for line.  Run first "invert" below to get 1 sigma value.

invert vis='070913_125512_B5IRS1_lsb' map='070913_125512_B5IRS1_lsb.map' beam='070913_125512_B5IRS1_lsb.beam' options=double,systemp slop=0.5 sup=0 line=vel,2423,11572.475,1.057,1.057
clean map='070913_125512_B5IRS1_lsb.map' beam='070913_125512_B5IRS1_lsb.beam' out='070913_125512_B5IRS1_lsb.cln' niters=10000 cutoff=0.11811
restor model='070913_125512_B5IRS1_lsb.cln' map='070913_125512_B5IRS1_lsb.map' beam='070913_125512_B5IRS1_lsb.beam' out='070913_125512_B5IRS1_lsb.cln.restor'

invert vis='070913_125512_B5IRS1_usb' map='070913_125512_B5IRS1_usb.map' beam='070913_125512_B5IRS1_usb.beam' options=double,systemp slop=0.5 sup=0 line=vel,2423,-1430.665,1.057,1.057
clean map='070913_125512_B5IRS1_usb.map' beam='070913_125512_B5IRS1_usb.beam' out='070913_125512_B5IRS1_usb.cln' niters=10000 cutoff=0.12555
restor model='070913_125512_B5IRS1_usb.cln' map='070913_125512_B5IRS1_usb.map' beam='070913_125512_B5IRS1_usb.beam' out='070913_125512_B5IRS1_usb.cln.restor'


#CONT: Image, clean, restore, and get statistics for continuum, both sidebands.  Run "invert" to get 1 sigma value.

invert vis='070913_125512_B5IRS1_usb,070913_125512_B5IRS1_lsb' map='070913_125512_B5IRS1_both.mapc' beam='070913_125512_B5IRS1_both.beamc' options=double,systemp slop=0.5 sup=0 line=vel,1,-1430.665,5123,5123
clean map='070913_125512_B5IRS1_both.mapc' beam='070913_125512_B5IRS1_both.beamc' out='070913_125512_B5IRS1_both.clnc' niters=10000 cutoff=0.0018036
restor model='070913_125512_B5IRS1_both.clnc' map='070913_125512_B5IRS1_both.mapc' beam='070913_125512_B5IRS1_both.beamc' out='070913_125512_B5IRS1_both.cln.restorc'


imstat in='070913_125512_B5IRS1_lsb.cln.restor'
# Total                  Sum      Mean      rms     Maximum   Minimum    Npoints
#                     -1.646E+03-6.237E-04 4.368E-02  4.27    -0.695     2638647
#####imstat in='070913_125512_B5IRS1_lsb.cln.restor' region='images(x,y)' --> z

imstat in='070913_125512_B5IRS1_usb.cln.restor'
# Total                  Sum      Mean      rms     Maximum   Minimum    Npoints
#                      3.667E+03 1.390E-03 0.106      13.2     -3.95     2638647
#####imstat in='070913_125512_B5IRS1_usb.cln.restor' region='images(1000,1500)' --> 1365 (xxx km/s)

imstat in='070913_125512_B5IRS1_both.cln.restorc'
# Total                  Sum      Mean      rms     Maximum   Minimum    Npoints
#                      0.759     6.967E-04 7.322E-03 6.411E-02-1.618E-02    1089


#------------------------------------------------------------
#Convert to fits format:
fits in='070913_125512_B5IRS1_lsb.cln.restor' op=xyout out='070913_125512_B5IRS1_lsb.cln.restor.fits'
fits in='070913_125512_B5IRS1_usb.cln.restor' op=xyout out='070913_125512_B5IRS1_usb.cln.restor.fits'
fits in='070913_125512_B5IRS1_both.cln.restorc' op=xyout out='070913_125512_B5IRS1_both.cln.restorc.fits'


#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# 1) Two Spectra
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

uvspec vis='070913_125512_B5IRS1_lsb,070913_125512_B5IRS1_usb' device=/ps interval=1E21 axis=freq,amp options=avall


#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# 2) Psuedo-continuum.  Add everything up
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#..........Make zeroeth moment map (integrated-intensity)
moment in='070913_125512_B5IRS1_both.cln.restorc' out='070913_125512_B5IRS1_both.cln.restorc.mom0' mom=0

#..........Make color map w/ contours and beam
cgdisp in='070913_125512_B5IRS1_both.cln.restorc.mom0,070913_125512_B5IRS1_both.cln.restorc.mom0' device=/cps type=p,c slev=p,10 levs1=1,2,3,4,5,6,7,8,9 range=0,0,0,8 beamtyp=b,l,3 options=blacklab,full,wedge labtyp=hms,dms cols1=0


#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
# 3) Channel map at frequency of peak emission
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#..........Map w/ contours
cgdisp in='070913_125512_B5IRS1_usb.cln.restor,070913_125512_B5IRS1_usb.cln.restor' device=/cps type=p,c slev=p,10 levs1=1,2,3,4,5,6,7,8,9 range=0,0,0,8 beamtyp=b,l,3 options=blacklab,full,wedge labtyp=hms,dms cols1=0 region='images(1365,1365)'