# This script describes the continuum imaging and self-calibration for
# the *COMBINED* ALMA LBC SV B3, 2.9mm continuum data on HL Tau.

# requires CASA 4.2.2 or higher. Meant to be run interactively.

############################################################
# Concat (you will need to make the paths agree with your
# data locations). 


os.system('rm -rf HLTau_B3cont16.ms')
concat(vis=['../HLTAU_B3_HCOp/HLTau_B3HCN_1_7.fixed.calavg16',
            '../HLTAU_B3_CO/HLTau_B3CO_1_7.fixed.calavg16'],
       concatvis='HLTau_B3cont16.ms')

################################################################
################################################################
# THIS IS THE SCIENCE PORTAL (SP) DATA PRODUCT CALLED "CalibratedData"
# If you have downloaded this product start below this section.
# To go directly to self-calibrated uv-data, run all flagdata
# commands below and using the calibration tables on the SP run
# the final applycal. 
################################################################
################################################################

# Inspect data
os.system('rm -rf HLTau_B3cont16.ms.listobs')
listobs(vis='HLTau_B3cont16.ms',listfile='HLTau_B3cont16.ms.listobs')

os.system('rm -rf HLTau_B3cont16.ms.ampuvdist.png')
plotms(vis='HLTau_B3cont16.ms',coloraxis='spw',yaxis='amp',
       xaxis='uvdist',avgchannel='16',avgtime='100s',avgscan=False,
       plotfile='HLTau_B3cont16.ms.ampuvdist.png')

##############################################################
# Continuum Imaging

# In the combined continuum data, there are 4 spws clustered around
# 102 GHz and 1 at 90 GHz. This spread requires nterms=2 due to the
# significant changes in spectral index across the source.
# This was attempted, but led to a poor spectral index result --
# The lowest frequency spw, with much less data than the combined
# cluster at 102 GHz simply does not have enough sensitivity to measure the
# spectral index in the out portions of the disk. Thus, to speed imaging,
# we limit the cleans below to spw='1~4' and use nterms=1. 

# There are 2 weak continuum sources in the outer parts of the primary
# beam.  To save time we do the self-calibration on just the HL Tau
# region, and only at the end make a large image. Note that this is
# not advisable if the outlying sources are a significant fraction of
# the total flux (without primary beam correction). In that case one
# can setup outlier fields or make the big images.

# This script is set up assuming it will be run interactively. The number of
# iterations given for the progression of interactive clean is for guidance
# and need not be followed exactly but is strongly recommended that care
# is taken to set the clean mask.

os.system('rm -rf HLTau_B3.contms.*')
clean(vis='HLTau_B3cont16.ms',
      imagename='HLTau_B3.contms',
      spw='1~4',
      mode='mfs',nterms=1,imagermode='csclean',
      imsize=640,cell='0.01arcsec',
      multiscale=[0,5,15],negcomponent=10,
      interactive=True,
      mask='',
      weighting='briggs',robust=0.0,
      niter=10000,threshold='0.1mJy')

# Iterations in interactive 100, 500, 1000, 3000

# Set time to get about several scans per solution

os.system('rm -rf b3cont_pcal1')
gaincal(vis='HLTau_B3cont16.ms',caltable='b3cont_pcal1',gaintype='T',
        refant='DV10',calmode='p',combine='spw,scan',spw='1~4',
        solint='300s',minsnr=3.0,minblperant=6)


plotcal(caltable='b3cont_pcal1',xaxis='time',yaxis='phase',
        timerange='<2014/10/15/00:00:00',
        iteration='antenna',subplot=421,plotrange=[0,0,-180,180])

plotcal(caltable='b3cont_pcal1',xaxis='time',yaxis='phase',
        timerange='2014/10/15/00:00:00~2014/10/31/00:00:00',
        iteration='antenna',subplot=421,plotrange=[0,0,-180,180])

plotcal(caltable='b3cont_pcal1',xaxis='time',yaxis='phase',
        timerange='2014/11/01/00:00:00~2014/11/03/00:00:00',
        iteration='antenna',subplot=421,plotrange=[0,0,-180,180])

