#!/usr/bin/env python
import sys

sys.path.append('/opt/share/casa/utils/science/analysis_scripts')
import analysisUtils as aU

asdmlist = ['uid___A002_X85c183_X36f',
            'uid___A002_X85c183_X60b',
            'uid___A002_X8602fa_X2ab',
            'uid___A002_X8602fa_X577',
            'uid___A002_X864236_X2d4',
            'uid___A002_X864236_X693',
            'uid___A002_X86fcfa_X664',
            'uid___A002_X86fcfa_X96c',
            'uid___A002_X86fcfa_Xd9']
vislist = map(lambda x: x + '.ms.bl', asdmlist)

fwhmfactor = 1.13
diameter = 12

refvis = asdmlist[0] + '.ms'

xSampling, ySampling, maxsize = aU.getTPSampling(refvis, showplot=False)

for spw in [3]: # this spectral window contains the 12CO(1-0) transition
    msmd.open(vislist[0])
    freq = msmd.meanfreq(spw)
    msmd.close()
    print "SPW %d: %.3f GHz" % (spw, freq*1e-9)
    
    theorybeam = aU.primaryBeamArcsec(frequency=freq*1e-9,
                                  fwhmfactor=fwhmfactor,
                                  diameter=diameter)
 
    cell = theorybeam/9.0
    imsize = int(round(maxsize/cell)*2)

    imagename = 'M100_TP_CO_cube.spw%s.image'%(spw)
    
    sdimaging(infiles=vislist,
              field='M100',
              spw='%s'%(spw),
              nchan=70,
              mode='velocity',
              start='1400km/s',
              width='5km/s',
              veltype='radio',
              outframe='lsrk',
              restfreq='115.271204GHz',
              gridfunction='SF',
              convsupport=6,
              stokes='I',
              phasecenter='J2000 12h22m54.9 +15d49m15',
              ephemsrcname='',
              imsize=imsize,
              cell='%sarcsec'%(cell),
              overwrite=True,
              outfile=imagename)

    # Put the updated units (Jy/beam) in the image header
    print 'Updating the brightness unit of the image ...'
    imhead(imagename=imagename,
           mode='put',
           hdkey='bunit',
           hdvalue='Jy/beam')
    
    print 'Subtracting baseline (manually) ...'
    blimagename = imagename + '.bl'
    os.system('rm -rf %s M100.ignorethis.image'%(blimagename))
    imcontsub(imagename=imagename,
              linefile=blimagename,
              contfile='M100.ignorethis.image',
              fitorder=1,
              chans='0~7;62~69')
    os.system('rm -rf M100.ignorethis.image')



    # Make an integrated intensity map
    os.system('rm -rf *M100.CO.cube.bl.image.mom0')
    integimagename = imagename + '.integ'
    os.system('rm -rf %s'%(integimagename))
    immoments(imagename=blimagename,
              moments=[0],
              axis='spectral',
              chans='8~61',
              outfile=integimagename)