######################################################################
#                                                                    #
# Use Case Script for reducing NGC 2403 HI line data                 #
#                                                                    #
# Original version GvM 2007-11-30 (Beta)                             #
# This version reads data from 4 VLA archive files                   #
#                                                                    #
######################################################################

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='n2403'

msfile = prefix + '.ms'

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

# This is where the UVFITS data will be
#fitsdata=pathname+'/data/nrao/VLA/planets_6cm.fits'
#fitsdata='/home/sandrock2/smyers/NAUG2/Data/VLA_CONT/FLUX99-6CM.CBAND'
#fitsdata='jupiter6cm.uvfits'

# The prefix to use for all output files
prefix='n2403'

# Clean up old files
os.system('rm -rf '+prefix+'*')

#=====================================================================
# Data Import
#=====================================================================
#
# Import the data from VLA archive files to MS
#
print "--importvla--"
print ""
print "Use importvla to read 4 VLA archive files and write the data"
print "into a Measurement Set (MS).  This will take a while ..."
print ""

# Set up the MS filename and save as new global variable
msfile = prefix + '.ms'

print "MS will be called "+msfile

default('importvla')

archivefiles=['AS649_1','AS649_2','AS649_3','AS649_4']
vis = msfile
importvla()

print ""
# Pause script if you are running in scriptmode
if scriptmode:
    user_check=raw_input('done with importvla, hit Return to continue script\n')


#=====================================================================
# List a summary of the MS
#=====================================================================
#
#
print "--listobs--"
print ""
print "Use listobs to print verbose summary to logger"
print "see the logger window for the listobs output"

# Don't default this one and make use of the previous setting of
# vis.  Remember, the variables are GLOBAL!

# You may wish to see more detailed information, in this case
# use the verbose = True option
#
verbose = True

listobs()

# The listobs output will be displayed in your logger window and in
# the casapy.log file

print ""
# Pause script if you are running in scriptmode
if scriptmode:
    user_check=raw_input('done with listobs, hit Return to continue script\n')


#=====================================================================
# Fill the model column for flux density calibrators
#=====================================================================
#
print "--setjy--"
print ""
print "find the flux of the flux calibrators, and write it to the"
print "column labeled MODEL_DATA"

default('setjy')

vis=msfile
field='1,3,4'     # note: field 1 is the source NGC 2403, and field 2 is
                  # the phase calibrator 0841+708'
spw='0:5~112'     # use spectral window 0 (which is the only one).  In that
                  # window, use channels 5 - 112 (ignoring edge channels)

setjy()

print ""
# Pause script if you are running in scriptmode
if scriptmode:
    user_check=raw_input('done with setjy, hit Return to continue script\n')

#=====================================================================
# Plot the source data
#=====================================================================
#
print "--plotxy--"
print ""
print "plots source visibilities, this time RR only"

default('plotxy')

vis=msfile
field='0'
spw='0:5~112'
correlation='RR'
average='chan'

plotxy()

print ""
print "Showing NGC 2403, RR"
print "clearly, there is one bad time interval"
print "click Mark Region, then draw a rectangle around the bad data"
print "click Locate to list the points, and Flag to flag them.  After"
print "Flag, the plot will be redrawn with the appropriate scale.  Do"
print "more flagging if necessary"
print ""

print ""
# Pause script if you are running in scriptmode
if scriptmode:
    user_check=raw_input('when happy with flagging, Return to continue script\n')

print ""
print "now NGC 2403, LL"
print ""

correlation='LL'

plotxy()

print ""
print "Now showing NGC 2403, LL"
print "one bad point is visible in the upper right corner"
print "click Mark Region, then draw a rectangle around it"
print "click Locate to list the point, then Flag to flag them"
print "the plot will be redrawn with an appropriate scale"
print "do more flagging if necessary"
print ""
# Pause script if you are running in scriptmode
if scriptmode:
    user_check=raw_input('when done with flagging, Return to continue script\n')


print "It is a good idea to save the flagging at certain times"
print "First, list the flag versions using flagmanager"

#=====================================================================
# Save flagging done up to this point
#=====================================================================
#
print ""
print "--flagmanager--"
print ""
print "first, list the current flag versions"
print ""

# first we list the current flagging tables

mode='list'

flagmanager()