plotcal(caltable='b3cont_pcal1',xaxis='time',yaxis='phase',
        timerange='2014/11/03/00:00:0~2014/11/12/00:00:00',
        iteration='antenna',subplot=421,plotrange=[0,0,-180,180])

plotcal(caltable='b3cont_pcal1',xaxis='time',yaxis='phase',
        timerange='>2014/11/12/00:00:00',
        iteration='antenna',subplot=421,plotrange=[0,0,-180,180])

# Note that spwmap must be set to the lowest spw number in the combine='spw'
# gaincal, and there must be an entry for every spw.
applycal(vis='HLTau_B3cont16.ms',spwmap=[1,1,1,1,1],
         gaintable=['b3cont_pcal1'],gainfield='',calwt=F,flagbackup=F)

os.system('rm -rf HLTau_B3.contms_p1.*')
clean(vis='HLTau_B3cont16.ms',
      imagename='HLTau_B3.contms_p1',
      spw='1~4',
      mode='mfs',nterms=1,imagermode='csclean',
      imsize=640,cell='0.01arcsec',
      multiscale=[0,5,15],negcomponent=10,
      interactive=True,
      mask='HLTau_B3.contms.mask',
      weighting='briggs',robust=0.0,
      niter=10000,threshold='0.03mJy')

# Iterations in interactive 1000, 3000


os.system('rm -rf b3cont_pcal2')
gaincal(vis='HLTau_B3cont16.ms',caltable='b3cont_pcal2',gaintype='T',
        refant='DV10',calmode='p',combine='spw',spw='1~4',
        solint='inf',minsnr=3.0,minblperant=6)


plotcal(caltable='b3cont_pcal2',xaxis='time',yaxis='phase',
        timerange='<2014/10/15/00:00:00',
        iteration='antenna',subplot=421,plotrange=[0,0,-180,180])

plotcal(caltable='b3cont_pcal2',xaxis='time',yaxis='phase',
        timerange='2014/10/15/00:00:00~2014/10/31/00:00:00',
        iteration='antenna',subplot=421,plotrange=[0,0,-180,180])

plotcal(caltable='b3cont_pcal2',xaxis='time',yaxis='phase',
        timerange='2014/11/01/00:00:00~2014/11/03/00:00:00',
        iteration='antenna',subplot=421,plotrange=[0,0,-180,180])

plotcal(caltable='b3cont_pcal2',xaxis='time',yaxis='phase',
        timerange='2014/11/03/00:00:0~2014/11/12/00:00:00',
        iteration='antenna',subplot=421,plotrange=[0,0,-180,180])

plotcal(caltable='b3cont_pcal2',xaxis='time',yaxis='phase',
        timerange='>2014/11/12/00:00:00',
        iteration='antenna',subplot=421,plotrange=[0,0,-180,180])


# *Optional* If you are concerned that too much data will be flagged
# by applying this self-cal you can back it up to make it easier to go
# back.
os.system('rm -rf HLTau_B3cont16.ms.beforep2')
os.system('cp -r HLTau_B3cont16.ms HLTau_B3cont16.ms.beforep2')

applycal(vis='HLTau_B3cont16.ms',spwmap=[1,1,1,1,1],
         gaintable=['b3cont_pcal2'],calwt=F,flagbackup=T)

# Try going directly to threshold this time (interactive=False), but then
# inspect residual. If more cleaning needed restart WITHOUT including
# os.system command
os.system('rm -rf HLTau_B3.contms_p2.*')
clean(vis='HLTau_B3cont16.ms',
      imagename='HLTau_B3.contms_p2',
      spw='1~4',field='',
      mode='mfs',nterms=1,imagermode='csclean',
      imsize=640,cell='0.01arcsec',
      multiscale=[0,5,15],negcomponent=10,
      interactive=False,
      mask='HLTau_B3.contms_p1.mask',
      weighting='briggs',robust=0.0,
      niter=10000,threshold='0.03mJy')

# this shows one could go a bit deeper, change threshold to 0.25 and
# restart with interactive. Note must remove the mask in the command or
# you will not be able to add new ones -- this forces it to use one on
# disk with name _p2.

