######################################################################
#                                                                    #
# Use Case Script for Jupiter 6cm VLA                                #
#                                                                    #
# Last Updated STM 2008-03-25 (Beta Patch 1.0)                       #
# IMAGING 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 for output from the previous scripts
prefix='jupiter6cm.usecase'

msfile = prefix + '.ms'

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

srcsplitms = prefix + '.split.ms'

# New prefix for output from this script
prefix='jupiter6cm.usecase.img'

clnimsize = [288,288]
clncell = [4.,4.]

imname1 = prefix + '.clean1'
clnimage1 = imname1+'.image'
clnmodel1 = imname1+'.model'
clnresid1 = imname1+'.residual'
clnmask1  = imname1+'.clean_interactive.mask'

selfcaltab1 = srcsplitms + '.selfcal1'

imname2 = prefix + '.clean2'
clnimage2 = imname2+'.image'
clnmodel2 = imname2+'.model'
clnresid2 = imname2+'.residual'
clnmask2  = imname2+'.clean_interactive.mask'

selfcaltab2 = srcsplitms + '.selfcal2'
smoothcaltab2 = srcsplitms + '.smoothcal2'

imname3 = prefix + '.clean3'
clnimage3 = imname3+'.image'
clnmodel3 = imname3+'.model'
clnresid3 = imname3+'.residual'
clnmask3  = imname3+'.clean_interactive.mask'

#
#=====================================================================
# FIRST CLEAN / SELFCAL CYCLE
#=====================================================================
#
# Now clean an image of Jupiter
#
print '--Clean 1--'
default('clean')

# Pick up our split source data
vis = srcsplitms

# Make an image root file name
imname1 = prefix + '.clean1'
imagename = imname1

print "Output images will be prefixed with "+imname1

# Set up the output continuum image (single plane mfs)
mode = 'mfs'
stokes = 'I'

print "Will be a single MFS continuum image"

# NOTE: current version field='' doesnt work
field = '*'

# Combine all spw
spw = ''

# This is D-config VLA 6cm (4.85GHz) obs
# Check the observational status summary
# Primary beam FWHM = 45'/f_GHz = 557"
# Synthesized beam FWHM = 14"
# RMS in 10min (600s) = 0.06 mJy (thats now, but close enough)

# Set the output image size and cell size (arcsec)
# 4" will give 3.5x oversampling
# 280 pix will cover to 2xPrimaryBeam
# clean will say to use 288 (a composite integer) for efficiency
clnalg = 'clark'
clnimsize = [288,288]

# double for CS Clean
#clnalg = 'csclean'
#clnimsize = [576,576]

clncell = [4.,4.]

alg = clnalg
imsize = clnimsize
cell = clncell

# NOTE: will eventually have an imadvise task to give you this
# information

# Standard gain factor 0.1
gain = 0.1

# Fix maximum number of iterations
niter = 10000

# Also set flux residual threshold (0.04 mJy)
# From our listobs:
# Total integration time = 85133.2 seconds
# With rms of 0.06 mJy in 600s ==> rms = 0.005 mJy
# Set to 10x thermal rms
threshold=0.05

# Note - we can change niter and threshold interactively
# during clean

# Set up the weighting
# Use Briggs weighting (a moderate value, on the uniform side)
weighting = 'briggs'
rmode = 'norm'
robust = 0.5

# No clean mask
mask = ''

# Use interactive clean mode
cleanbox = 'interactive'

# Moderate number of iter per interactive cycle
npercycle = 100

clean()

# When the interactive clean window comes up, use the right-mouse
# to draw rectangles around obvious emission double-right-clicking
# inside them to add to the flag region.  You can also assign the
# right-mouse to polygon region drawing by right-clicking on the
# polygon drawing icon in the toolbar.  When you are happy with
# the region, click 'Done Flagging' and it will go and clean another
# 100 iterations.  When done, click 'Stop'.

# Set up variables
clnimage1 = imname1+'.image'
clnmodel1 = imname1+'.model'
clnresid1 = imname1+'.residual'
clnmask1  = imname1+'.clean_interactive.mask'

print ""
print "----------------------------------------------------"
print "Clean"
print "Final clean model is "+clnmodel1
print "Final restored clean image is "+clnimage1
print "The clean residual image is "+clnresid1
print "Your final clean mask is "+clnmask1