print ""
print "then, we save the flagging we just did"

# then we save the flagging we just did

mode='save'
versionname='afterplotxy'
comment='flags after running plotxy'
merge='replace'

flagmanager()

# and now we list one more time to show the changes made

print ""
print "then, list one more time to show the changes"


mode='list'
flagmanager()

print "now we will determine the antenna gains"
print ""
# Pause script if you are running in scriptmode
if scriptmode:
    user_check=raw_input('Return to continue script\n')

#=====================================================================
# Determine antenna gains
#=====================================================================
#
print "--gaincal--"
print ""
print "creates user defined table containing gain solutions"
print "once for each calibrator since uv ranges are different."
print "Note: append is False first, then True"
print ""

default('gaincal')
vis=msfile
caltable= prefix + '.gcal'
field='1'
spw='0:5~112'
selectdata=True
uvrange='0~40kl'
append=False

gaincal()

field='2'
uvrange='0~20kl'
append=True

print "starting field 1"
print ""

gaincal()


field='3'
uvrange='0~50kl'

print "starting field 3"
print ""

gaincal()

field='4'
uvrange=''

print "starting field 4"
print ""

gaincal()


print ""
print "next, we will plot some of the antenna gains"
print ""
# Pause script if you are running in scriptmode
if scriptmode:
    user_check=raw_input('Return to continue script\n')

#=====================================================================
# Plot antenna gains
#=====================================================================
#
print "--plotcal--"

print "first we plot the amplitude gains for antennas 9 - 12"

default('plotcal')
caltable= prefix + '.gcal'
xaxis='time'
yaxis='amp'
field='2'
antenna='9~12'

plotcal()

print ""
# Pause script if you are running in scriptmode
if scriptmode:
    user_check=raw_input('Return to continue script\n')

print "then, we plot R and L just for antenna 9 in separate plots on"
print "the same page.  Note use of the subplot parameter"

yaxis='phase'
antenna='9'
poln='R'
subplot=211
plotcal()
poln='L'
subplot=212
plotcal()

print ""
print "Next, we will determine the flux of 0841+708"

print ""
# Pause script if you are running in scriptmode
if scriptmode:
    user_check=raw_input('Return to continue script\n')


#=====================================================================
# Bootstrap flux of 0841+708
#=====================================================================
#
print "--fluxscale--"
print ""
print "determines flux based on gains and flux calibrator fluxes"
print "see Log window for flux value found"
default('fluxscale')
vis=msfile
caltable= prefix + '.gcal'
fluxtable=prefix + '.fcal'
reference='1,3~4'
transfer='2'

fluxscale()

print ""
print "Next, we will determine bandpasses for the three bandpass calibrators"

print ""
# Pause script if you are running in scriptmode
if scriptmode:
    user_check=raw_input('Return to continue script\n')


#=====================================================================
# Solves for bandpass, writes it to table
#=====================================================================
#
print "--bandpass--"
print ""
print "determine bandpass"

default('bandpass')

vis=msfile
caltable= prefix + '.bcal'
field='1,3~4'
solnorm=True

bandpass()

print "Next we plot bandpass solutions for a few selected antennas;"
print ""
# Pause script if you are running in scriptmode
if scriptmode:
    user_check=raw_input('Return to continue script\n')

#=====================================================================
# Plot bandpass
#=====================================================================
#
print "--plotcal--"
print ""
print "First we plot solutions for antennas 9-12 for field 1 only"

default('plotcal')

vis=msfile
caltable= prefix + '.bcal'
xaxis               =     'chan'
yaxis               =      'amp'
field               =        '1'
antenna             =     '9~12'
plotrange           = [-1, -1, 0.9, 1.15]

plotcal()

print ""
# Pause script if you are running in scriptmode
if scriptmode:
    user_check=raw_input('Return to continue script\n')
print ""

field='1,3~4'
antenna='25'
iteration='field'

plotcal()

print "we iterate over all three fields, just for antenna 25 using"
print "the iteration parameter"
print ""

print "Make sure to click Next to go through the fields before hitting return"
print ""
print "note the galactic absorption in the first two fields"
print "only field 4 (1331+305) is free of absorption"
print "for now, we will use only bandpass solutions for field 4"
print ""
# Pause script if you are running in scriptmode
if scriptmode:
    user_check=raw_input('Return to run applycal\n')


