#! /bin/bash
# JPM - 12/31/99 - initial
# JPM - 01/04/00 - add cases
#
# script assumes that you have the miriad environment enabled and that
# there are three data directories in the cwd:
# 0609-157, 3c273, N2023MM1
case $# in 
0) 
	echo 'Usage: test command';
	echo '   listobs  - write file of observations and correlator setup';
	echo '   linecal  - linecalibrate all observations';
	echo '   flag     - flag all observations';
	echo '   fluxcal  - bootstrap flux for phase cal';
	echo '   calib    - do p/a calibration - apply to source and phasecal';
	echo '   mappha   - map phase calibrator';
	echo '   mapsouc  - map source - continuum';
	echo '   mapsouda - map source deutero-ammonia';
	echo '   full     - do full reduction (all of above)';
	echo '   cleanup  - remove all reduction files (start over!)';
	echo '   cleansou - remove all source images';

esac;

case "$1" in 

listobs)
	# preliminaries
	listobs in=N2023MM1,3c273,0609-157 log=n2023.listobs
	corset freq=86.230 iffreq=-150 device=corset.ps/ps cormode=4 \
	    corf=296.24634,311.25488,465.00244,605.75562 corbw=25,25,25,25
;;

linecal)
	# phasem variable stores measurements of line lengths - compute 
	# corrections over course of observations
  	linecal vis=0609-157
  	linecal vis=3c273
  	linecal vis=N2023MM1

	# linecal writes a gains table with the corrections - 
        # uvcat applies those corrections
  	uvcat vis=0609-157 out=0609_ll
  	uvcat vis=3c273 out=3c273_ll
  	uvcat vis=N2023MM1 out=n2023_ll # this crashes my laptop!
	;;

flag)
	# flag data edge, shadowed and spike data
  	uvflag vis=0609_ll edge=2,2 flagval=flag
  	uvflag vis=3c273_ll edge=2,2 flagval=flag
  	uvflag vis=n2023_ll edge=2,2 flagval=flag
  	uvflag vis=0609_ll 'select=shadow(6.1)' flagval=flag
  	uvflag vis=3c273_ll 'select=shadow(6.1)' flagval=flag
  	uvflag vis=n2023_ll 'select=shadow(6.1)' flagval=flag
	# any clipping? - use uvplt,uvspec(t) to view
  	uvflag vis=0609_ll 'select=-amp(0,55)' flagval=flag
  	uvflag vis=3c273_ll 'select=-amp(0,20)' flagval=flag
  	uvflag vis=n2023_ll 'select=-amp(0,130)' flagval=flag
	# reconstruct the wideband channels
  	uvwide vis=0609_ll out=0609.nw
  	uvwide vis=3c273_ll out=3c273.nw
  	uvwide vis=n2023_ll out=n2023.nw
;;

fluxcal)
	# flux calibration: selfcal phase and flux calib then use bootflux 
	# to bootstrap the flux of the phase/amp calibrator
	# use refant=7 since it's closest antenna to array center
  	selfcal vis=0609.nw interval=10 refant=7 line=wide,1,1
  	selfcal vis=3c273.nw interval=10 refant=7 line=wide,1,1
	# Note: it's not finding the flux table so must specify flux of 3c273
	# flux from website is 10.8 Jy (1.4 Jy rms) from 01Oct99
	# use primary=3c273,6.9
  	bootflux vis=0609_ll,3c273_ll line=wide,1,1 taver=10 \
		primary=3c273,10.8 log=bootflux.log
	# doesn't look too good
	# call it 8.8 from el>~30
;;

calib)
	# phase and amplitude calibration
	# options=amp(do amp and phase self-cal),apriori(no input model-
	# causes data to be scaled by the known flux of the source),noscale; 
  	gmake vis=0609.nw interval=1 options=amp,apriori,noscale refant=7 \
		fluxes=0609-157,8.8 line=wide,2,1,1,1 out=0609.gvis
  	gfiddle vis=0609.gvis device=/xs out=0609.gfit fit=poly
	# apply fit to calibrator to see how well we've done
	# map calibrator
   	gapply vis=0609.nw gvis=0609.gfit out=0609.cal
  	gapply vis=n2023.nw gvis=0609.gfit out=n2023.cal
;;

