# Antennae band 6 imaging script for CASA 3.3

#Zeroing pointing table

tb.open('Antennae_North.cal.ms/POINTING',nomodify=False)
a = tb.rownumbers()
tb.removerows(a)
tb.close()

tb.open('Antennae_South.cal.ms/POINTING',nomodify=False)
a = tb.rownumbers()
tb.removerows(a)
tb.close()

# Check the line
plotms(vis='Antennae_North.cal.ms',xaxis='channel',yaxis='amp',
      avgtime='1e8',avgscan=T,plotfile='Antennae_North-AMPvsCH.png')

dummy_string = raw_input("Hit <Enter> to see next data set")

plotms(vis='Antennae_South.cal.ms',xaxis='channel',yaxis='amp',
      avgtime='1e8',avgscan=T,plotfile='Antennae_South-AMPvsCH.png')

dummy_string = raw_input("Hit <Enter> to see next data set")

#Northern mosaic. Center of the mosaic: 7           
#Southern mosaic. Center of the mosaic: 9         

# Continuum. North
os.system('rm -rf Antennae_North.Cont.Dirty*')
clean(vis='Antennae_North.cal.ms',imagename='Antennae_North.Cont.Dirty',
     field='',phasecenter='7',
     mode='mfs',restfreq='230.538GHz',
     spw='0:1~50;120~164',
     imagermode='mosaic',
     imsize=300,cell='0.3arcsec',
     interactive=T,niter=0)

os.system('rm -rf Antennae_North.Cont.Clean*')
clean(vis='Antennae_North.cal.ms',imagename='Antennae_North.Cont.Clean',
     field='',phasecenter='7',
     mode='mfs',restfreq='230.538GHz',
     spw='0:1~50;120~164',
     imagermode='mosaic',
     imsize=300,cell='0.3arcsec',
     interactive=T, 
     niter=100,threshold='0.1mJy')

# Continuum. South
os.system('rm -rf Antennae_South.Cont.Clean*')
clean(vis='Antennae_South.cal.ms',imagename='Antennae_South.Cont.Clean',
     field='',phasecenter='9',
     mode='mfs',restfreq='230.538GHz',
     spw='0:1~40;130~164',
     imagermode='mosaic',
     imsize=400,cell='0.3arcsec',
     interactive=T,niter=0)

os.system('rm -rf Antennae_South.Cont.Clean*')
clean(vis='Antennae_South.cal.ms',imagename='Antennae_South.Cont.Clean',
     field='',phasecenter='9',
     mode='mfs',restfreq='230.538GHz',
     spw='0:1~40;130~164',
     imagermode='mosaic',
     imsize=400,cell='0.3arcsec',
     interactive=T,
     niter=100,threshold='0.1mJy')

#Some statistics on the continuum maps.
print 'Statistics in the Northern mosaic:'
imstat('Antennae_North.Cont.Clean.image')

print 'Statistics in the Southern mosaic:'
imstat('Antennae_South.Cont.Clean.image')

print 'Plotting Northern and Southern mosaics:'

imview (
  raster={'file': 'Antennae_North.Cont.Clean.image',
  'colorwedge':T,'range':[-0.0004,0.004]},
  zoom=1, out='Antennae_North.Cont.Clean.image.png')
imview (
  raster={'file': 'Antennae_South.Cont.Clean.image',
  'colorwedge':T,'range':[-0.0004,0.004]},
  zoom=1, out='Antennae_South.Cont.Clean.image.png')


#Continuum subtraction

uvcontsub2(vis='Antennae_North.cal.ms',fitspw='0:1~50;120~164',fitorder = 1)
uvcontsub2(vis='Antennae_South.cal.ms',fitspw='0:1~40;130~164',fitorder = 1)

plotms(vis='Antennae_North.cal.ms.contsub',xaxis='velocity',yaxis='amp',
       avgtime='1e8',avgscan=T,transform=T,freqframe='LSRK',
       restfreq='230.538GHz',plotfile='North_CO2_1_vel.png')
dummy_string = raw_input("Hit <Enter>")

plotms(vis='Antennae_South.cal.ms.contsub',xaxis='velocity',yaxis='amp',
      avgtime='1e8',avgscan=T,transform=T,freqframe='LSRK',
      restfreq='230.538GHz',plotfile='South_CO2_1_vel.png')
dummy_string = raw_input("Hit <Enter>")

os.system('rm -rf Antennae_North.CO2_1Line*')
clean(vis='Antennae_North.cal.ms.contsub',
     imagename='Antennae_North.CO2_1Line.Clean',
     spw='0',field='',phasecenter='7',
     mode='velocity',outframe='LSRK',restfreq='230.538GHz',
     nchan=35,start='1200km/s',width='20km/s',
     imagermode='mosaic',
     imsize=300,cell='0.3arcsec',minpb=0.2,
     interactive=T,
     weighting='briggs',robust=0.5,
     niter=10000, threshold='1.0mJy')

viewer('Antennae_North.CO2_1Line.Clean.image')
dummy_string = raw_input("Hit <Enter> to clean the CO(2-1) emission of the Southern mosaic")

