# This script describes the line imaging for the ALMA LBC SV B3 data
# on HL Tau that includes HCO+ and HCN (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_B3HCN_1_7.fixed.cal')
concat(vis=['../Xc7f_1/HLTau_Xc7f.ms.split.fixed.cal',
            '../Xfd2_2/HLTau_Xfd2.ms.split.fixed.cal',
            '../X1f4_3/HLTau_X1f4.ms.split.fixed.cal',
            '../X444_4/HLTau_X444.ms.split.fixed.cal',
            '../Xc_5/HLTau_Xc.ms.split.fixed.cal',
            '../X22f_6/HLTau_X22f.ms.split.fixed.cal',
            '../X461_7/HLTau_X461.ms.split.fixed.cal'],
       concatvis='HLTau_B3HCN_1_7.fixed.cal')

os.system('rm -rf HLTau_B3HCN_1_7.fixed.cal.listobs')
listobs(vis='HLTau_B3HCN_1_7.fixed.cal',
        listfile='HLTau_B3HCN_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_B3HCN_1_7.fixed.calavg16')
split(vis='HLTau_B3HCN_1_7.fixed.cal',
      outputvis='HLTau_B3HCN_1_7.fixed.calavg16',
      width='16',datacolumn='data',spw='0~2')

os.system('rm -rf HLTau_B3HCN_1_7.fixed.calHCN')
split(vis='HLTau_B3HCN_1_7.fixed.cal',
      outputvis='HLTau_B3HCN_1_7.fixed.calHCN',
      width='',datacolumn='data',spw='3')

os.system('rm -rf HLTau_B3HCN_1_7.fixed.calHCOp')
split(vis='HLTau_B3HCN_1_7.fixed.cal',
      outputvis='HLTau_B3HCN_1_7.fixed.calHCOp',
      width='',datacolumn='data',spw='4')

################################################################
################################################################
# 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_B3HCN_1_7.fixed.calavg16',coloraxis='spw',yaxis='amp',
       xaxis='uvdist',avgchannel='16',avgtime='300s',avgscan=True)

####################################################
# Extra flagging based on plotms

data=['HLTau_B3HCN_1_7.fixed.calavg16',
      'HLTau_B3HCN_1_7.fixed.calHCN','HLTau_B3HCN_1_7.fixed.calHCOp']

for i in data:
    flagdata(i,mode='manual',scan='299,181')  
    flagdata(i,mode='manual',
         observation='4',antenna='DA49&PM04',flagbackup=False)
    # remove one stranded longest baseline
    flagdata(i,mode='manual',antenna='DV07&DA55',flagbackup=False)

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

plotms(vis='HLTau_B3HCN_1_7.fixed.calHCN',yaxis='amp',xaxis='frequency',
       iteraxis='spw',avgtime='1e8',avgscan=True)

plotms(vis='HLTau_B3HCN_1_7.fixed.calHCOp',yaxis='amp',xaxis='frequency',
       iteraxis='spw',avgtime='1e8',avgscan=True)

# HCOp restfreq='89.18853GHz'
# HCN 1-0
# F=1-1 88.63042 offset +4.837 km/s
# F=2-1 88.63185
# F=0-1 88.63394 offset -7.069 km/s
# channel width 61.035 kHz ~0.206 km/s

plotms(vis='HLTau_B3HCN_1_7.fixed.calHCOp',yaxis='amp',xaxis='velocity',
       spw='',avgtime='1e8',avgscan=True,
       transform=True,restfreq='89.18853GHz',freqframe='LSRK')

plotms(vis='HLTau_B3HCN_1_7.fixed.calHCN',yaxis='amp',xaxis='velocity',
       spw='',avgtime='1e8',avgscan=True,
       transform=True,restfreq='88.63185GHz',freqframe='LSRK')


uvcontsub(vis='HLTau_B3HCN_1_7.fixed.calHCN',
          fitspw='0:10~700;1200~1900',fitorder=1,
          spw='0:700~1200')


uvcontsub(vis='HLTau_B3HCN_1_7.fixed.calHCOp',
          fitspw='0:10~700;1200~1900',fitorder=1,
          spw='0:700~1200')

# Note in thes uvcontsubs we limit the output to the
# central channels as this speeds up the imaging. However,
# it is important to stay well away from the edges of the
# cubes you want to make when using this option or the
# Doppler regridding can run into problems.


#####################################################
# Line Imaging
# The 200 klambda taper used for the CO/CN data makes a
# very elliptical beam. Here we use a the full uvtaper
# specification to make it more circular.


