######################################################################
#                                                                    #
# Use Case Script for Jupiter 6cm VLA                                #
#                                                                    #
# Last Updated STM 2008-03-25 (Beta Patch 1.0)                       #
# CALIBRATION ONLY                                                   #
#                                                                    #
######################################################################

import time
import os

# 
#=====================================================================
#
# This script has some interactive commands: scriptmode = True
# if you are running it and want it to stop during interactive parts.

scriptmode = True

#=====================================================================
#
# Set up some useful variables - these will be set during the script
# also, but if you want to restart the script in the middle here
# they are in one place:

pathname=os.environ.get('AIPSPATH').split()[0]
prefix='jupiter6cm.usecase'

msfile = prefix + '.ms'

gtable = prefix + '.gcal'
ftable = prefix + '.fluxscale'
atable = prefix + '.accum'

srcsplitms = prefix + '.split.ms'

#
#=====================================================================
#
# Get to path to the CASA home and stip off the name
pathname=os.environ.get('AIPSPATH').split()[0]

print "Assumes that flagged ms file is "+msfile
#
#=====================================================================
# Calibration
#=====================================================================
#
# Set the fluxes of the primary calibrator(s)
#
print '--Setjy--'
default('setjy')

print "Use setjy to set flux of 1331+305 (3C286)"

vis = msfile

#
# 1331+305 = 3C286 is our primary calibrator
field = '1331+305'     

# Setjy knows about this source so we dont need anything more

setjy()

#
# You should see something like this in the logger and casapy.log file:
#
# 1331+305  spwid=  0  [I=7.462, Q=0, U=0, V=0] Jy, (Perley-Taylor 99)
# 1331+305  spwid=  1  [I=7.51, Q=0, U=0, V=0] Jy, (Perley-Taylor 99)
# 

print "Look in logger for the fluxes (should be 7.462 and 7.510 Jy)"

#
#=====================================================================
#
# Initial gain calibration
#
print '--Gaincal--'
default('gaincal')

print "Solve for antenna gains on 1331+305 and 0137+331"
print "We have 2 single-channel continuum spw"
print "Do not want bandpass calibration"

vis = msfile

# set the name for the output gain caltable
gtable = prefix + '.gcal'
caltable = gtable

print "Output gain cal table will be "+gtable

# Gain calibrators are 1331+305 and 0137+331 (FIELD_ID 7 and 0)
# We have 2 IFs (SPW 0,1) with one channel each

# selection is via the field and spw strings
field = '1331+305,0137+331'
spw = ''

# a-priori calibration application
# atmospheric optical depth (turn off)
gaincurve = True
opacity = 0.0

# scan-based G solutions for both amplitude and phase
gaintype = 'G'
solint = 0.
calmode = 'ap'

# reference antenna 11 (11=VLA:N1)
refant = '11'

# minimum SNR 3
minsnr = 3

gaincal()

#
#=====================================================================
#
# Bootstrap flux scale
#
print '--Fluxscale--'
default('fluxscale')

print "Use fluxscale to rescale gain table to make new one"

vis = msfile

# set the name for the output rescaled caltable
ftable = prefix + '.fluxscale'
fluxtable = ftable

print "Output scaled gain cal table is "+ftable

# point to our first gain cal table
caltable = gtable

# we will be using 1331+305 (the source we did setjy on) as
# our flux standard reference
reference = '1331+305'

# we want to transfer the flux to our other gain cal source 0137+331
# to bring its gain amplitues in line with the absolute scale
transfer = '0137+331'

fluxscale()

# You should see in the logger something like:
#Flux density for 0137+331 in SpW=0 is:
#   5.42575 +/- 0.00285011 (SNR = 1903.7, nAnt= 27)
#Flux density for 0137+331 in SpW=1 is:
#   5.46569 +/- 0.00301326 (SNR = 1813.88, nAnt= 27)

#=====================================================================
#
# Interpolate the gains onto Jupiter (and others)
#
print '--Accum--'
default('accum')

print "This will interpolate the gains onto Jupiter"

vis = msfile

tablein = ''
incrtable = ftable
calfield = '1331+305, 0137+331'

# set the name for the output interpolated caltable
atable = prefix + '.accum'
caltable = atable

print "Output cumulative gain table will be "+atable

# linear interpolation
interp = 'linear'

# make 10s entries
accumtime = 10.0

accum()

#=====================================================================
#
# Correct the data
# (This will put calibrated data into the CORRECTED_DATA column)
#
print '--ApplyCal--'
default('applycal')

print "This will apply the calibration to the DATA"
print "Fills CORRECTED_DATA"

vis = msfile

# Start with the interpolated fluxscale/gain table
gaintable = atable

# Since we did gaincurve=True in gaincal, we need it here also
gaincurve = True
opacity=0.0

# select the fields
field = '1331+305,0137+331,JUPITER'
spw = ''
selectdata = False

# do not need to select subset since we did accum
# (note that correct only does 'nearest' interp)
gainfield = ''

applycal()

#
#=====================================================================
#
# Now split the Jupiter target data
#
print '--Split Jupiter--'
default('split')

vis = msfile

# Now we write out the corrected data to a new MS

# Make an output vis file
srcsplitms = prefix + '.split.ms'
outputvis = srcsplitms

print "Split Jupiter data into new ms "+srcsplitms

# Select the Jupiter field
field = 'JUPITER'
spw = ''

# pick off the CORRECTED_DATA column
datacolumn = 'corrected'

split()

#=====================================================================
#
# Export the Jupiter data as UVFITS
# Start with the split file.
#
print '--Export UVFITS--'
default('exportuvfits')

srcuvfits = prefix + '.split.uvfits'

print "Writing split Jupiter data to UVFITS file "+srcuvfits

vis = srcsplitms
fitsfile = srcuvfits

# Since this is a split dataset, the calibrated data is
# in the DATA column already.
datacolumn = 'data'

# Write as a multisource UVFITS (with SU table)
# even though it will have only one field in it
multisource = True

# Run asynchronously so as not to interfere with other tasks
# (BETA: also avoids crash on next importuvfits)
async = True

exportuvfits()

print '--Done with calibration--'
#
#=====================================================================