####################################################################
# THIS TUTORIAL SHOULD BE PERFORMED IN CASA 4.0.1                  #
# ANALYSIS UTILS WILL BE NEEDED                                    #
####################################################################
#                BAND 9 HANDS ON TUTORIAL                          #
####################################################################
#                                                                  #
####################################################################
#           THINGS YOU NEED TO KNOW ABOUT THE DATASET              #
#                                                                  #
# This data is an execution block of a Cycle 0 dataset, with the   #
# science target data deleted from it, used with kind permission   #
# of the PI. Hence only the calibrator data has been provided. A   #
# secondary phase-reference soource was observed, we will pretend  #
# that this is the target.                                         #
#                                                                  #
# The data were taken in the Cycle 0 extended array.               #
#                                                                  #
# The edges of the spectral windows have already been flagged      #
# and a model for Titan has been loaded. All standard apriori      #
# calibration and flagging has already been performed. The         #
# calibration here can be started at the bandpass calibration      #
# stage.                                                           # 
#                                                                  #
# The science goal of the observation was to image CO(J=6-5)       #
# emission. For the purposes of this tutorial PI is:               #
# (1) somewhat worried about adding noise to the science target    #
# from the bandpass calibration; (2) wants to have a recommendation#
# from the reducer about whether to pursue self-calibration;       #
# (3) wants to quantify the absolute position accuracy of the      #
# observation; (4) is not too worried about the accuracy of the    #
# absolute flux calibration. PI also intends to user-smooth the    #
# spectral resolution of the spectral line quite considerably.     #
#                                                                  #
# Bearing these things in mind, try to solve the problems in this  #
# dataset marked PROBLEM. Commands can either be cut and pasted    #
# into CASA, or run the script in steps e.g. in CASA               # 
# mysteps = [0]                                                    #
# execfile('band9handson_withproblems.py')                         #
####################################################################
# PROBLEM 1: Determine calibrator intents from plotms (mystep=0)   #
# PROBLEM 2: Determine spectral channel ranges that are not        #
# contaminated by atmospheric lines (mystep=3&4)                   #
# PROBLEM 3: Pre-bandpass phase-up channel selection (mystep=5)    #
# PROBLEM 4: Examine the gain table and decide if changes to       #
# the pre-bandpass phase-up are needed                             #
# PROBLEM 5: Create 2 bandpass tables and decide which to use      #
# PROBLEM 6: Fluxscale baseline selection method. Use ANTENNA or   #
# UVRANGE selection and choose which gain table to use.            #
# PROBLEM 7: Titan CO(6-5) emission: decide which spectral window  #
# to replace in scaling the flux from Titan, and which spectral    #
# window to replace it with                                        #
# PROBLEM 8: Fill in the quasar fluxes                             #
# PROBLEM 9: If this were the science target, would you recommend  #
# pursuing self-calibration?                                       #
# PROBLEM 10: Determine the absolute position uncertainty of the   #
# observation                                                      #
####################SET THE STEPS###################################
thesteps = []
step_title = {0: 'Listobs',
              1: 'Data inspection: amp vs time',
              2: 'Data inspection: Titan model',
              3: 'Data inspection: atmospheric lines',
              4: 'Selecting channels without atm lines',
              5: 'Pre-bandpass phase-up',
              6: 'Check the cal table:phases',
              7: 'Check the cal table: solution SNR',
              8: 'Create bandpass table 1: high spectral resolution',
              9: 'Examine the caltable (spw 3 only) (BP table 1)',
              10: 'Create bandpass table 2: smoothed spectral resolution',
              11: 'Examine bandpass table: smoothed resolution (BP table 2)',
              12: 'Overplot BPOLY bandpass on BP table 1',
              13: 'Choose the bandpass table',
              14: 'Fluxscale model and baseline selection 1',
              15: 'Plot cal table for baseline selection 1',
              16: 'Plot Titan amp vs. uvdistance data',
              17: 'Select uvrange for baseline selection 2',
              18: 'Fluxscale model and baseline selection 2',
              19: 'Plot the cal table for baseline selection 2',
              20: 'Gain table for fluxscale',
              21: 'Examining cal table (commented out for speed)',
              22: 'Plotting Titans spectral windows',
              23: 'Perform the fluxscaling',
              24: 'Set BP cal and phase cal fluxes for whole array',
              25: 'Gain calibration using all antennas',
              26: 'Apply the calibration!',
              27: 'Split out the data',
              28: 'Imaging',
              29: 'View image and image statistics',
              30: 'Fit a Gaussian to the Check Source'}

