# This script describes the line imaging for the ALMA LBC SV B3 data
# on HL Tau that includes CO and CN (1-0). The combined B3
# imaging is described in a separate script.

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

#########################################################
# If working from the individual datasets
# Concat and split off continuum and line spws

os.system('rm -rf HLTau_B3CO_1_7.fixed.cal')
concat(vis=['X3b2/HLTau_X3b2.ms.split.fixed.cal',
            'X5fa/HLTau_X5fa.ms.split.fixed.cal',
            'X845/HLTau_X845.ms.split.fixed.cal',
            'X220/HLTau_X220.ms.split.fixed.cal',
            'X5b2/HLTau_X5b2.ms.split.fixed.cal',
            'X693/HLTau_X693.ms.split.fixed.cal',
            'X306/HLTau_X306.ms.split.fixed.cal'],
       concatvis='HLTau_B3CO_1_7.fixed.cal')

os.system('rm -rf HLTau_B3CO_1_7.fixed.cal.listobs')
listobs(vis='HLTau_B3CO_1_7.fixed.cal',listfile='HLTau_B3CO_1_7.fixed.cal.listobs')

# Note more continum channel averaging will cause bandwidth smearing
# when imaging the full primary beam
os.system('rm -rf HLTau_B3CO_1_7.fixed.calavg16')
split(vis='HLTau_B3CO_1_7.fixed.cal',
      outputvis='HLTau_B3CO_1_7.fixed.calavg16',
      width='16',datacolumn='data',spw='0~1')

os.system('rm -rf HLTau_B3_7exec.ms.cal.co')
split(vis='HLTau_B3CO_1_7.fixed.cal',spw='2',datacolumn='data',
      outputvis='HLTau_B3_7exec.ms.cal.co')

os.system('rm -rf HLTau_B3_7exec.ms.cal.cn')
split(vis='HLTau_B3CO_1_7.fixed.cal',spw='3',datacolumn='data',
      outputvis='HLTau_B3_7exec.ms.cal.cn')

################################################################
################################################################
# THESE ARE THE SCIENCE PORTAL (SP) DATA PRODUCT CALLED
# "CalibratedData" FOR THE LINE DATA. If you have downloaded these
# products start below this section for Line imaging. The combined B3
# imaging is described in a separate script.
################################################################
################################################################

# amplitude vs uv-distance, time-averaging within 5min

plotms(vis='HLTau_B3CO_1_7.fixed.calavg16',coloraxis='spw',yaxis='amp',
       xaxis='uvdist',avgchannel='16',avgtime='300s',avgscan=False)


#################################################### #   
# UVCONTSUB

plotms(vis='HLTau_B3HCN_1_7.fixed.calavg16',coloraxis='spw',yaxis='amp',
       xaxis='uvdist',avgchannel='4',avgtime='300s',avgscan=True)

plotms(vis='HLTau_B3_7exec.ms.cal.co',yaxis='amp',xaxis='channel',
       avgtime='1e8',avgscan=True)

uvcontsub(vis='HLTau_B3_7exec.ms.cal.cn',
          fitorder=1,fitspw='0:0~1700;2500~3839')

uvcontsub(vis='HLTau_B3_7exec.ms.cal.co',
          fitorder=1,fitspw='0:0~1700;2500~3839')

plotms(vis='HLTau_B3_7exec.ms.cal.co.contsub',yaxis='amp',xaxis='channel',
       avgtime='1e8',avgscan=True)

plotms(vis='HLTau_B3_7exec.ms.cal.co.contsub.split',yaxis='amp',xaxis='velocity',
       avgtime='1e8',avgscan=True)


##########################################################
# CO Line Imaging

# Reduce size of data for faster imaging
os.system('rm -rf HLTau_B3_7exec.ms.cal.co.contsub.split')
split(vis='HLTau_B3_7exec.ms.cal.co.contsub',
      outputvis='HLTau_B3_7exec.ms.cal.co.contsub.split',
      spw='0:1900~2300',datacolumn='data')

# Flag some bad data
flagdata(vis='HLTau_B3_7exec.ms.cal.co.contsub.split',
         observation='0,1',antenna='DV04')

flagdata(vis='HLTau_B3_7exec.ms.cal.co.contsub.split',
         observation='3,4,6',antenna='DV22')

name='HLTau_CO1_0.tap200klR0.5_0.25kms_ms'
os.system('rm -rf '+name+'.*')
# restart
clean('HLTau_B3_7exec.ms.cal.co.contsub.split',
      imagename=name,
      imagermode='csclean',weighting='briggs',robust=0.5,
      imsize=500, cell='0.2arcsec',multiscale=[0,5,15,45],
      negcomponent=10,
      mode='velocity',start='-6.0km/s',width='0.25km/s',nchan=136,
      outframe='LSRK',restfreq='115.27120GHz',
      interactive=True,
      mask='',
      uvtaper=True,outertaper='200klambda',
      niter=100000,threshold='30.0mJy',usescratch=True)

name='HLTau_CO1_0.tap200klR0.5_0.25kms_ms'
exportfits(name+'.image',name+'.image.fits')
exportfits(name+'.mask',name+'.mask.fits')
exportfits(name+'.flux',name+'.flux.fits')



############################################################
# Image CN

# Reduce size of data for faster imaging
os.system('rm -rf HLTau_B3_7exec.ms.cal.cn.contsub.split')
split(vis='HLTau_B3_7exec.ms.cal.cn.contsub',
      outputvis='HLTau_B3_7exec.ms.cal.cn.contsub.split',
      spw='0:1371~2680',datacolumn='data')