#=====================================================================
#Apply calibration - results go to corrected_data column
#=====================================================================
print "--applycal--"
print ""
print "applies supplied calibration tables (gain and bandpass) and"
print "writes corrected data to column labeled CORRECTED_DATA.  This"
print "may take a while ..."

default('applycal')

vis=msfile
spw='0:5~112'
gaintable=[prefix+'.fcal',prefix+'.bcal']
gainfield=['1', '4']

applycal()

print ""
print "Next, we will Fourier transform the corrected data using invert"

# Pause script if you are running in scriptmode
if scriptmode:
    user_check=raw_input('Return to continue script\n')


#=====================================================================
# Image a few channels
#=====================================================================
#
print "--invert--"
print "Fourier transforms data in column labeled CORRECTED_DATA"
print "and creates an output image and beam"
print ""
print "For now, we won't clean, and to make things go faster we image"
print "two channels only"

default('invert')

vis=msfile
imagename=prefix
mode='channel'
nchan=2
step=5
start=50
imsize=[512,512]
cell=['4arcsec','4arcsec']
field='0'
spw='0:5~112'

invert()


# Pause script if you are running in scriptmode
if scriptmode:
    user_check=raw_input('Return to continue script\n')


#=====================================================================
# Inspect the result of invert
#=====================================================================
#
print "--viewer--"
print "allows viewing both visibilities and images. Here: image"
print ""
print "Here we inspect the two channels we just created"
print ""

default('viewer')

infile=prefix+'.dirty'

viewer()

print "note that HI emission is visible, but that there is a structure"
print "superposed suggesting bad data.  Maybe using the channel average for"
print "flagging was not enough.  Sometimes data can be bad in individual"
print "channels but averaged out if the phases are random"
print ""
print "please exit viewer first!"
print ""
print "next, we will run plotxy again but now using only one channel"

print ""
# Pause script if you are running in scriptmode
if scriptmode:
    user_check=raw_input('Return to continue script\n')


#=====================================================================
# Use plotxy again, but now on a single channel only
#=====================================================================
#
print "--plotxy  --"
print ""
print "plots visibility data"
print ""

execfile('plotxy.last')

spw='0:50'
correlation=''

plotxy()

print "Note that bad data not flagged previously now show up.  Don't"
print "flag these points interactively; instead, use the Locate button"
print "to see what's going on.  Turns out, there are some bad time stamps"
print "in both RR and LL, and a few in just LL.  In addition, antenna 0"
print "(LL) is suspect throughout the observation"
print ""
print "When done playing around with the CASA plotter, we will use the"
print "non-interactive task flagdata to flag these data"
print ""
print "first: the bad timeranges:"
print ""
# Pause script if you are running in scriptmode
if scriptmode:
    user_check=raw_input('Return to continue script\n')

#=====================================================================
# Flag data non-interactively
#=====================================================================
#
print "--flagdata--"
print "flags data based on a number of criteria"
print ""
print "first we flag a narrow time range for both correlations"
print ""

default('flagdata')

vis=msfile
spw='0'
timerange='03:52:44~03:52:46'
mode='manualflag'

flagdata()

print ""
print "then, we flag antenna 0 for correlation LL over the whole time range"
print ""


antenna='0'
correlation='LL'
timerange=''

flagdata()

print ""
print "done with flagdata, run invert again to see if things improved"

print ""
# Pause script if you are running in scriptmode
if scriptmode:
    user_check=raw_input('Return to continue script\n')

# first delete previously made images
os.system('rm -rf '+prefix+'.dirty')
os.system('rm -rf '+prefix+'.beam')

#=====================================================================
# Image selected channels after latest flagging
#=====================================================================
#
print "--invert--"

execfile('invert.last')

nchan=55
start=32

invert()

print ""
print "invert finished; now run the viewer again to see if things look better"
print ""
# Pause script if you are running in scriptmode
if scriptmode:
    user_check=raw_input('Return to continue script\n')

#=====================================================================
# Inspect to see if the latest flagging made things better
#=====================================================================
#
print "--viewer--"

execfile('viewer.last')

viewer()

print ""
print "This marks the end of the script.  Thanks for hanging in there."
print "Soon we will add continuum subtraction in the uv plane and imaging"
print "with cleaning"