mappha)
   	invert vis=0609.cal map=0609.mp beam=0609.bm imsize=512 cell=2 \
		line=wide,1,2,1,1 options=systemp sup=0
  	clean map=0609.mp beam=0609.bm out=0609.cln
  	restor model=0609.cln beam=0609.bm map=0609.mp out=0609.res
  	cgdisp in=0609.res,0609.res 'region=b(129,129,394,394)' \
		type=grey,contour device=/xs
;;

mapsouc)
  	invert vis=n2023.cal map=n2023.mp beam=n2023.bm imsize=512 cell=2 \
		sup=0 options=systemp line=wide,1,2,1,1
  	clean map=n2023.mp beam=n2023.bm out=n2023.cln
  	restor model=n2023.cln beam=n2023.bm map=n2023.mp out=n2023.res
	cgdisp in=n2023.res,n2023.res type=grey,contour device=/xw \
		'region=b(129,129,394,394)'
;;

mapsouda)
	invert vis=n2023.cal map=nh2d.mp beam=nh2d.bm imsize=512 cell=2 \
		sup=0 'select=window(3)' line=velocity,20,2.6,0.68,0.68
#	contsub in=nh2d.mp out=nh2d_contsub.mp contchan='(15,20)' mode=mean
# 	this fails with: ### Fatal Error:  Bad dimension ???
#	don't worry for now (other ways to do this) - only a 5% effect
        clean map=nh2d.mp beam=nh2d.bm out=nh2d.cln \
		'region=b(225,240,280,280)'
        restor model=nh2d.cln beam=nh2d.bm map=nh2d.mp out=nh2d.res
        cgdisp in=nh2d.res,nh2d.res type=grey,contour device=/xw \
        'region=b(129,129,394,394)' nxy=2,2 options=beambl,3value\
	 labtyp=hms,dms, slev=a,.25
;;

mapsouso)
        invert vis=n2023.cal map=so.mp beam=so.bm imsize=512 cell=2 \
                sup=0 'select=window(1)' line=velocity,20,2.6,0.68,0.68
#       contsub in=so.mp out=so_contsub.mp contchan='(15,20)' mode=mean
#       this fails with: ### Fatal Error:  Bad dimension ???
#       don't worry for now (other ways to do this) - only a 5% effect
        clean map=so.mp beam=so.bm out=so.cln \
                'region=b(225,240,280,280)'
        restor model=so.cln beam=so.bm map=so.mp out=so.res
        cgdisp in=so.res,so.res type=grey,contour device=/xw \
        'region=b(129,129,394,394)' nxy=2,2 options=beambl,3value\
         labtyp=hms,dms, slev=a,.25
;;

cleanup)
	# remove everything that's there
	if [ -e "0609_ll" ]; then
		rm -r 0609_ll
	fi
	if [ -e "0609.nw" ]; then
		rm -r 0609.nw
	fi
	if [ -e "0609.bm" ]; then
		rm -r 0609.bm
	fi
	if [ -e "0609.mp" ]; then
		rm -r 0609.mp
	fi
	if [ -e "0609.gvis" ]; then
		rm -r 0609.gvis
	fi
	if [ -e "0609.gfit" ]; then
		rm -r 0609.gfit
	fi
	if [ -e "3c273_ll" ]; then
		rm -r 3c273_ll
	fi
	if [ -e "3c273.nw" ]; then
		rm -r 3c273.nw
	fi
;;

cleansou) 
        if [ -e "n2023.mp" ]; then
                rm -r n2023.mp
        fi
        if [ -e "n2023.bm" ]; then
                rm -r n2023.bm
        fi
        if [ -e "n2023.cln" ]; then
                rm -r n2023.cln
	fi
  	if [ -e "n2023.res" ]; then
		rm -r n2023.res
	fi
	if [ -e "nh2d.res" ]; then
		rm -r nh2d.res
	fi
	if [ -e "nh2d.mp" ]; then
		rm -r nh2d.mp
	fi
	if [ -e "nh2d.cln" ]; then
		rm -r nh2d.cln
	fi
	if [ -e "nh2d.bm" ]; then
		rm -r nh2d.bm
	fi
;;