plotms(vis='HLTau_B3_7exec.ms.cal.cn.contsub.split',yaxis='amp',xaxis='velocity',
       avgtime='1e8',avgscan=True)


name='HLTau_CN1_0.tap200klR0.5_0.25kms_113.49097'
os.system('rm -rf '+name+'.*')
# restart
clean('HLTau_B3_7exec.ms.cal.cn.contsub.split',
      imagename=name,
      imagermode='csclean',weighting='briggs',robust=0.5,
      imsize=500, cell='0.2arcsec',
      negcomponent=10,
      mode='velocity',start='2.5km/s',width='0.25km/s',nchan=32,
      outframe='LSRK',restfreq='113.49097GHz',
      interactive=True,
      mask='',
      uvtaper=True,outertaper='200klambda',
      niter=10000,threshold='15.0mJy',usescratch=True)


name='HLTau_CN1_0.tap200klR0.5_0.25kms_113.48812'
os.system('rm -rf '+name+'.*')
# restart
clean('HLTau_B3_7exec.ms.cal.cn.contsub.split',
      imagename=name,
      imagermode='csclean',weighting='briggs',robust=0.5,
      imsize=500, cell='0.2arcsec',
      negcomponent=10,
      mode='velocity',start='2.5km/s',width='0.25km/s',nchan=32,
      outframe='LSRK',restfreq='113.48812GHz',
      interactive=False,
      mask='HLTau_CN1_0.tap200klR0.5_0.25kms_113.49097.mask',
      uvtaper=True,outertaper='200klambda',
      niter=10000,threshold='15.0mJy',usescratch=True)

name='HLTau_CN1_0.tap200klR0.5_0.25kms_113.49964'
os.system('rm -rf '+name+'.*')
# restart
clean('HLTau_B3_7exec.ms.cal.cn.contsub.split',
      imagename=name,
      imagermode='csclean',weighting='briggs',robust=0.5,
      imsize=500, cell='0.2arcsec',
      negcomponent=10,
      mode='velocity',start='2.5km/s',width='0.25km/s',nchan=32,
      outframe='LSRK',restfreq='113.49964GHz',
      interactive=False,
      mask='HLTau_CN1_0.tap200klR0.5_0.25kms_113.49097.mask',
      uvtaper=True,outertaper='200klambda',
      niter=10000,threshold='15.0mJy',usescratch=True)

name='HLTau_CN1_0.tap200klR0.5_0.25kms_113.50891'
os.system('rm -rf '+name+'.*')
# restart
clean('HLTau_B3_7exec.ms.cal.cn.contsub.split',
      imagename=name,
      imagermode='csclean',weighting='briggs',robust=0.5,
      imsize=500, cell='0.2arcsec',
      negcomponent=10,
      
      outframe='LSRK',restfreq='113.50891GHz',
      interactive=False,
      mask='HLTau_CN1_0.tap200klR0.5_0.25kms_113.49097.mask',
      uvtaper=True,outertaper='200klambda',
      niter=10000,threshold='15.0mJy',usescratch=True)

image='HLTau_CN1_0.tap200klR0.5_0.25kms_'
lines=['113.49097','113.48812','113.49964','113.50891']

for i in lines:
    exportfits(image+i+'.image',image+i+'.image.fits')

exportfits(image+'113.49097.mask',image+'113.49097.mask.fits')

#########################################################
# Example Moment Maps

# channel rms 8 mJy/beam line-free, considerably higher in
# line-full channels due to dyanmic range limitations.

# 1st blue channel ch24 0.0 km/s
# Last red channel ch104 20 km/s

# cut for systemic blue < ch 48, 6 km/s
# cut for systemic red   >ch 53, 7.25 km/s

# Moment zero all signal channels

name='HLTau_CO1_0.tap200klR0.5_0.25kms_ms.image'
rms=0.008

os.system('rm -rf '+name+'.mom0')
immoments(imagename=name,outfile=name+'.mom0',
          chans='24~104',moments=[0])


os.system('rm -rf '+name+'.mom0')
immoments(imagename=name,outfile=name+'.mom0mask',
          chans='24~104',moments=[0],
          includepix=[rms*1.0,1000])


# Set masked regions to zero if needed for some plotting purposes.
ia.open(name+'.mom0bluemask')
ia.replacemaskedpixels(0.0,update=True)
ia.close

ia.open(name+'.mom0redmask')
ia.replacemaskedpixels(0.0,update=True)
ia.close

# How to primary beam correct moment map if desired, though this
# mostly makes a big mess because the emission extends far outside the
# sensitive part of FOV.

imsubimage(imagename='HLTau_CO1_0.tap200klR0.5_0.25kms_ms.flux',
           outfile='HLTau_CO1_0.tap200klR0.5_0.25kms_ms.flux_1ch',
           chans='70')

impbcor(imagename=name+'.mom0mask',
        pbimage='HLTau_CO1_0.tap200klR0.5_0.25kms_ms.flux_1ch',
        outfile=name+'.mom0mask.pbcor',
        cutoff=0.25)

name='HLTau_CO1_0.tap200klR0.5_0.25kms_ms.image'
exportfits(name+'.mom0mask',name+'.mom0mask.fits')


##########################################################
# Moment map for strongest CN hyperfine

name='HLTau_CN1_0.tap200klR0.5_0.25kms_113.49097.image'
os.system('rm -rf '+name+'.mom0')
immoments(imagename=name,outfile=name+'.mom0',
          chans='12~20',moments=[0])

exportfits(name+'.mom0',name+'.mom0.fits')