print ""
print "This is the final restored clean image in the viewer"
print "Zoom in and set levels to see faint emission"
print "Use rectangle drawing tool to box off source"
print "Double-click inside to print statistics"
print "Move box on-source and get the max"
print "Calcualte DynRange = MAXon/RMSoff"
print "I got 1.060/0.004 = 270"
print "Still not as good as it can be - lets selfcal"
print "Close viewer panel when done"

#
#---------------------------------------------------------------------
#
# If you did not do interactive clean, bring up viewer manually
viewer(clnimage1,'image')

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

# You can use the right-mouse to draw a box in the lower right
# corner of the image away from emission, the double-click inside
# to bring up statistics.  Use the right-mouse to grab this box
# and move it up over Jupiter and double-click again.  You should
# see stuff like this in the terminal:
#
# jupiter6cm.usecase.clean1.image     (Jy/beam)
# 
# n           Std Dev     RMS         Mean        Variance    Sum
# 4712        0.003914    0.003927    0.0003205   1.532e-05   1.510     
# 
# Flux        Med |Dev|   IntQtlRng   Median      Min         Max
# 0.09417     0.002646    0.005294    0.0001885   -0.01125    0.01503   
#
#
# On Jupiter:
#
# n           Std Dev     RMS         Mean        Variance    Sum
# 3640        0.1007      0.1027      0.02023     0.01015     73.63     
# 
# Flux        Med |Dev|   IntQtlRng   Median      Min         Max
# 4.592       0.003239    0.007120    0.0001329   -0.01396    1.060     
#
# Estimated dynamic range = 1.060 / 0.003927 = 270 (poor)
#
# Note that the exact numbers you get will depend on how deep you
# take the interactive clean and how you draw the box for the stats.

#
#---------------------------------------------------------------------
#
# Self-cal using clean model
#
# Note: clean will have left FT of model in the MODEL_DATA column
# If you've done something in between, can use the ft task to
# do this manually.
#
print '--SelfCal 1--'
default('gaincal')

vis = srcsplitms

print "Will self-cal using MODEL_DATA left in MS by clean"

# New gain table
selfcaltab1 = srcsplitms + '.selfcal1'
caltable = selfcaltab1

print "Will write gain table "+selfcaltab1

# Don't need a-priori cals
selectdata = False
gaincurve = False
opacity = 0.0

# This choice seemed to work
refant = '11'

# Lets do phase-only first time around
gaintype = 'G'
calmode = 'p'

# Do scan-based solutions with SNR>3
solint = 0.0
minsnr = 3.0

# Do not need to normalize (let gains float)
solnorm = False

gaincal()

#
#---------------------------------------------------------------------
#
# Correct the data (no need for interpolation this stage)
#
print '--ApplyCal--'
default('applycal')

vis = srcsplitms

print "Will apply self-cal table to over-write CORRECTED_DATA in MS"

gaintable = selfcaltab1

gaincurve = False
opacity = 0.0
field = ''
spw = ''
selectdata = False

calwt = True

applycal()

# Self-cal is now in CORRECTED_DATA column of split ms
#
#=====================================================================
# SECOND CLEAN / SELFCAL CYCLE
#=====================================================================
#
print '--Clean 2--'
default('clean')

print "Now clean on self-calibrated data"

vis = srcsplitms

imname2 = prefix + '.clean2'
imagename = imname2

field = '*'
spw = ''
mode = 'mfs'
gain = 0.1
niter = 10000
threshold=0.04

alg = clnalg
imsize = clnimsize
cell = clncell

weighting = 'briggs'
rmode = 'norm'
robust = 0.5

cleanbox = 'interactive'
npercycle = 100

clean()

# Set up variables
clnimage2 = imname2+'.image'
clnmodel2 = imname2+'.model'
clnresid2 = imname2+'.residual'
clnmask2  = imname2+'.clean_interactive.mask'

print ""
print "----------------------------------------------------"
print "Clean"
print "Final clean model is "+clnmodel2
print "Final restored clean image is "+clnimage2
print "The clean residual image is "+clnresid2
print "Your final clean mask is "+clnmask2

