#!/bin/csh
# reduce.csh for L1157 continuum and 2-1 CO  29sep2011

  # notes: 3 wide windows, one 31 MHz with 12CO in USB, C18O in LSB
  # fluxcal = MWC349, passcal = 3c454.3, gaincal = 1927+739
  # grid data taken approx every 20 min (3 x 10sec each); purpose(P) = grids in
  # tau ~ 0.45, rmspath ~ 200 um
  # script crash midway through apparently caused by radio pointing

  # known problems: bands 5,7 on C2 have unstable xyphase

set RAW1 = ct002.L1157.3.miriad
set RAW2 = ct002.L1157.5.miriad
set REFANT = 1
set SHORT = L1157
goto start

# -------- (0) untar the data -------- #

rm -r $RAW1
tar xvf $RAW1.tar.gz
rm -r $RAW2
tar xvf $RAW2.tar.gz

# -------- (1) initial data inspection -------- #

uvindex vis=$RAW1 > $SHORT.3.log
uvindex vis=$RAW2 > $SHORT.5.log
varplt vis=$RAW1 yaxis=rmspath device=/xs		# ~ 200 um
varplt vis=$RAW1 yaxis=tau230 device=/xs	    # ~ 0.45
uvlist options=spectra vis=$RAW                 # review correlator setup

# plot tsys; must separate file into LL and RR or risk overflowing buffers
rm -r win1LL
uvcat vis=$RAW1,$RAW2 out=win1LL select='win(1),pol(LL),-source(noise)' \
	options=nocal,nopol,nopass
varplt vis=win1LL yaxis=systemp device=/xs yrange=0,2000	# all stable, about 700 K
rm -r win1RR
uvcat vis=$RAW1,$RAW2 out=win1RR select='win(1),pol(RR),-source(noise)' \
	options=nocal,nopol,nopass
varplt vis=win1RR yaxis=systemp device=/xs yrange=0,2000	# C7,C12 unstable! why?

# -------- (2) fit XYphase on 10-m antennas -------- #

/o/plambeck/1mmDualPol/Calibration/mfcal vis=$RAW1 interval=10000 \
    select='auto,-ant(7,8,9,10,11,12,13,14,15),purpose(P)' stokes=LR
smagpplt vis=$RAW1 options=bandpass,nofit,wrap device=/xs yrange=-180,180 \
    xaxis=chan yaxis=phase

# plot record by record to search for glitches
uvspec vis=$RAW1 select='auto,-ant(7,8,9,10,11,12,13,14,15),purpose(P)' \
	axis=chan,phase device=/xs  nxy=3,2 yrange=-60,60 options=all
  # C2 win(5,7) are unstable

# -------- (3) rewrite data to apply XYphase -------- #

rm -r $SHORT.xy
uvcat vis=$RAW1,$RAW2 out=$SHORT.xy options=nocal,nopol,unflagged \
  select='-source(noise)' 

  # note: at this point phase(LL) = phase(RR) on 10-m, but not 6-m, telescopes
  # we have not yet solved for a passband - there will be band to band offsets

# --- flag windows 5,7,13,15 for C2 because of xyphase instability --- #
uvflag vis=$SHORT.xy flagval=flag select='ant(2)' line=chan,238,95,1,1
uvflag vis=$SHORT.xy flagval=flag select='ant(2)' line=chan,238,427,1,1

# -------- (4) strip out wideband windows, fit and apply passband to them -------- #

rm -r $SHORT.xy2
uvcat vis=$SHORT.xy out=$SHORT.xy2 options=nocal,nopass,nopol \
	select='win(1,3,5,9,11,13),-auto,-source(noise),-purpose(P)' 

mfcal vis=$SHORT.xy2 select='source(3c454.3),-purpose(P)' interval=.1 refant=$REFANT
smagpplt vis=$SHORT.xy2 options=bandpass,nofit,wrap device=/xs yrange=-180,180 \
    xaxis=chan yaxis=phase
$<
gpplt vis=$SHORT.xy2 options=wrap yaxis=phase device=/xs yrange=-180,180 nxy=5,3
$<
gpplt vis=$SHORT.xy2 yaxis=amp device=/xs yrange=0,3 nxy=5,3

# rewrite data to apply passband, but not gains
rm -r $SHORT.pb
uvcat vis=$SHORT.xy2 out=$SHORT.pb options=nocal,nopol 

# -------- (5) fit gains to calibrators -------- #

mfcal vis=$SHORT.pb refant=$REFANT interval=12 options=nopass,nopassol \
	select='-source(L1157),pol(LL,RR)'
gpplt vis=$SHORT.pb yaxis=phase options=wrap device=/xs yrange=-180,180 nxy=5,3
gpplt vis=$SHORT.pb yaxis=amp device=/xs yrange=0,3 nxy=5,3

# -------- (6) solve for leakages on 1927+739  -------- #

gpcal vis=$SHORT.pb options=circular,qusolve,noxy interval=0.1 \
   refant=$REFANT select='source(1927+739),-purpose(P)' line=chan,1,1,282 

# -------- (8) dust continuum map -------- #

set MAP = dust
rm -r  $MAP.I.mp $MAP.Q.mp $MAP.U.mp $MAP.V.mp $MAP.bm
invert vis=$SHORT.pb line=chan,6,1,47,47 \
  map=$MAP.I.mp,$MAP.Q.mp,$MAP.U.mp,$MAP.V.mp beam=$MAP.bm stokes=I,Q,U,V sup=0 \
  'select=source(L1157)' options=systemp,mfs cell=0.25 imsize=415
	# using imsize=415 to match co.mom1.cm below
implot in=$MAP.I.mp device=/xs

rm noiseList
foreach MP ($MAP.I $MAP.Q $MAP.U $MAP.V)
  rm -r $MP.sl
  clean map=$MP.mp beam=$MAP.bm out=$MP.sl niters=2000
  rm -r $MP.cm
  restor map=$MP.mp beam=$MAP.bm model=$MP.sl out=$MP.cm
  cgdisp in=$MP.cm device=/xs region=quarter labtyp=arcsec
  echo " " >> noiseList
  echo $MP".cm" >> noiseList
  imlist options=stat in=$MP.cm region='arcsec,box(20,-20,5,20)' | tail -2 >> noiseList
  $<
end
tail -20 noiseList

rm -r $MAP.poli.cm $MAP.polm.cm $MAP.pa.cm
impol poli=$MAP.poli.cm polm=$MAP.polm.cm pa=$MAP.pa.cm sigma=.0012,.023 \
    in=$MAP.Q.cm,$MAP.U.cm,$MAP.I.cm sncut=2,1.

cgdisp device=/xs in=$MAP.I.cm,$MAP.poli.cm,$MAP.pa.cm type=contour,amp,angle \
  region='arcsec,box(-8,-8,8,8)' options=noepoch,full,rot90 labtyp=arcsec \
  slev=a,.006 cols1=1 levs1=-2,2,4,6,8,10,12,14,16,18,20,22,24,26,28,30,32,34,36,38,40,42,44 \
  vecfac=1.2,4,4 beamtyp=b,l,4 lines=1,1,10 device=L1157D.ps/vcps 
end: