# Script to image ALMA spectral line data for disk around T-Tauri star TW Hya

# Split out continuum data

os.system('rm -rf TWHydra_cont.ms*')
split(vis='TWHydra_corrected.ms',
      outputvis='TWHydra_cont.ms',
      spw='0~3:7~1273', width=30,
      datacolumn='data')

# Plot the continuum data amplitude with channel

plotms(vis='TWHydra_cont.ms', spw='0~3',
       xaxis='channel', yaxis='amp', 
       avgtime='1e8', avgscan=T, 
       coloraxis='spw', iteraxis='spw', 
       xselfscale=T)

# Flag the line channels in the continuum dataset

flagdata(vis='TWHydra_cont.ms', mode='manual', 
         spw='0:18~18, 2:23~24, 3:33~42')

# Run listobs on the full dataset

listobs('TWHydra_corrected.ms', listfile='TWHydra_corrected.ms.listobs')

# Run script to determine time on source (for determining sensitivity)

execfile('time_on_source.py')

# Initially clean the continuum data

os.system('rm -rf TWHydra_contall.*')
clean(vis='TWHydra_cont.ms',
      imagename='TWHydra_contall',
      mode='mfs', imagermode='csclean',
      imsize=100, cell=['0.3arcsec'], spw='',
      weighting='briggs', robust=0.5, 
      mask='', usescratch=False, interactive=T,
      threshold='0.6mJy', niter=10000)

# Split out the 12CO(3-2) data

os.system('rm -rf TWHydra_CO3_2.ms*')
split(vis='TWHydra_corrected.ms',
      outputvis='TWHydra_CO3_2.ms', 
      datacolumn='data', spw='2')

# Split out the HCO+ data

os.system('rm -rf TWHydra_HCOplus.ms*')
split(vis='TWHydra_corrected.ms',
      outputvis='TWHydra_HCOplus.ms', 
      datacolumn='data', spw='0')

# Find line free-channels for both datasets

os.system('rm -rf CO3_2_channel.png')
plotms(vis='TWHydra_CO3_2.ms', 
       spw='0',xaxis='channel', yaxis='amp', 
       avgtime='1e8', avgscan=T, coloraxis='spw', 
       plotfile='CO3_2_channel.png')

uvcontsub(vis='TWHydra_CO3_2.ms', 
          fitorder=1, fitspw='0:6~630,0:800~1265')

os.system('rm -rf HCOplus_channel.png')
plotms(vis='TWHydra_HCOplus.ms', 
       spw='0',xaxis='channel', yaxis='amp', 
       avgtime='1e8', avgscan=T, coloraxis='spw', 
       plotfile='HCOplus_channel.png')

uvcontsub(vis='TWHydra_HCOplus.ms', 
          fitorder=1, fitspw='0:6~630,0:800~1265')

# Plot the continuum subtracted data as a function of velocity

os.system('rm -rf CO3_2_vel.png')
plotms(vis='TWHydra_CO3_2.ms.contsub',
       xaxis='velocity', yaxis='amp', avgtime='1e8',
       avgscan=T, transform=T, freqframe='LSRK',
       restfreq='345.79599GHz', plotrange=[-20,23,0,0], 
       plotfile='CO3_2_vel.png')

os.system('rm -rf HCOplus_vel.png')
plotms(vis='TWHydra_HCOplus.ms.contsub',
       xaxis='velocity', yaxis='amp', avgtime='1e8',
       avgscan=T, transform=T, freqframe='LSRK',
       restfreq='356.7342GHz', plotrange=[-20,23,0,0], 
       plotfile='HCOplus_vel.png')

# 12CO(3-2) imaging

os.system('rm -rf TWHydra_CO3_2line.*')
clean(vis='TWHydra_CO3_2.ms.contsub',  
      imagename='TWHydra_CO3_2line', imagermode='csclean',
      spw='', imsize=100, cell=['0.3arcsec'], 
      mode='velocity', start='-4km/s', width='0.32km/s',
      nchan=40, restfreq='345.79599GHz', outframe='LSRK',
      weighting='briggs', robust=0.5, 
      mask='', usescratch=False, interactive=T, 
      threshold='33mJy', niter=100000)

# HCO+ imaging

