#! /bin/csh -f

#######################
# warm up with MIRIAD #
#######################

#admin: getting the data

# set tar=../SgrA.mir.tar.gz
# tar -zxf $tar
# # set vis=$tar:t:r:r
# set vis=SgrA.mir

# get a listing of what happened

# listobs vis=$vis log=listobs.log
listobs vis=SgrA.mir log=listobs.log

# define sources S (science source), B (passband calibrator) and G (phase/gain calibrator) as in listobs.log
# set FILE=SgrA.mir
set sname=GCRING
set bname=3C345
set gname=1733-130

# you can list stuff like

# uvlist vis=$FILE options=spectra

# set an antpos file for baseline correction. Look this up in listobs.log
# two are available for SgrA.mir: 090225 and 090331
# ask Dan how we know this?

# set antpos=$MIRCAT/baselines/carma/antpos.090331

# reference antennas

set bref=2
set pref=2

# work starts here------------------------------------------------------------------------------------------

# select only useful data by removing NOISE and auto-correlations

rm -rf all.vis
uvcat vis=SgrA.mir select='-source(NOISE),-auto' out=all.vis

########################
# look at data quality #
########################

# 1 check phase vs. time; there should be a smooth trend
# smauvplt vis=all.vis device=/xs axis=time,phase options=nocal,nopass
# leave out phase(gain) and bandpass calibration with nocal,nopass
# notice that ant1 has some issue remove it by flagging

uvflag vis=all.vis flagval=flag select="ant(1)"

# 2 check amp vs. time; it should be constant and the same for all antennas
# smauvplt vis=all.vis device=/xs axis=time,amp options=nocal,nopass
# option line=$WIDE
# ask Dan what this does exactly?

# 3 check source phase vs. time;
# smauvplt vis=all.vis device=/xs line=$WIDE axis=time,phase

# not sure if this actually refers to source, why can I not use GCRING?

##########################
# ready for calibrations #
##########################

#1 baseline calibration

rm -rf all_1.vis
uvedit vis=all.vis out=all_1.vis apfile=$MIRCAT/baselines/carma/antpos.090331

#2 linelength calibration

rm -rf all_2.vis
linecal vis=all_1.vis
uvcat vis=all_1.vis out=all_2.vis
# smauvplt vis=all_2.vis device=/xs axis=time,phase options=nocal,nopass

#3 bandpass calibration

# use a short observation of a bright quasar 3C345 to calibrate the data
# such that the bandpass has zero phase and constant amplitude

rm -rf all_3.vis
#mfcal vis=all_2.vis refant=2 interval=999 select="source(3C345)"
mfcal vis=all_2.vis refant=$bref interval=999 select="source($bname)"
# gpplt vis=all_2.vis yaxis=phase nxy=5.3 device=/xs options=bandpass
# gpplt vis=all_2.vis yaxis=amp nxy=5.3 device=/xs options=bandpass
uvcat vis=all_2.vis out=all_3.vis options=nocal

#4 amplitude/phase(gain)calibration
# mselfcal vis=all_3.vis refant=2 select="source(1733-130)" options=amp,apriori,noscale interval=5
mselfcal vis=all_3.vis refant=$pref select="source($gname)" options=amp,apriori,noscale interval=5
# widen the validity interval(days) to something big enough to safely interpolate intop the source data

puthd in=all_3.vis/interval value=0.1
rm -rf all_4.vis

# apply gains

uvcat vis=all_3.vis out=all_4.vis

#5 phase-only calibration

# above we calibrated both phase and amplitude, another option is to calibrate phase only and use a 
# planet for amplitude scaling

###################################
# sanity check after calibrations #  probably need to check this
###################################

# 1 

# after linelength calibration plot antenna gain, phase gain should be smooth,amp should be about 1
# gpplt vis=all_1.vis yaxis=phase nxy=5.3 device=/xs options=wrap
# gpplt vis=all_1.vis yaxis=amp nxy=5.3 device=/xs options=wrap

#2

# after bandpass calibration check for zero phase and constant amplitude
# gpplt vis=all_2.vis yaxis=phase nxy=5.3 device=/xs options=bandpass
# gpplt vis=all_2.vis yaxis=amp nxy=5.3 device=/xs options=bandpass

#3

# after selfcal check phase vs. time; should be centred around zero
# gpplt vis=all_3.vis yaxis=phase nxy=5.3 device=/xs options=gains
# gpplt vis=all_3.vis yaxis=amp nxy=5.3 device=/xs options=gains

#4

# check phase vs baseline length to asses atmospheric decorrelation - should be flaring at 
# longer baselines but why!?
# smauvplt vis=all_2.vis device=/xs line=$WIDE axis=uvdist,phase options=nobase

#5

# check amplitude vs time also
# smauvplt vis=all_2.vis device=/xs axis=time,amp


# start cleaning, imaging phase calibrators, mapping and deconvolution

rm -rf beam0 map0 model0 clean0

# map the calibrator

#invert vis=all_4.vis map=map0 beam=beam0 imsize=128 line=chan,1,1,2,2 select="source(1733-130)"
# invert vis=all_4.vis map=map0 beam=beam0 imsize=128 line=chan,1,1,2,2 select="source($gname)"
invert vis=all_4.vis map=map0 beam=beam0 imsize=128 line=wide,1,1,2,2 select="source($gname)"
clean map=map0 beam=beam0 out=model0
restor map=map0 beam=beam0 model=model0 out=clean0

# map the source

rm -rf beam1 map1 model1 clean1
invert vis=all_4.vis map=map1 beam=beam1 line=wide,1,1,2,2 "select=source(GCRING)" options=mosaic imsize=129
mossdi map=map1 beam=beam1 out=model1
restor map=map1 beam=beam1 model=model1 out=clean1