# This script describes the combined imaging of the ALMA LBC SV data
# for Band 6 + Band 7

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

# It requires that you have either downloaded the Band 6 and 7 
# CalibratedData from the ALMA Science Portal SV page and applied 
# the flagdata commands in their individual imaging scripts along 
# with the self-calibration tables supplied at the Science Portal.

###########################################################
# Make an nterms=2 version that takes into account the
# changing spectral index in image. Note this adds a signifcant
# free parameter to the imaging and should only be attempted in
# presence of high S/N on at all sampled frequencies. If the
# spectral index map comes out mostly flagged, the nterms=2 is not
# well constrained.

# We have found nterms=2 with multiscale to require careful interactive 
# cleaning, in order to force major cycles more frequently than the 
# code would do automatically. It is recommended to clean slowly, 
# with about 1000 iterations per interactive cycle. If you have already 
# produced it, the Band 6 clean mask can be used. If you don't have it,
# comment out that line. 

os.system('rm -rf HLTau_B6B7cont_mscale_n2_ap.*')
clean(vis=['HLTau_B6cont.calavg','HLTau_B7cont.calavg'],
      imagename='HLTau_B6B7cont_mscale_n2_ap',
      mode='mfs',nterms=2,imagermode='csclean',
      imsize=2048,cell='0.0035arcsec',
      multiscale=[0,5,15],negcomponent=10,
      mask='HLTau_B6cont_mscale_ap.mask',
      weighting='briggs',robust=0.0,
      interactive=False,niter=1000000,threshold='0.03mJy')

# This mode produces more image products than others, the primary
# image at the mean combined frequency is called .tt0. The spectral
# index .alpha, and .alpha.error are also created. The .alpha map is
# automatically masked at the 4 sigma level.


exportfits('HLTau_B6B7cont_mscale_n2_ap.image.tt0',
           'HLTau_B6B7cont_mscale_n2_ap.image.tt0.fits')

exportfits('HLTau_B6B7cont_mscale_n2_ap.image.alpha',
           'HLTau_B6B7cont_mscale_n2_ap.image.alpha.fits')


exportfits('HLTau_B6B7cont_mscale_n2_ap.image.alpha.error',
           'HLTau_B6B7cont_mscale_n2_ap.image.alpha.error.fits')