try:
  print 'List of steps to be executed ...', mysteps
  thesteps = mysteps
except:
  print 'global variable mysteps not set.'
if (thesteps==[]):
  thesteps = range(0,len(step_title))
  print 'Executing all steps: ', thesteps

####################Importing analysis utils#########################

sys.path.append('../analysis_scripts/') 

import analysisUtils as aU

#####################################################################
# PROBLEM 1: Use listobs to determine the primary calibration intent#  
# of 3c279, Titan and J1625-254. Ignore J1733-130 for the moment,   #
# which we will treat as a check source/science target since        #
# the real science target was removed.                              #
#                                                                   #
# Note: All Band 9 SBs at Cycle 1 should include a Check Source     #
# within 15 degrees of Science target. Calibration can be applied   #
# to this to check that the observation went OK.                    #
#####################################################################

######
mystep = 0

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  listobs('band9_calibration_hands-on.ms')

######
mystep = 1

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

# Plot amp vs time to get familiar with the data (spw 2 only
# to save plotting time)

  plotms(vis='band9_calibration_hands-on.ms',spw='2',xaxis='time',
      yaxis='amp',field='',avgchannel='3840',
      coloraxis='field', ydatacolumn='data',symbolsize=4)

######
mystep = 2

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]
  
  plotms(vis='band9_calibration_hands-on.ms',spw='',xaxis='uvdist',
    yaxis='amp',field='1',avgchannel='3840',
    coloraxis='spw', ydatacolumn='model',symbolsize=4)

#####################################################################
#PROBLEM 2: After some quick data inspection with the following     #
#commands, assess where terrestrial atmospheric lines are           #
#contaminating the dataset in spectral windows 0 and 3.             #
#####################################################################  

######
mystep = 3

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

# Check for the atmospheric lines using 2 strongest sources
# For iter plots use the green arrows at the bottom of plotms
# It might be useful to use zoom in plotms/change channel averaging
# etc

  plotms(vis='band9_calibration_hands-on.ms',spw='0,1,2,3',xaxis='channel',
         yaxis='amp',field='0,1',antenna='*&*',avgtime='1e8',avgscan=T,
         coloraxis='field',iteraxis='spw', xselfscale=T,ydatacolumn='data',
         symbolsize=4)

######
mystep = 4

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

# fill in the spw0 and spw3 atmospheric line-free regions below, needed later

  spw_without_atm='PROBLEM,1:200~3400,2:200~3640,PROBLEM'

######################################################################
# PROBLEM 3: Pre-bandpass phase-up.                                  #
# Take a look at the spectral channel selection in the command       #
# Do you see anything wrong with it? If so, edit to give similar     # 
# bandwidth for the moment.                                          #
######################################################################

######
mystep = 5

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

    
  spw_for_bp_phaseup='0:1536~2304,1:1536~2304,2:1536~2304,3:1536~2304' 

  os.system('rm -rf band9_calibration_hands-on.ms.ap_pre_bandpass') 

  gaincal(vis='band9_calibration_hands-on.ms',
    caltable = 'band9_calibration_hands-on.ms.ap_pre_bandpass',
    field = '0', # bandpass calibrator
    spw = spw_for_bp_phaseup,
    solint = 'int',
    refant = 'DV04',
    calmode = 'ap')

######################################################################
# PROBLEM 4: Examine the pre-bandpass gain table                     #
# Assess whether changes need to be made e.g. increase solint,       #
# reduce minsnr, increase bandwidth selected                         #
######################################################################

