#######################################################################
# Imaging Script for SgrA* Band 6 ALMA data                           #
#######################################################################

"""


At the end of the Calibration script you should have SgrA_B6.ms
which contains corrections derived from the calibrators and flagging.
It has been averaged up every 30s and every 5 channels, this can be
changed at the end of the calibration script.

This measurement set can also be found here:
SgrA_Band6_CalibratedData.tgz

As of mid-Sept, the central SgrA* field is well-calibrated but the
other mosaic fields have delay errors on some baselines and some
discrepancies between XX and YY. Nonetheless, the expected noise (~5.5
mJy/beam) is achieved in line-only cubes made averaging 5 channels,
showing that the data set is basically good.

Lines are visible in the visibility spectrum the central field (4) on
shorter baselines, e.g. <65 m in spw 1 and 2, and in some of the outer
fields (to be taken with caution) on baselines <50 m.

No amplitude self-cal is performed in case of variability. There is a
small variation in amplitude with time, which might be astrophysical,
although some change is seen on shorter baselines. The phase-reference
calibrated data is flat with time. It might be worth doing amplitude
self-calibration.

The H30alpha line at 231.900094 GHz (in spw 2) is imaged; spw 2 and 3
are used for continuum subtraction.

The user-created variable contspwf0 in this script defines the
line-free channels used for continuum imaging and as a basis for
selection for continuum subtraction.  It was found by inspecting the
data for the central field, in plotms, change the selection if
desired.

If no prior channel averaging were done, this would be
'0:86~1630;1856~3839,1:1~300;1641~2865,2:1~1450;2466~3839,3'

Using the averaging every 5 channels this is
contspwf0='0:19~326;371~767,1:0~60;328~573,2:0~290;493~767,3'

For continuum subtraction around H30alpha:
contchans='2:0~280;540~767;3'

For imaging the line
linechans='2:280~540'

- if you change the channel averaging interval, change the above variables

"""

#------------------------#
#------------------------#
#------ Script Sections ------#
#-- use value T to activate --#

#-- User-defined variables for selection of channels  --#
# set after inspection in plotms on all and on short baselines
# make sure all 3 sets are consistent

# Continuum-only channels
contspwf0='0:19~326;371~767,1:0~60;328~573,2:0~290;493~767,3'

# channels to use as continuum around H30alpha
contchans='2:0~280;540~767;3'     # selection consistent with contspwf0

# H30alpha line channels
linechans='2:280~540' # selection consistent with contspwf0. contchans

#-- User-defined variables for imaging --#
# Interactive cleaning of continuum central field? Set to T or F

#intclean= T
intclean= F

# Mask (clean box) - this should be set if intclean = F
#cbox = ''
cbox = [65,65,192,192]

#- INITIAL CHECKS-#

test=F			# quick test to see if script looking in right
			# place			
do_listobs=F            # List imput measurement set

#- IMAGE AND PHASE-ONLY SELF-CAL CENTRAL FIELD CONTINUUM-#
# Using contspwf0 defined above
# Default is non-interactive clean, see intclean and cbox variables
# above for interactive clean

target_image0=F         # First image of target 
target_self_cal1=F      # 1st phase-only self-cal of continuum
			# central field, applied to all fields/channels
target_self_cal2=F      # 2nd self-cal with improved model
target_self_cal3=F      # 3rd self-cal with improved model

#- IMAGE CALIBRATED CONTINUUM MOSAIC -#			
cont_mosaic=F           # Make mosaic image of target (interactive= T)
# set mask yourself if you want non-interactive

#- IMAGE CALIBRATED CONTINUUM CENTRAL SPEC INDEX -#	

# Default is non- interactive clean, see intclean and cbox variables
# above for non-interactive clean

cont_spindex=T          # Spectral index clean central field

#- SUBTRACT CONTINUUM FROM H30alpha LINE, MAKE CUBE -#

# using contchans defined above
Halpha_contsub=T	# Subtract continuum from spw 2

# using linechans defined above
line_split=T            # Split out H30alpha line
H30a_cube=T             # Image the H30alpha cube and make moments
# Whole region above c. 20% mosaic PB is masked for non-interactive
# clean, change if you wish
                      