os.system('rm -rf TWHydra_HCOplusline.*')
clean(vis='TWHydra_HCOplus.ms.contsub',  
      imagename='TWHydra_HCOplusline', imagermode='csclean',
      spw='', imsize=100, cell=['0.3arcsec'], 
      mode='velocity', start='-4km/s', width='0.32km/s',
      nchan=40, restfreq='356.7342GHz', outframe='LSRK',
      weighting='briggs', robust=0.5, 
      mask='', usescratch=False, interactive=T, 
      threshold='33mJy', niter=100000)

# Determine synthesised beam for images using imhead 
# (if you didn't remember to note it down from logger during clean)

imhead("TWHydra_CO3_2line.image")
imhead("TWHydra_HCOplusline.image")

# Make spectrum using spectral profile tool in viewer, e.g. for 12CO
# (use elliptical region)

viewer("TWHydra_CO3_2line.image")

# Also fit a gaussian to the line (interactively)

# Save region using View -> Regions -> File tab
# (Regions panel may already be open)

# Fit gaussian to spectrum in region, e.g.

specfit(imagename='TWHydra_CO3_2line.image',
        region='spectrum_region.crtf', poly=-1,
        logresults=True)

# Open image in viewer to estimate spectral extent
# of the emission

viewer("TWHydra_CO3_2line.image")

# (need to use regions to determine the rms first)

results = imstat("TWHydra_CO3_2line.image", chans="7")
print results
print "  s.d. ", results['sigma']
print "  RMS ", results['rms']

os.system('rm -rf TWHydra_CO3_2line.image.mom0')
immoments(imagename='TWHydra_CO3_2line.image',  
          outfile='TWHydra_CO3_2line.image.mom0', 
          moments=[0], chans='13~32')

# View moment zero maps in viewer

imview( raster= {'file':'TWHydra_CO3_2line.image.mom0',
                  'range':[-1.,10.]}, 
         contour={'file':'TWHydra_contall.image', 
                  'base':0, 'unit':0.0025, 
                  'levels':[3,100]} )

# Make first moment map for 12CO

os.system('rm -rf TWHydra_CO3_2line.image.mom1')
immoments(imagename='TWHydra_CO3_2line.image',moments=[1],
          outfile='TWHydra_CO3_2line.image.mom1',
          chans='13~32',includepix=[0.5,100])

# Make second moment map for 12CO

os.system('rm -rf TWHydra_CO3_2line.image.mom2')
immoments(imagename='TWHydra_CO3_2line.image',moments=[2],
          outfile='TWHydra_CO3_2line.image.mom2',
          chans='13~32',includepix=[0.5,100])

# View the moment maps using viewer

imview( raster=[ {'file':'TWHydra_CO3_2line.image.mom0'},
                 {'file':'TWHydra_CO3_2line.image.mom1'},
                 {'file':'TWHydra_CO3_2line.image.mom2'} ], 
        contour={'file':'TWHydra_contall.image', 
                 'base':0, 'unit':0.0025, 
                 'levels':[3,100]} )

# Export an image using exportfits

os.system('rm -rf TWHydra_CO3_2line.image.fits')
exportfits(imagename='TWHydra_CO3_2line.image',
           fitsimage='TWHydra_CO3_2line.image.fits')

# Correct for the primary beam response

impbcor(imagename='TWHydra_contall.image',
       pbimage='TWHydra_contall.flux',
       mode='divide', 
       outfile='TWHydra_contall.pbcor')

# Fit the continuum emission with a 2D gaussian

imfit(imagename="TWHydra_contall.image",
      box="40,40,60,60", logfile = "contin_fit.log", 
      residual="TWHydra_contall.fitresid")

# Make pv diagram using task impv

os.system('rm -rf TWHydra_CO3_2line.image.pv')
impv(imagename='TWHydra_CO3_2line.image', 
               mode='length', center=[50,49], 
               length='10arcsec', width='8arcsec',
               pa='-35deg', chans='12~30',
               outfile='TWHydra_CO3_2line.image.pv')

# Reproject an image to Galactic coordinates

imregrid(imagename='TWHydra_CO3_2line.image',
         template='GALACTIC',
         output='TWHydra_CO3_2line.Galactic')