# 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')