###################################

#------------------------#
#------------------------#
#--------test------------#

if (test):
	os.system('ls -lh uid*')

#------------------------#
#------------------------#
#-----List Obs Data------#

if (do_listobs):
	default(listobs)		#show MS information
	verbose= F
	vis='SgrA_all_B6.ms'
	listobs
		

# Sgr A* fields now numbered 0~6, 0 being centre.

#---------------------------------------#
#---------------------------------------#
#-- Chose line-free channels of SgrA* --#

# Use plotms, inspecting all baselines and  0~50 m baselines.
# Assume that flanking fields have lines in same places
# (this corresponds to a +/- 300 km/s span excluded around H30 alpha
# line, in spw 2). 

#-----------------------------------------------------#
#-----------------------------------------------------#
#--- First continuum image of target central field ---#

if (target_image0):
	print "<< Make first image of central field continuum\n"
	print "<< using interactive = ", intclean
	print "<< and (if noninteractive) mask = ",cbox
	default(clean)
	vis='SgrA_B6.ms'
	field='0'
	imagename = 'SgrA_f0-0.clean'
	psfmode='hogbom'
	interactive= intclean
	mask = cbox
	cell=['0.25arcsec']
	spw=contspwf0
	niter=5000
	threshold='2mJy'
        clean()

# I used elliptical boxes increasing with cycle for up to 5 x 1000 iter

        default(imstat)
	imagename  =  'SgrA_f0-0.clean.image'
	box= '21,21,70,230'
	sgrastats=imstat()
	sgrarms=sgrastats['rms']
	box= '51,51,200,200'
	sgrastats=imstat()
	sgramax=sgrastats['max']
	snr = sgramax/sgrarms
	print 'SgrA*  snr', snr[0], 'peak', sgramax[0], 'rms', sgrarms[0]
# 'SgrA* snr', 247, 'peak', 4.423, 'rms', 0.018

#-------------------------------------------#
#-------------------------------------------#
#--- Self-calibrate target central field ---#

# Amplitudes may be variable but position i.e. phase should not be, so
# I used several rounds of phase-only self-calibration of the
# continuum-selected channels of the central field, applying the
# corrections to all channels and fields.
# The cleans are interactive but you could save a mask and re-use it
# if you are convinced of what is real.
# Checking in plotcal and plotms not included in script

if (target_self_cal1):
	print "<< Self-calibrate and image central field cycle 1\n"
	default(gaincal)
	vis='SgrA_B6.ms'
	field='0'
	caltable='SgrA_f0-1.phcal'
	calmode='p'
	spw=contspwf0
	solint='inf'
	refant='8'
	gaincal()
	print "<< Generated SgrA_f0-1.phcal\n"
	
	default(applycal)
	vis='SgrA_B6.ms'
	gaintable='SgrA_f0-1.phcal'
	field='0~6'
	interp='nearest'   # to ensure all fields are covered
	applycal()
	print "<< Applied SgrA_f0-1.phcal to all fields\n"

	default(clean)
	vis='SgrA_B6.ms'
	field='0'
	imagename = 'SgrA_f0-1.clean'
	psfmode='hogbom'
	interactive= intclean
	mask = cbox
	cell=['0.25arcsec']
	spw=contspwf0
	niter=5000
	threshold='1mJy'
	clean()
	print "<< Created SgrA_f0-1.clean.image (central field)\n"

	default(imstat)
	imagename  = 'SgrA_f0-1.clean.image'
	box= '21,21,70,230'
	sgrastats=imstat()
	sgrarms=sgrastats['rms']
	box= '51,51,200,200'
	sgrastats=imstat()
	sgramax=sgrastats['max']
	snr = sgramax/sgrarms
	print 'SgrA*  snr', snr[0], 'peak', sgramax[0], 'rms', sgrarms[0]
#SgrA*  snr', 645, 'peak', 4.5111, 'rms', 0.0070