clean(vis='HLTau_B3cont16.ms',
      imagename='HLTau_B3.contms_p2',
      spw='1~4',field='',
      mode='mfs',nterms=1,imagermode='csclean',
      imsize=640,cell='0.01arcsec',
      multiscale=[0,5,15],negcomponent=10,
      interactive=True,
      mask='',
      weighting='briggs',robust=0.0,
      niter=10000,threshold='0.025mJy')

# Ignore "model image on disk" warning, adjust mask if needed and clean
# 1000 iterations

# Again optional backup
os.system('rm -rf HLTau_B3cont16.ms.beforeap')
os.system('cp -r HLTau_B3cont16.ms HLTau_B3cont16.ms.beforeap')

# Get one solution per dataset
os.system('rm -rf b3cont_apcal')
gaincal(vis='HLTau_B3cont16.ms',caltable='b3cont_apcal',gaintype='T',
        refant='DV10',calmode='a',combine='spw,scan',spw='1~4',
        gaintable='b3cont_pcal2',spwmap=[1,1,1,1,1],
        solint='inf',minsnr=4.0,minblperant=6)

plotcal(caltable='b3cont_apcal',xaxis='time',yaxis='amp',timerange='',
        spw='1~4',subplot=111)

# Large gains for DA55 and DV07 (longest baselines), see what effect is
# on data (i.e. check for outliers after applying.

applycal(vis='HLTau_B3cont16.ms',spwmap=[[1,1,1,1,1],[1,1,1,1,1]],
         gaintable=['b3cont_pcal2','b3cont_apcal'],calwt=F,flagbackup=T)

plotms(vis='HLTau_B3cont16.ms',coloraxis='spw',yaxis='amp',avgchannel='16',
       xaxis='uvdist',avgtime='100',ydatacolumn='corrected')

# Flag a couple of prominent scans
flagdata(vis='HLTau_B3cont16.ms',scan='533,611')

os.system('rm -rf HLTau_B3.contms_ap.*')
clean(vis='HLTau_B3cont16.ms',
      imagename='HLTau_B3.contms_ap',
      spw='1~4',field='',
      mode='mfs',nterms=1,imagermode='csclean',
      imsize=640,cell='0.01arcsec',
      multiscale=[0,5,15],negcomponent=10,
      interactive=False,
      mask='HLTau_B3.contms_p2.mask',
      weighting='briggs',robust=0.0,
      niter=10000,threshold='0.02mJy')

########################################################################
# Make a biggger image covering whole primary beam

os.system('rm -rf HLTau_B3.contms_ap_big.*')
clean(vis='HLTau_B3cont16.ms',
      imagename='HLTau_B3.contms_ap_big',
      spw='1~4',
      mode='mfs',nterms=1,imagermode='csclean',
      imsize=8000,cell='0.01arcsec',
      multiscale=[0,5,15],negcomponent=10,
      interactive=True,
      mask='',
      weighting='briggs',robust=0.0,
      niter=10000,threshold='0.02mJy')

# Do primary beam correction

os.system('rm -rf HLTau_B3.contms_ap_big.image.pbcor')
impbcor(imagename='HLTau_B3.contms_ap_big.image',
        pbimage='HLTau_B3.contms_ap_big.flux',
        outfile='HLTau_B3.contms_ap_big.image.pbcor',
        cutoff=0.3)

exportfits('HLTau_B3.contms_ap_big.image.pbcor',
           'HLTau_B3.contms_ap_big.image.pbcor.fits')

########################################################################
# Make tappered image for comparison with tapered line data.
# Since the beam is not much bigger than the disk in this case,
# not using multiscale.

os.system('rm -rf HLTau_B3.cont_ap_200kl.*')
clean(vis='HLTau_B3cont16.ms',
      imagename='HLTau_B3.cont_ap_200kl',
      spw='1~4',field='',
      mode='mfs',nterms=1,imagermode='csclean',
      imsize=600,cell='0.1arcsec',
      interactive=True,
      mask='',
      uvtaper=True,outertaper='200klambda',
      weighting='briggs',robust=0.5,
      niter=10000,threshold='0.02mJy')