print ""
print "This is the final restored clean image in the viewer"
print "Zoom in and set levels to see faint emission"
print "Use rectangle drawing tool to box off source"
print "Double-click inside to print statistics"
print "Move box on-source and get the max"
print "Calcualte DynRange = MAXon/RMSoff"
print "This time I got 1.076 / 0.001389 = 775 (better)"
print "Still not as good as it can be - lets selfcal again"
print "Close viewer panel when done"

#
#---------------------------------------------------------------------
#
# If you did not do interactive clean, bring up viewer manually
viewer(clnimage2,'image')

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

# jupiter6cm.usecase.clean2.image     (Jy/beam)
# 
# n           Std Dev     RMS         Mean        Variance    Sum
# 5236        0.001389    0.001390    3.244e-05   1.930e-06   0.1699    
# 
# Flux        Med |Dev|   IntQtlRng   Median      Min         Max
# 0.01060     0.0009064   0.001823    -1.884e-05  -0.004015   0.004892  
# 
# 
# On Jupiter:
# 
# n           Std Dev     RMS         Mean        Variance    Sum
# 5304        0.08512     0.08629     0.01418     0.007245    75.21     
# 
# Flux        Med |Dev|   IntQtlRng   Median      Min         Max
# 4.695       0.0008142   0.001657    0.0001557   -0.004526   1.076     
#
# Estimated dynamic range = 1.076 / 0.001389 = 775 (better)
#
# Note that the exact numbers you get will depend on how deep you
# take the interactive clean and how you draw the box for the stats.
#
#---------------------------------------------------------------------
#
# Next self-cal cycle
#
print '--SelfCal 2--'
default('gaincal')

print "Self-cal again using MODEL_DATA from second clean"

vis = srcsplitms

selfcaltab2 = srcsplitms + '.selfcal2'
caltable = selfcaltab2

print "Write self-cal table "+selfcaltab2

selectdata = False
gaincurve = False
opacity = 0.0
refant = '11'

# This time amp+phase on 10s timescales SNR>1
gaintype = 'G'
calmode = 'ap'
solint = 10.0
minsnr = 1.0
solnorm = False

gaincal()

#
# It is useful to put this up in plotcal
#
#---------------------------------------------------------------------
#
print '--PlotCal--'
default('plotcal')

caltable = selfcaltab2
multiplot = True
yaxis = 'amp'

plotcal()

# Use the Next button to iterate over antennas

print ""
print "-------------------------------------------------"
print "Plotcal"
print "Looking at amplitude in self-cal table "+caltable
print "Note coherence of solutions"

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

yaxis = 'phase'

plotcal()

#
# You can see it is not too noisy.

print ""
print "Now look at phases"
print "These are not too noisy"

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

# Lets do some smoothing anyway.
#
#---------------------------------------------------------------------
#
# Smooth calibration solutions
#
print '--Smooth--'
default('smoothcal')

print "Lets do smoothcal anyway as a demnstration"

vis = srcsplitms

tablein = selfcaltab2

smoothcaltab2 = srcsplitms + '.smoothcal2'
caltable = smoothcaltab2

print "Will write smoothed table "+smoothcaltab2

# Do a 30s boxcar average
smoothtype = 'mean'
smoothtime = 30.0

smoothcal()

# If you put into plotcal you'll see the results
# For example, you can grap the inputs from the last
# time you ran plotcal, set the new tablename, and plot!
#run plotcal.last
#tablein = smoothcaltab2
#plotcal()

#
#---------------------------------------------------------------------
#
# Correct the data
#
print '--ApplyCal--'
default('applycal')

print "Calibrate data using the smoothed table"

vis = srcsplitms

gaintable = smoothcaltab2

gaincurve = False
opacity = 0.0
field = ''
spw = ''
selectdata = False
calwt = True

applycal()

#
#=====================================================================
# THIRD CLEAN / SELFCAL CYCLE
#=====================================================================
#
print '--Clean 3--'
default('clean')

print "This is the 3rd clean using smoothed 2nd self-cal"

vis = srcsplitms

imname3 = prefix + '.clean3'
imagename = imname3

field = '*'
spw = ''
mode = 'mfs'
gain = 0.1
niter = 10000
threshold=0.04