if (target_self_cal2):
	print "<< Self-calibrate and image central field cycle 2\n"
        default(gaincal)
	vis='SgrA_B6.ms'
	field='0'
	caltable='SgrA_f0-2.phcal'
	calmode='p'
	spw=contspwf0
	solint='inf'
	refant='8'
	gaincal()
	print "<< Generated SgrA_f0-2.phcal\n"
	
	default(applycal)
	vis='SgrA_B6.ms'
	gaintable='SgrA_f0-2.phcal'
	field='0~6'
	interp='nearest'
	applycal()
	print "<< Applied SgrA_f0-2.phcal to all fields\n"

# Try multiscale clean
# Largest spatial scale in the data probably <10"

        default(clean)
	vis='SgrA_B6.ms'
	field='0'
	imagename = 'SgrA_f0-2.clean'
	psfmode='hogbom'
	interactive= intclean
	mask = cbox
	cell=['0.25arcsec']
	spw=contspwf0
	niter=10000
	threshold='1mJy'
	multiscale=[0,5,10]
	negcomponent   = 100
        clean()
	print "<< Created SgrA_f0-2.clean.image (using multi-scale)\n"	
	
	default(imstat)
	imagename  =  'SgrA_f0-2.clean.image'
	box= '21,21,70,230'
	sgrastats=imstat()
	sgrarms=sgrastats['rms']
	box= '51,51,200,200'
	sgrastats=imstat()
	sgramax=sgrastats['max']
	snr = sgramax/sgrarms
	print 'SgrA*  snr', snr[0], 'peak', sgramax[0], 'rms', sgrarms[0]
# SNR is slightly improved, large scale flux looks good
# I got SgrA*  snr', 746, 'peak', 4.5004, 'rms', 0.0060

if (target_self_cal3):
	print "<< Self-calibrate and image central field cycle 3\n"
        default(gaincal)
	vis='SgrA_B6.ms'
	field='0'
	caltable='SgrA_f0-3.phcal'
	calmode='p'
	spw=contspwf0
	solint='inf'
	refant='8'
	gaincal()
	print "<< Generated SgrA_f0-3.phcal\n"
	
	default(applycal)
	vis='SgrA_B6.ms'
	gaintable='SgrA_f0-3.phcal'
	field='0~6'
	interp='nearest'
	applycal()
	print "<< Applied SgrA_f0-3.phcal to all fields\n"

        default(clean)
	vis='SgrA_B6.ms'
	field='0'
	imagename = 'SgrA_f0-3.clean'
	psfmode='hogbom'
	interactive= intclean
	mask = cbox
	cell=['0.25arcsec']
	spw=contspwf0
	niter=20000
	threshold='1mJy'
	multiscale=[0,5,10]
	negcomponent   = 100
        clean()
	print "<< Created SgrA_f0-3.clean.image (using multi-scale)\n"	
	
	default(imstat)
	imagename  =  'SgrA_f0-3.clean.image'
	box= '21,21,70,230'
	sgrastats=imstat()
	sgrarms=sgrastats['rms']
	box= '51,51,200,200'
	sgrastats=imstat()
	sgramax=sgrastats['max']
	snr = sgramax/sgrarms
	print 'SgrA*  snr', snr[0], 'peak', sgramax[0], 'rms', sgrarms[0]
#SgrA*  snr', 791, 'peak', 4.5014, 'rms', 0.0057
# Small improvement, I stopped here, but you could do further self-cal

#------------------------------#
#------------------------------#
#--- Image continuum mosaic ---#

# NB cannot also solve for spectral index, but outlying stuff is
# probably too faint for that to matter

if (cont_mosaic):
	print "<< Cleaning whole mosaic continuum, interactively\n"
	default(clean)
	vis='SgrA_B6.ms'
	field='0~6'
	imagename = 'SgrA_mosaic_cont.clean'
	psfmode='hogbom'
	interactive= T
	cell=['0.25arcsec']
	spw=contspwf0
	niter=20000
	npercycle=50  # clean peak first, then increase for diffuse emission
	threshold='2mJy'
	multiscale=[0,5,10]
	negcomponent   = 100
	imagermode='mosaic'
	cyclefactor = 1
	cyclespeedup=50
	phasecenter='0'
	pbcor= T
	clean()
	print "<< Created SgrA_mosaic_cont.clean.image\n"

#-------------------------------#
#-------------------------------#
#--- Image continuum spindex ---#