os.system('rm -rf Antennae_South.CO2_1Line*')
clean(vis='Antennae_South.cal.ms.contsub',
     imagename='Antennae_South.CO2_1Line.Clean',
     spw='0',field='',phasecenter='9',
     mode='velocity',outframe='LSRK',restfreq='230.538GHz',
     nchan=35,start='1200km/s',width='20km/s',
     imagermode='mosaic',
     imsize=400,cell='0.3arcsec',minpb=0.2,
     interactive=T,
     weighting='briggs',robust=0.5,
     niter=10000, threshold='2.0mJy')

viewer('Antennae_South.CO2_1Line.Clean.image')
dummy_string = raw_input("Hit <Enter>")

#Immoments

os.system('rm -rf Antennae_North.CO2_1Line.Clean.image.mom* ')
os.system('rm -rf Antennae_South.CO2_1Line.Clean.image.mom* ')

immoments('Antennae_North.CO2_1Line.Clean.image', 
          moments=[0],chans='15~32',
          outfile='Antennae_North.CO2_1Line.Clean.image.mom0')

immoments('Antennae_South.CO2_1Line.Clean.image',
          moments=[0], chans='6~32', 
          outfile='Antennae_South.CO2_1Line.Clean.image.mom0')

immoments('Antennae_North.CO2_1Line.Clean.image', 
          moments=[1,2], chans='15~32', 
          mask='Antennae_North.CO2_1Line.Clean.image>0.001*5',
          outfile='Antennae_North.CO2_1Line.Clean.image.mom')

immoments('Antennae_South.CO2_1Line.Clean.image', 
          moments=[1,2], chans='6~32',
          mask=' Antennae_South.CO2_1Line.Clean.image>0.002*5',
          outfile='Antennae_South.CO2_1Line.Clean.image.mom')

os.system('rm -f Antennae_North.CO2_1Line.Clean.image.mom0.png')
imview (raster={'file': 'Antennae_North.CO2_1Line.Clean.image.mom0',
                 'colorwedge':T,'scaling': -0.5},
         zoom=1,out='Antennae_North.CO2_1Line.Clean.image.mom0.png')
dummy_string = raw_input("Hit <Enter>")

os.system('rm -f Antennae_North.CO2_1Line.Clean.image.mom.weighted_coord.png')
imview (raster={'file': 'Antennae_North.CO2_1Line.Clean.image.mom.weighted_coord',
                 'colorwedge':T},
         zoom=1,out='Antennae_North.CO2_1Line.Clean.image.mom.weighted_coord.png')


os.system('rm -f Antennae_South.CO2_1Line.Clean.image.mom0.png')
imview (raster={'file': 'Antennae_South.CO2_1Line.Clean.image.mom0',
                 'colorwedge':T,},contour={'file': 'Antennae_South.CO2_1Line.Clean.image.mom0'}, zoom=1,out='Antennae_South.CO2_1Line.Clean.image.mom0.png')

os.system('rm -f Antennae_South.CO2_1Line.Clean.image.mom.weighted_coord.png')
imview (raster={'file': 'Antennae_South.CO2_1Line.Clean.image.mom.weighted_coord',
                 'colorwedge':T},zoom=1,out='Antennae_South.CO2_1Line.Clean.image.mom.weighted_coord.png')


#Exporting fits
exportfits(imagename='Antennae_North.Cont.Clean.image',
     fitsimage='Antennae_North.Cont.Clean.image.fits')
exportfits(imagename='Antennae_North.CO2_1Line.Clean.image',
     fitsimage='Antennae_North.CO2_1Line.Clean.image.fits')
exportfits(imagename='Antennae_North.CO2_1Line.Clean.image.mom0',
     fitsimage='Antennae_North.CO2_1Line.Clean.image.mom0.fits')
exportfits(imagename='Antennae_North.CO2_1Line.Clean.image.mom.weighted_coord',
     fitsimage='Antennae_North.CO2_1Line.Clean.image.mom.weighted_coord.fits')
exportfits(imagename='Antennae_North.CO2_1Line.Clean.image.mom.weighted_dispersion_coord',
     fitsimage='Antennae_North.CO2_1Line.Clean.image.mom.weighted_dispersion_coord.fits')

exportfits(imagename='Antennae_South.Cont.Clean.image',
     fitsimage='Antennae_South.Cont.Clean.image.fits')
exportfits(imagename='Antennae_South.CO2_1Line.Clean.image',
     fitsimage='Antennae_South.CO2_1Line.Clean.image.fits')
exportfits(imagename='Antennae_South.CO2_1Line.Clean.image.mom0',
     fitsimage='Antennae_South.CO2_1Line.Clean.image.mom0.fits')
exportfits(imagename='Antennae_South.CO2_1Line.Clean.image.mom.weighted_coord',
     fitsimage='Antennae_South.CO2_1Line.Clean.image.mom.weighted_coord.fits')
exportfits(imagename='Antennae_South.CO2_1Line.Clean.image.mom.weighted_dispersion_coord',
     fitsimage='Antennae_South.CO2_1Line.Clean.image.mom.weighted_dispersion_coord.fits')