######
mystep = 6

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

#Examine the cal table, 1 polarization only for speed

  plotcal(caltable='band9_calibration_hands-on.ms.ap_pre_bandpass',
          xaxis='time',yaxis='phase',spw='',antenna='', 
          iteration='antenna',subplot=421,
          plotrange=[0,0,-180,180],poln='X')

######
mystep = 7

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

#Examine the cal table solution SNR, 1 polarization only for speed

  plotcal(caltable='band9_calibration_hands-on.ms.ap_pre_bandpass',
          xaxis='time',yaxis='snr',spw='',antenna='', 
          iteration='antenna',subplot=421,poln='X')

######################################################################
# PROBLEM 5: Create 2 bandpass tables and decide which to apply.     #
# Creation of a BPOLY bandpass table is also given as an example.    #
######################################################################

######
mystep = 8

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]
  
# creating bandpass table at high spectral resolution of observation

  os.system('rm -rf band9_calibration_hands-on.ms.bandpass')

  bandpass(vis = 'band9_calibration_hands-on.ms',
    caltable = 'band9_calibration_hands-on.ms.bandpass',
    field = '0', # bandpass calibrator
    solint = 'inf',
    combine = 'scan',  
    refant = 'DV04',
    solnorm = T,
    bandtype = 'B',
    gaintable = 'band9_calibration_hands-on.ms.ap_pre_bandpass')

######
mystep = 9

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]
  
# Use analysis utils to plot bandpass table: high spectral resolution

  aU.plotbandpass('band9_calibration_hands-on.ms.bandpass',
            xaxis='freq',yaxis='both', spw='3',
            antenna='', subplot=421, showatm=T,
            interactive=True,figfile='bandpass1.png')

######
mystep = 10

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

# creating bandpass table at a smoothed resolution

  os.system('rm -rf band9_calibration_hands-on.ms.bandpass_smooth20ch')
   
  bandpass(vis = 'band9_calibration_hands-on.ms',
    caltable = 'band9_calibration_hands-on.ms.bandpass_smooth20ch',
    field = '0', # bandpass calibrator
    solint = 'inf,20ch',
    combine = 'scan',
    refant = 'DV04',
    solnorm = T,
    bandtype = 'B',
    gaintable = 'band9_calibration_hands-on.ms.ap_pre_bandpass')

######
mystep = 11

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]
  
# use analysis utils to plot bandpass table: smoothed resolution

  aU.plotbandpass('band9_calibration_hands-on.ms.bandpass_smooth20ch',
            xaxis='freq',yaxis='both', spw='3',
            antenna='', subplot=421, showatm=T,
            interactive=True,figfile='bandpass2.png')

######
mystep = 12

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]
  
# create a BPOLY bandpass. This take too long to run, so intead
# of running the command, this cal table is provided for examination
  
#  os.system('rm -rf band9_calibration_hands-on.ms.bandpass_bpoly')
#  bandpass(vis='band9_calibration_hands-on.ms',
#    caltable='band9_calibration_hands-on.ms.bandpass_bpoly',
#    field='0',
#    spw='',
#    combine='scan',
#    refant='DV04',
#    solint='inf',solnorm=T,minblperant=4,fillgaps=5,
#    bandtype='BPOLY',degamp=7,degphase=13,
#    gaintable='band9_calibration_hands-on.ms.ap_pre_bandpass')

# overplot the provided BPOLY solution on bandpass table: high spectral
# resolution

  aU.plotbandpass(caltable='band9_calibration_hands-on.ms.bandpass',
                caltable2='band9_calibration_hands-on.ms.bandpass_bpoly',
                field='0',xaxis='freq',yaxis='both', spw='3',
                interactive=True,subplot=42)


###### FILL IN THE BANDPASS TABLE CHOSEN BELOW
mystep = 13

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  bandpasstable='PROBLEM'