alg = clnalg
imsize = clnimsize
cell = clncell

weighting = 'briggs'
rmode = 'norm'
robust = 0.5

cleanbox = 'interactive'
npercycle = 100

clean()

# Cleans alot deeper
# You can change the npercycle to larger numbers
# (like 250 or so) as you get deeper also.

# Set up variables
clnimage3 = imname3+'.image'
clnmodel3 = imname3+'.model'
clnresid3 = imname3+'.residual'
clnmask3  = imname3+'.clean_interactive.mask'

print ""
print "-------------------------------------------------"
print "Clean"
print "Final clean model is "+clnmodel3
print "Final restored clean image is "+clnimage3
print "The clean residual image is "+clnresid3
print "Your final clean mask is "+clnmask3

print ""
print "This is the final restored clean image in the viewer"
print "Zoom in and set levels to see faint emission"
print "Use rectangle drawing tool to box off source"
print "Double-click inside to print statistics"
print "Move box on-source and get the max"
print "Calcualte DynRange = MAXon/RMSoff"
print "I got 1.076 / 0.001015 = 1060 (even better!)"
print "Greg Taylor got 1600:1 - can you do better than me?"
print ""
print "Close viewer panel when done"

#
#---------------------------------------------------------------------
#
# If you did not do interactive clean, bring up viewer manually
viewer(clnimage3,'image')

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

# jupiter6cm.usecase.clean3.image     (Jy/beam)
# 
# n           Std Dev     RMS         Mean        Variance    Sum
# 5848        0.001015    0.001015    -4.036e-06  1.029e-06   -0.02360  
# 
# Flux        Med |Dev|   IntQtlRng   Median      Min         Max
# -0.001470   0.0006728   0.001347    8.245e-06   -0.003260   0.003542  
# 
# 
# On Jupiter:
# 
# n           Std Dev     RMS         Mean        Variance    Sum
# 6003        0.08012     0.08107     0.01245     0.006419    74.72     
# 
# Flux        Med |Dev|   IntQtlRng   Median      Min         Max
# 4.653       0.0006676   0.001383    -1.892e-06  -0.002842   1.076     
# 
# Estimated dynamic range = 1.076 / 0.001015 = 1060 (even better!)
#
# Note that the exact numbers you get will depend on how deep you
# take the interactive clean and how you draw the box for the stats.
#
# Greg Taylor got 1600:1 so we still have some ways to go
# This will probably take several more careful self-cal cycles.

# Set up final variables
clnimage = clnimage3
clnmodel = clnmodel3
clnresid = clnresid3
clnmask  = clnmask3

print ""
print "--------------------------------------------------"
print "After this script is done you can continue on with"
print "more self-cal, or try different cleaning options"

#
#=====================================================================
# Image Analysis
#=====================================================================
#
# Can do some image statistics if you wish
print '--Final Imstat--'
default('imstat')

imagename = clnimage
imstat()

on_statistics = xstat

# Now do stats in the lower right corner of the image
# remember clnimsize = [288,288]
box = '216,1,287,72'
imstat()

off_statistics = xstat

#=====================================================================
#
# Print results and regression versus previous runs
# WARNING: currently requires toolkit
#
print ""
print ' Jupiter results '
print ' =============== '
print ''
# Pull the max and rms from the clean image
thistest_immax=on_statistics['max'][0]
oldtest_immax = 1.07732224464
print ' Clean image ON-SRC max should be ',oldtest_immax
print ' Found : Max in image = ',thistest_immax
diff_immax = abs((oldtest_immax-thistest_immax)/oldtest_immax)
print ' Difference (fractional) = ',diff_immax

print ''
thistest_imrms=off_statistics['rms'][0]
oldtest_imrms = 0.0010449
print ' Clean image OFF-SRC rms should be ',oldtest_imrms
print ' Found : rms in image = ',thistest_imrms
diff_imrms = abs((oldtest_imrms-thistest_imrms)/oldtest_imrms)
print ' Difference (fractional) = ',diff_imrms

print ''
print ' Final Clean image Dynamic Range = ',thistest_immax/thistest_imrms
print ''
print ' =============== '

print ''
print '--- Done ---'

#
#=====================================================================