full)
	# preliminaries
	listobs in=N2023MM1,3c273,0609-157 log=n2023.listobs
	corset freq=86.230 iffreq=-150 device=corset.ps/ps cormode=4 \
		corf=296.24634,311.25488,465.00244,605.75562 corbw=25,25,25,25
	# phasem variable stores measurements of line lengths - compute 
	# corrections over course of observations
  	linecal vis=0609-157
  	linecal vis=3c273
  	linecal vis=N2023MM1
	# linecal writes a gains table with the corrections - 
        # uvcat applies those corrections
  	uvcat vis=0609-157 out=0609_ll
  	uvcat vis=3c273 out=3c273_ll
  	uvcat vis=N2023MM1 out=n2023_ll # this crashes my laptop!
	# flag data edge, shadowed and spike data
  	uvflag vis=0609_ll edge=2,2 flagval=flag
  	uvflag vis=3c273_ll edge=2,2 flagval=flag
  	uvflag vis=n2023_ll edge=2,2 flagval=flag
  	uvflag vis=0609_ll 'select=shadow(6.1)' flagval=flag
  	uvflag vis=3c273_ll 'select=shadow(6.1)' flagval=flag
  	uvflag vis=n2023_ll 'select=shadow(6.1)' flagval=flag
	# any clipping? - use uvplt,uvspec(t) to view
  	uvflag vis=0609_ll 'select=-amp(0,55)' flagval=flag
  	uvflag vis=3c273_ll 'select=-amp(0,20)' flagval=flag
  	uvflag vis=n2023_ll 'select=-amp(0,130)' flagval=flag
	# reconstruct the wideband channels
  	uvwide vis=0609_ll out=0609.nw
  	uvwide vis=3c273_ll out=3c273.nw
  	uvwide vis=n2023_ll out=n2023.nw
	# flux calibration: selfcal phase and flux calib then use bootflux 
	# to bootstrap the flux of the phase/amp calibrator
	# use refant=7 since it's closest antenna to array center
  	selfcal vis=0609.nw interval=10 refant=7 line=wide,1,1
  	selfcal vis=3c273.nw interval=10 refant=7 line=wide,1,1
	# Note: it's not finding the flux table so must specify flux of 3c273
	# flux from website is 10.8 Jy (1.4 Jy rms) from 01Oct99
	# use primary=3c273,6.9
  	bootflux vis=0609_ll,3c273_ll line=wide,1,1 taver=10 \
		primary=3c273,10.8 log=bootflux.log
	# doesn't look too good
	# call it 8.8 from el>~30
	# phase and amplitude calibration
	# options=amp(do amp and phase self-cal),apriori(no input model-
	# causes data to be scaled by the known flux of the source),noscale; 
  	gmake vis=0609.nw interval=1 options=amp,apriori,noscale refant=7 \
		fluxes=0609-157,8.8 line=wide,2,1,1,1 out=0609.gvis
  	gfiddle vis=0609.gvis device=/xs out=0609.gfit fit=poly
	# apply fit to calibrator to see how well we've done
	# map calibrator
   	gapply vis=0609.nw gvis=0609.gfit out=0609.cal
  	gapply vis=n2023.nw gvis=0609.gfit out=n2023.cal
   	invert vis=0609.cal map=0609.mp beam=0609.bm imsize=512 cell=2 \
		line=wide,1,2,1,1 options=systemp sup=0
  	clean map=0609.mp beam=0609.bm out=0609.cln
  	restor model=0609.cln beam=0609.bm map=0609.mp out=0609.res
  	cgdisp in=0609.res,0609.res 'region=b(129,129,394,394)' \
		type=grey,contour device=/xs
  	invert vis=n2023.cal map=n2023.mp beam=n2023.bm imsize=512 cell=2 \
		sup=0 options=systemp line=wide,1,2,1,1
  	clean map=n2023.mp beam=n2023.bm out=n2023.cln
  	restor model=n2023.cln beam=n2023.bm map=n2023.mp out=n2023.res
	cgdisp in=n2023.res,n2023.res type=grey,contour device=/xw \
		'region=b(129,129,394,394)'
	invert vis=n2023.cal map=n2023.mp beam=n2023.bm imsize=512 cell=2 \
		sup=0 options=systemp line=wide,1,2,1,1
        clean map=n2023.mp beam=n2023.bm out=n2023.cln
        restor model=n2023.cln beam=n2023.bm map=n2023.mp out=n2023.res
        cgdisp in=n2023.res,n2023.res type=grey,contour device=/xw \
                'region=b(129,129,394,394)'
;;

esac;