#######################################################################
#PROBLEM 6: Highly-resolved flux calibrator                           #
# Create 2 different gain tables using a subset of antennas/uvrange   #
# in which only reasonably high signal to noise data are used. When   #
# limiting by uvrange, go out to about the 20% power point before     #
# the first null. Decide which gain table to use                      #
# #####################################################################

######
mystep = 14

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

# The automatic script generator has chosen these antennas selecting
# short baselines. Try making the cal table this way

  os.system('rm -rf band9_calibration_hands-on.ms.phase_short_int_antenna')
  
  gaincal(vis = 'band9_calibration_hands-on.ms',
    caltable = 'band9_calibration_hands-on.ms.phase_short_int_antenna',
    field = '0,1,3', # bandpasscal, fluxcal, phasecal
    selectdata = T,
    antenna = 'DA41,DA42,DV04,DV08,DV09,DV12&',                   
    solint = '30s',  # increased from solint='int'
    minsnr = 2.0,    # decreased from minsnr=3.0
    minblperant = 2, # decreased from minblperant=4
    spw=spw_without_atm, 
    refant = 'DV04',
    gaintype = 'G',  
    calmode = 'p',
    gaintable = bandpasstable)

  
######
mystep = 15

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  
# examine the solutions obtained by selecting by antenna

  plotcal(caltable='band9_calibration_hands-on.ms.phase_short_int_antenna',
          xaxis='time',yaxis='phase',spw='',antenna='', 
          iteration='antenna',subplot=421,plotrange=[0,0,-180,180], poln='X',
          figfile='baselineselection1_caltable.png')

######
mystep = 16

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

# plotting Titan data for spw 2

  plotms(vis='band9_calibration_hands-on.ms',field='1',
         xaxis='uvdist',yaxis='amp',spw='2',coloraxis='baseline',
         ydatacolumn='data',avgchannel='3840',symbolsize=4)

######
mystep = 17

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

# decide which uvrange to use in the gain cal instead
# of antenna selection

  uvrange1='PROBLEM'

######
mystep = 18

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

# Make gain table with restricted uv range as alternative

  os.system('rm -rf band9_calibration_hands-on.ms.phase_short_int_uvrange1')

  gaincal(vis = 'band9_calibration_hands-on.ms',
    caltable = 'band9_calibration_hands-on.ms.phase_short_int_uvrange1',
    field = '0,1,3', # bandpasscal, fluxcal, phasecal
    selectdata = T,
    uvrange = uvrange1,                   
    solint = '30s',  # increased from solint='int'
    minsnr = 2.0,    # decreased from minsnr=3.0
    minblperant = 2, # decreased from minblperant=4
    spw= spw_without_atm,
    refant = 'DV04',
    gaintype = 'G',  
    calmode = 'p',
    gaintable = bandpasstable)

######
mystep = 19

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  
# plotting the solutions for case 2

  plotcal(caltable='band9_calibration_hands-on.ms.phase_short_int_uvrange1',
          xaxis='time',yaxis='phase',spw='',antenna='', 
          iteration='antenna',subplot=421,plotrange=[0,0,-180,180], poln='X',
          figfile='baselineselection2_caltable.png')

##################FILL IN THE TABLE YOU WANT TO USE
##################DEPENDING ON WHICH YOU CHOOSE, THEN  
##################COMMENT OUT EITHER ANTENNA= OR UVRANGE=
################## IN MYSTEP 20  
mystep = 20

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  phaseshortinttable='PROBLEM'

  os.system('rm -rf band9_calibration_hands-on.ms.ampli_short_inf')

  gaincal(vis = 'band9_calibration_hands-on.ms',
    caltable = 'band9_calibration_hands-on.ms.ampli_short_inf',
    field = '0,1,3', # bandpass cal, flux cal, phase cal
    selectdata = T,
    antenna = 'DA41,DA42,DV04,DV08,DV09,DV12&', #PROBLEM: comment this out? 
    uvrange =  uvrange1,  # PROBLEM: or comment this out?               
    solint = 'inf',
    spw=spw_without_atm, 
    minblperant = 2, # decreased from minblperant = 2
    minsnr = 2.0,    # decreased from minsnr=3.0
    refant = 'DV04',
    gaintype = 'T',
    calmode = 'ap',
    gaintable = [bandpasstable, phaseshortinttable])

