# The following assumes standard calibration of the amplitude calibrator and
# science target via the script generator.
# Examples of each are provided in:
# Amplitude Calibrator: uid___A002_X86fcfa_X3ae.ms.scriptForSDCalibration.py
# Science target: uid___A002_X86fcfa_X664.ms.scriptForSDCalibration.py
#
# Below is the procedure to achieve a fully calibrated TP Science image,
# in Jansky/beam units, taking as input a set of calibrated Science and
# Amplitude Calibrator datasets (in units of Kelvin).
#
# 1. Image the amplitude calibrator and measure the value of Jy/K from the image
# (values should be obtained for each night data is taken and for each spw and antenna)
#
# 2. Convert the science target units from Kelvin to Jansky
# (treating each night's data and each spw separately)
#
# 3. Image the science target
# (combining data from all antennas)
#
# 4. Optionally subtract a residual background from the image
# (currently, this is likely necessary in most cases; in future, improvements
# in the sdbaseline task should remove the need for this step)
#
# 5. Apply the restoring beam to the science image

### scriptForImagingScienceTarget.py deals with STEPs 3, 4, and 5. ###
### Tested in CASA Version 4.3.0                                   ###

# Prior to running this script, scriptForImagingAmpCalAndDerivingJyPerK.py
# and scriptForJyPerKConversion.py must have been run.

### 3. Image the science target ###

# Image the science target for the chosen spectral lines
print 'Imaging the science target...'

# give the list of the science datasets (output of script generator calibration)
sciencedata = ['uid___A002_X85c183_X36f.ms.cal.jy',   # 2014-07-01
               'uid___A002_X85c183_X60b.ms.cal.jy',   # 2014-07-01
               'uid___A002_X8602fa_X2ab.ms.cal.jy',   # 2014-07-05
               'uid___A002_X8602fa_X577.ms.cal.jy',   # 2014-07-05
               'uid___A002_X864236_X2d4.ms.cal.jy',   # 2014-07-07
               'uid___A002_X864236_X693.ms.cal.jy',   # 2014-07-07
               'uid___A002_X86fcfa_Xd9.ms.cal.jy',    # 2014-07-17
               'uid___A002_X86fcfa_X664.ms.cal.jy',   # 2014-07-17
               'uid___A002_X86fcfa_X96c.ms.cal.jy']   # 2014-07-17

sciencespw = [23]

# Define some properties of the beam
fwhmfactor = 1.13
diameter = 12

# Determine from the pointing table the sampling used to observe the data
# (just use the first dataset in the list - they must all be the same)
xSampling, ySampling, maxsize = aU.getTPSampling(sciencedata[0], showplot=False)

for spw in sciencespw:
    # Determine the frequency of each science spw
    msmd.open(sciencedata[0])
    freq = msmd.meanfreq(spw)
    msmd.close()
    print "SPW %d: %.3f GHz" % (spw, freq*1e-9)

    # Determine cell size (using a Gaussian approximation)
    theorybeam = aU.primaryBeamArcsec(frequency=freq*1e-9,
                                      fwhmfactor=fwhmfactor,
                                      diameter=diameter)

    # Set up the imaging parameters
    # Note that the parameters must be the same as for the amplitude calibrator
    cell = theorybeam/9.0
    imsize = int(round(maxsize/cell)*2)
    # For M100 we want an image 500 arcsec by 500 arcsec

    sdimaging(infiles=sciencedata,
              field='M100',
              spw='%d' % spw,  #sciencespw
              nchan=70,
              mode='velocity',
              start='1400km/s',
              width='5km/s',
              veltype='radio',
              outframe='lsrk',
              restfreq='115.271204GHz',
              gridfunction='SF',
              convsupport=6,
              stokes='',
              phasecenter='J2000 12h22m54.9 +15d49m15',
              ephemsrcname='', 
              imsize=imsize,
              cell=str(cell)+'arcsec',
              overwrite=True,
              outfile='M100.CO.cube.image') # modify outfile if more than one spw being imaged

    # Put the updated units (Jy/beam) in the image header
    print 'Updating the brightness unit of the image ...'
    imhead(imagename='M100.CO.cube.image',
           mode='put',
           hdkey='bunit',
           hdvalue='Jy/beam')

    
    ### 4. Optionally subtract a residual background from the image ###

    print 'Subtracting baseline (manually) ...'
    os.system('rm -rf M100.CO.cube.bl.image M100.ignorethis.image')
    imcontsub(imagename='M100.CO.cube.image',
              linefile='M100.CO.cube.bl.image',
              contfile='M100.ignorethis.image',
              fitorder=1,
              chans='0~7;62~69')
    os.system('rm -rf M100.ignorethis.image')


    ### 5. Apply the restoring beam to the science image ###
    minor, major, fwhmsfBeam, sfBeam = aU.sfBeam(frequency=freq*1e-9,
            pixelsize=cell,
            convsupport=6,
            img=None, #to use Gaussian theorybeam
            stokes='both',
            xSamplingArcsec=xSampling,
            ySamplingArcsec=ySampling,
            fwhmfactor=fwhmfactor,
            diameter=diameter)

    ia.open('M100.CO.cube.bl.image')
    ia.setrestoringbeam(major=str(sfBeam)+'arcsec', minor=str(sfBeam)+'arcsec', pa='0deg')
    ia.done()

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