if (cont_spindex):
	print "<< Make spectral index image for central field \n"
	print "<< using interactive = ", intclean
	print "<< and (if noninteractive) mask = ",cbox
	default(clean)
	vis='SgrA_B6.ms'
	field='0'
	nterms=2
	imagename = 'SgrA_spindex.clean'
	psfmode='hogbom'
	interactive= intclean
	mask = cbox
	cell=['0.25arcsec']
	spw=contspwf0
	niter=500
	npercycle = 50     
	threshold='1mJy'
	multiscale=[0,5,10]
	negcomponent   = 100
	clean()
	print "<< Created SgrA_spindex.clean.image.tt0 (total intensity)\n"
	print "<< and SgrA_spindex.clean.image.alpha (spectral index)\n"

#----------------------------------------#
#----------------------------------------#
#--- Subtract continuum in Halpha spw ---#
# Only use spw 2~3 as 0~1 have other lines and are far away in frequency
# Use default per-integration (now 30s) solint

if (Halpha_contsub):
	print "<< Subtract continuum in spw 2, 3, from the H30a line\n"
	print "<< Using continuum channels ", contchans
	default(uvcontsub)
	vis='SgrA_B6.ms'
	field  = ''
	fitspw  = contchans
	combine =  'spw'
	fitorder =1
	uvcontsub()
	print "<< Created SgrA_B6.ms.contsub\n"

# This produces SgrA_B6.ms.contsub
# You could check in plotms

#-----------------------#
#-----------------------#
#--- Split out line  ---#

if (line_split):
	print "<< Split out line channels ", linechans
	default(split)
	vis='SgrA_B6.ms.contsub'
	datacolumn='data'
	field=''
	outputvis='SgrA_B6_fall_H30a.ms'
	spw=linechans
	split()
	print "<< Created SgrA_B6_fall_H30a.ms\n"

#-----------------------------------#
#-----------------------------------#
#--- Make H30alpha spectral cube ---#

# Make line cube.  I made a test cube of the central field and an
# integrated flux (0th moment) map first, then used the 0th moment map
# to guide setting a mask for the whole mosaic cube.  It might be
# better to mask separately for different velocity ranges.
# However, the initial map made by including everything down to 20% PB
# was not very different.
# From sensitivity calculator, if two fields contribute to each
# position, expect noise of about 5.5 mJy/beam which is exactly what
# we get in quiet channels (varies from 5-6 mJy/beam/channel, more
# near bright peaks).

if (H30a_cube):
	print "<< About to clean cube interactively \n"
	default(clean)
	vis='SgrA_B6_fall_H30a.ms'
	field='0~6'
	imagename = 'SgrA_fall_H30a.clean'
	psfmode='hogbom'
	interactive= T
	cell=['0.25arcsec']
	niter=50000
	threshold='5mJy'
	mode='channel'
	multiscale=[0,5,10]
	negcomponent   = 100
	restfreq='231.900094GHz'
	outframe='LSRK'
	imagermode='mosaic'
	phasecenter='0'
	pbcor= T
	clean()
	print "<< Created SgrA_fall_H30a.clean.image \n"

	default(immoments)
	imagename = 'SgrA_fall_H30a.clean.image'
	chans='31~230'
	moments=[0,1]
	excludepix=[-0.03,0.03]
	outfile='SgrA_fall_H30a.mom'
	immoments()
	print "<< Created total intensity and velocity images\n"
	print "<< SgrA_fall_H30a.mom.intergrated, SgrA_fall_H30a.mom.weighted_coord\n"


	exportfits(imagename='SgrA_fall_H30a.mom.weighted_coord',
		   fitsimage='SgrA_fall_H30a.mom.weighted_coord.fits')
	exportfits(imagename='SgrA_fall_H30a.mom.integrated',
		   fitsimage='SgrA_fall_H30a.mom.integrated.fits')
	exportfits(imagename='SgrA_mosaic_cont.clean.image',
		   fitsimage='SgrA_mosaic_cont.clean.fits')
	exportfits(imagename='SgrA_fall_H30a.clean.image',
		   fitsimage='SgrA_fall_H30a.clean.fits')