######
mystep = 21

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

# examine the amp solutions obtained using the chosen restriction
# Commented out here for speed

#  plotcal(caltable='band9_calibration_hands-on.ms.ampli_short_inf',
#          xaxis='time',yaxis='amp',spw='',antenna='', 
#          iteration='antenna',subplot=421, poln='X')

############################################################# 
# PROBLEM 7: TITAN CO(J=6-5) EMISSION                       #
# One of the spw of Titan has a broad contaminating CO line #
# at about 691.5 GHz. Determine the window with the problem #
# from the following plots, and decide which number spw     #
# to replace it with in the fluxscale                       #
#############################################################

######
mystep = 22

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

# Titan spectral windows
  
  plotms(vis='band9_calibration_hands-on.ms',spw='0,1,2,3',xaxis='frequency',
       yaxis='amp',field='1', antenna='DA41&DA42',avgtime='1e8',avgscan=T,
       coloraxis='spw', xselfscale=T,ydatacolumn='corrected',avgchannel='100',
       symbolsize=4)

###### FILL IN THE NUMBERS HERE FOR USING IN REFSPWMAP
mystep = 23

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]


  fscale_spwmap=[9,9,9,9] # PROBLEM: replace '9' with suitable spw numbers

  os.system('rm -rf band9_calibration_hands-on.ms.flux_short_inf')

  fluxscale(vis = 'band9_calibration_hands-on.ms',
    caltable = 'band9_calibration_hands-on.ms.ampli_short_inf',
    fluxtable = 'band9_calibration_hands-on.ms.flux_short_inf',
    refspwmap =  fscale_spwmap,
    reference = '1') # Flux calibrator

############################################################# 
# PROBLEM 8: SETTING BANDPASS AND PHASE CAL FLUXES          #
# For speed in this tutoral, take a rough, average flux     #
# density for all spw from the logger values for each source#.
# For greater accuracy, you can set each spw individually.In#
# cases like J1625-254, where the values seem to have large #
# uncertainties, you could make an error-weighted fit       #
# assuming a single spectral index and estimate consistent  #
# values for each spw.                                      #
#############################################################


###### 
mystep = 24

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  delmod('band9_calibration_hands-on.ms', field='0,3')

  flux3c279=[PROBLEM,0,0,0]
  flux1625=[PROBLEM,0,0,0]

  setjy(vis='band9_calibration_hands-on.ms',field='3c279',fluxdensity=flux3c279,usescratch=F)
  setjy(vis='band9_calibration_hands-on.ms',field='J1625-254',fluxdensity=flux1625,usescratch=F)


###### 
mystep = 25

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

# perform the gain calibration using all antennas

  
  os.system('rm -rf band9_calibration_hands-on.ms.phase_int')

  gaincal(vis = 'band9_calibration_hands-on.ms',
    caltable = 'band9_calibration_hands-on.ms.phase_int',
    field = '0,1,3', # bandpass cal, flux cal, phase cal
    spw= spw_without_atm,
    solint = '30s', # increased from solint='int'
    minsnr = 2.0, # decreased from minsnr='3.0'
    refant = 'DV04',
    gaintype = 'G',
    calmode = 'p',
    gaintable = bandpasstable)
  
  os.system('rm -rf band9_calibration_hands-on.ms.flux_inf')

  gaincal(vis = 'band9_calibration_hands-on.ms',
    caltable = 'band9_calibration_hands-on.ms.flux_inf',
    field = '0,1,3', # bandpass cal, flux cal, phase cal
    spw= spw_without_atm,
    solint = 'inf',
    refant = 'DV04',
    gaintype = 'T',
    calmode = 'ap',
    gaintable = [bandpasstable, 'band9_calibration_hands-on.ms.phase_int'])
  
  os.system('rm -rf band9_calibration_hands-on.ms.phase_inf')
  
  gaincal(vis = 'band9_calibration_hands-on.ms',
    caltable = 'band9_calibration_hands-on.ms.phase_inf',
    field = '0,1,3', # bandpass cal, flux cal, phase cal
    spw= spw_without_atm,
    solint = 'inf',
    refant = 'DV04',
    gaintype = 'G',
    calmode = 'p',
    gaintable = bandpasstable)