name='HLTau_HCOp1_0.tap200klR0.5_0.25kms_ms'
os.system('rm -rf '+name+'.*')
# restart here
clean('HLTau_B3HCN_1_7.fixed.calHCOp.contsub',
      imagename=name,
      imagermode='csclean',weighting='briggs',robust=0.5,
      imsize=500, cell='0.2arcsec',multiscale=[0,5,15],
      negcomponent=10,
      mode='velocity',start='-6.0km/s',width='0.25km/s',nchan=136,
      outframe='LSRK',restfreq='89.18853GHz',
      interactive=False,
      mask='',
      uvtaper=True,outertaper=['200klambda','175klambda','45.0deg'],
      niter=10000,threshold='10.0mJy',usescratch=True)


name='HLTau_HCN1_0.tap200klR0.5_0.25kms'
os.system('rm -rf '+name+'.*')
# restart here
clean('HLTau_B3HCN_1_7.fixed.calHCN.contsub',
      imagename=name,
      imagermode='csclean',weighting='briggs',robust=0.5,
      imsize=500, cell='0.2arcsec',
      negcomponent=10,
      mode='velocity',start='-6.0km/s',width='0.25km/s',nchan=136,
      outframe='LSRK',restfreq='88.63185GHz',
      interactive=False,
      mask='',
      uvtaper=True,outertaper=['200klambda','175klambda','45.0deg'],
      niter=10000,threshold='10.0mJy',usescratch=True)


image=['HLTau_HCN1_0.tap200klR0.5_0.25kms',
       'HLTau_HCOp1_0.tap200klR0.5_0.25kms_ms']

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

####################################################
# Example Moment maps

# blue start 1.25 km/s ch 29
# red ends at 22.5 km/s 114
# red drops significantly "close-in" at 12.75 km/s ch 75
# The rms in line channels ranges from about 2.5 to 2.8 mJy/beam

name='HLTau_HCOp1_0.tap200klR0.5_0.25kms_ms.image'
rms=0.0025

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

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

os.system('rm -rf '+name+'.mom.*')
immoments(imagename=name,outfile=name+'.mom',
          chans='29~114',moments=[1,2],
          includepix=[rms*5.5,1000])

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

exportfits(name+'.mom0mask',name+'.mom0mask.fits')
exportfits(name+'.mom.weighted_coord',name+'.mom.weighted_coord.fits')


name='HLTau_HCN1_0.tap200klR0.5_0.25kms.image'
rms=0.0025

os.system('rm -rf '+name+'.mom0')
immoments(imagename=name,outfile=name+'.mom0',
          chans='46~53',moments=[0])
exportfits(name+'.mom0',name+'.mom0.fits')


#######################################################################
# Make a higher resolution cube that resolves out the outflow emission.


name='HLTau_HCOp_taper0.25_R0_0.25kms'
os.system('rm -rf '+name+'.*')
clean('HLTau_B3HCN_1_7.fixed.calHCOp.contsub',
      imagename=name,
      imagermode='csclean',weighting='briggs',robust=0.0,
      imsize=500, cell='0.04arcsec',
      mode='velocity',start='-2.0km/s',width='0.25km/s',nchan=72,
      outframe='LSRK',restfreq='89.18853GHz',
      interactive=True,
      uvtaper=True,outertaper='0.25arcsec',
      mask='combined.mask',
      niter=100000,threshold='3.8mJy',usescratch=True)

image=['HLTau_HCOp_taper0.25_R0_0.25kms']
for i in image:
    exportfits(i+'.image',i+'.image.fits')
    exportfits(i+'.mask',i+'.mask.fits')
    exportfits(i+'.flux',i+'.flux.fits')


#######################################################################
# Example spectral smoothing and moment maps   

#If desired spectral smoothing can be applied in CASA 4.3 but a
# bug means you have to reinsert the beam if you want to make
# an measurements on it.
specsmooth(imagename='HLTau_HCOp_taper0.25_R0_0.25kms.image',
           outfile='HLTau_HCOp_taper0.25_R0_0.25kms.image.specsmooth',
           function='boxcar',width=2)

ia.open('HLTau_HCOp_taper0.25_R0_0.25kms.image.specsmooth')
ia.setrestoringbeam(major='0.2765arcsec',minor='0.220632arcsec',pa='40.53deg')
ia.close()

# Example moment maps
name='HLTau_HCOp_taper0.25_R0_0.25kms.image.specsmooth'
rms=0.0019
chans='5~24'

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

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

# see CO script for example of primary beam correcting moment maps.