# If you have time, could generate plot commands to examine the tables  


###### 
mystep = 26

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  for i in ['0', '1']: # bandpass cal, flux cal
   applycal(vis = 'band9_calibration_hands-on.ms',
    field = i,
    gaintable = [bandpasstable, 'band9_calibration_hands-on.ms.phase_int', 'band9_calibration_hands-on.ms.flux_inf'],
    gainfield = ['', i, i],
    interp = 'linear,linear',
    calwt = F,
    flagbackup = F)
  
  applycal(vis = 'band9_calibration_hands-on.ms',
    field = '2,3', # Check Source, phase cal
    gaintable = [bandpasstable, 'band9_calibration_hands-on.ms.phase_inf','band9_calibration_hands-on.ms.flux_inf'],
    gainfield = ['','3', '3'], # phase cal
    interp = 'linear,linear',
    calwt = F,
    flagbackup = F)

# use plotms to see if everything looks reasonable

###### 
mystep = 27

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

# Image the check source. First split the data out at lower spectral resolution to make the imaging go faster

# It isn't possible to select more than one channel range per spectral window, so use this

  os.system('rm -rf checksource.cont.ms')

  split(vis = 'band9_calibration_hands-on.ms',
        outputvis = 'checksource.cont.ms',
        spw = '0:200~3639,1:200~3639,2:200~3639,3:200~3639',
        field='2',width = 100,
        datacolumn='corrected')

###### 
mystep = 28

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]


# If you have time, set interactive to True and comment out the
# 'mask=' line, to set the clean mask and perform imaging interactively.
# If you do this, the imsize can be set to [216,216] to speed things up.
# However, then check the calstat parameters in the next step

  os.system('rm -rf checksource.cont.clean*')
 
  clean(vis = 'checksource.cont.ms',
     imagename = 'checksource.cont.clean',
     mode = 'mfs',
     interactive = F,
     imsize = [360,360],
     mask =  [160,165,195,190],
     cell = '0.05arcsec',
     weighting = 'briggs',
     niter = 10000,
     threshold = '10mJy',
     robust = 0.5)


###### 
mystep = 29

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

  imview('checksource.cont.clean.image')

  calstat=imstat(imagename='checksource.cont.clean.image', region='', box='234,65,326,312')
  rms=(calstat['rms'][0])
  print '>> rms in continuum image: '+str(rms)
  calstat=imstat(imagename='checksource.cont.clean.image', region='')
  peak=(calstat['max'][0])
  print '>> Peak in continuum image: '+str(peak)
  print '>> Dynamic range in continuum image: '+str(peak/rms)

########################################################################
# PROBLEM 9: if this were the real science target would you recommend  #
# that PI pursue self-calibration?                                     #
########################################################################

###### 
mystep = 30

if(mystep in thesteps):
  casalog.post('Step '+str(mystep)+' '+step_title[mystep],'INFO')
  print 'Step ', mystep, step_title[mystep]

# Going back to using this as a check source, use the check source as a way
# to determine the absolute position uncertainty of the observation

  imfit(imagename='checksource.cont.clean.image', box='160,165,195,190')

########################################################################
# PROBLEM 10: compare the RA and dec given in the logger with the      # 
# position for J1733-130 (nrao530) in the SMA calibrator catalogue.    #
# What is the absolute position uncertainty of this observation?       #
#                                                                      #
# http://sma1.sma.hawaii.edu/callist/callist.html                      #
#                                                                      #
########################################################################  
#                                                                      #
#                         THE END                                      #
#                                                                      #
########################################################################