###############################################
#                                             #
# Regression/Benchmarking Script for L02D     #
#             26 JAN               (U Cam)    #
###############################################
import time
import os
# Data must be pre-filled
# l.ms, l2.ms

pathname=os.environ.get('AIPSPATH')
datapath=pathname+'/data/regression/ATST1/L02D/l*.ms'

restore()

startTime=time.time()
startProc=time.clock()

print '--Copy MS--'
# Fill - really copy over pristine MS
os.system('rm -rf l.* l2.* all.* l1.* l02d.*')
copystring='cp -r '+datapath+' .'
os.system(copystring)
clearcal('l.ms')
clearcal('l2.ms')
copytime=time.time()

print '--Flag data--'
default('flagdata')
flagdata('l.ms',timerange=['2002/01/27/05:45:47.0', '2002/01/27/07:00:00.0'],fieldid=-1)
flagtime=time.time()

#
# l1 calibration 
#
print '--Calibration phase/bandpass (3mm)--'
#setjy
default('setjy')
setjy('l.ms',fieldid=0,spwid=6,fluxdensity=[10.87,0.,0.,0.])
setjy('l.ms',fieldid=0,spwid=14,fluxdensity=[10.87,0.,0.,0.])
setjy('l.ms',fieldid=0,spwid=18,fluxdensity=[10.87,0.,0.,0.])
setjy('l.ms',fieldid=0,spwid=7,fluxdensity=[10.87,0.,0.,0.])
setjy('l.ms',fieldid=0,spwid=15,fluxdensity=[10.87,0.,0.,0.])
setjy('l.ms',fieldid=0,spwid=19,fluxdensity=[10.87,0.,0.,0.])
setjy('l.ms',fieldid=0,spwid=3,fluxdensity=[10.87,0.,0.,0.])
setjy('l.ms',fieldid=0,spwid=11,fluxdensity=[10.87,0.,0.,0.])
setjy('l.ms',fieldid=0,spwid=26,fluxdensity=[6.0,0,0,0])
setjy('l.ms',fieldid=0,spwid=30,fluxdensity=[6.0,0,0,0])
setjy('l.ms',fieldid=0,spwid=27,fluxdensity=[6.0,0,0,0])
setjy('l.ms',fieldid=0,spwid=31,fluxdensity=[6.0,0,0,0])

# 3mm USB
#preliminary time-dependent phase solutions to improve coherent average for bandpass solution
default('gaincal')
gaincal('l.ms','l.3mmUSB.gcal0','channel',nchan=110,start=9,step=1,msselect='FIELD_ID==0 && DATA_DESC_ID IN [7,15,19]',gaintype='GSPLINE',calmode='p',splinetime=10000.,refant='1',phasewrap=260,gaincurve=False,opacity=False,preavg=0)

print '--bandpass (3mm)--'
#derive bandpass calibration for 3mm (C34S J=2-1) using USB continuum phase solns
default('bandpass')
bandpass('l.ms','l.3mmC34S.bpoly','none',msselect='FIELD_ID==0 && DATA_DESC_ID==3',bandtype='BPOLY',degamp=2,degphase=2,bpnorm=False,maskcenter=2,maskedge=0,refant='1',gaintable='l.3mmUSB.gcal0',gaincurve=False,opacity=False,solint=0.)

# Band pass for CH3OH,3mm continuum USB/LSB
default('bandpass')
bandpass('l.ms','l.3mmch3oh.bpoly','none',msselect='FIELD_ID==0 && DATA_DESC_ID==11',bandtype='BPOLY',degamp=2,degphase=4,bpnorm=False,maskcenter=2,maskedge=0,refant='1',gaintable='l.3mmUSB.gcal0',gaincurve=False,opacity=False)
default('bandpass')
bandpass('l.ms','l.3mmcont.bpoly','none',msselect='FIELD_ID==0 && DATA_DESC_ID IN [7,15,19]',bandtype='BPOLY',degamp=10,degphase=25,bpnorm=False,maskcenter=2,maskedge=0,refant='1',gaintable='l.3mmUSB.gcal0',gaincurve=False,opacity=False)
default('bandpass')
bandpass('l.ms','l.3mmcont.bpoly','none',msselect='FIELD_ID==0 && DATA_DESC_ID IN [6,14,18]',bandtype='BPOLY',degamp=10,degphase=25,bpnorm=False,maskcenter=2,maskedge=0,refant='1',gaintable='l.3mmUSB.gcal0',gaincurve=False,opacity=False,append=True)

print '--gaincal phase (3mm)--'
#derive new and better phase solutions for 3mm LSB
default('gaincal')
gaincal('l.ms','l.3mmcont.gcal','none',msselect='FIELD_ID IN [0,1,2] && DATA_DESC_ID IN [6,14,18,7,15,19]',gaintype='GSPLINE',calmode='p',splinetime=3000.,refant='1',phasewrap=260,gaincurve=False,opacity=False,bptable='l.3mmcont.bpoly',preavg=0.)

#Apply all solutions derived so far, determine calibrator's flux densities by solving for T 
#and use fluxscale
default('gaincal')
gaincal('l.ms','l.3mmcont.temp','none',msselect='FIELD_ID IN [0,1,2] && DATA_DESC_ID IN [7,15,19,6,14,18]',solint=600.,refant='1',gaintype='T',opacity=False,gaincurve=False,gaintable='l.3mmcont.gcal',bptable='l.3mmcont.bpoly')

#fluxscale
default('fluxscale')
fluxscale('l.ms','l.3mmcont.temp','l.3mmcont.flux',reference=['3C273   '],transfer=['MWC349  ','2013+370'])
calphase3mmtime=time.time()
#MWC349 1.14
#2013+370 3.25

#use old values
print '--Set fluxscale (setjy)--'
default('setjy')
setjy('l.ms',fieldid=1,spwid=6,fluxdensity=[1.01,0.,0.,0.])
setjy('l.ms',fieldid=1,spwid=7,fluxdensity=[1.01,0.,0.,0.])
setjy('l.ms',fieldid=1,spwid=14,fluxdensity=[1.01,0.,0.,0.])
setjy('l.ms',fieldid=1,spwid=15,fluxdensity=[1.01,0.,0.,0.])
setjy('l.ms',fieldid=1,spwid=18,fluxdensity=[1.01,0.,0.,0.])
setjy('l.ms',fieldid=1,spwid=19,fluxdensity=[1.01,0.,0.,0.])
setjy('l.ms',fieldid=2,spwid=6,fluxdensity=[2.79,0.,0.,0.])
setjy('l.ms',fieldid=2,spwid=7,fluxdensity=[2.79,0.,0.,0.])
setjy('l.ms',fieldid=2,spwid=14,fluxdensity=[2.79,0.,0.,0.])
setjy('l.ms',fieldid=2,spwid=15,fluxdensity=[2.79,0.,0.,0.])
setjy('l.ms',fieldid=2,spwid=18,fluxdensity=[2.79,0.,0.,0.])
setjy('l.ms',fieldid=2,spwid=19,fluxdensity=[2.79,0.,0.,0.])
setjy3mmtime=time.time()

## Amplitude calibration of 3mm continuum:
##
##  phase solutions will be pre-applied as well as carried forward 
##   to the output solution table.
print '--gaincal amp (3mm)--'
default('gaincal')
gaincal('l.ms','l.3mmcont.amp.gcal','none',msselect='FIELD_ID IN [0,1,2] && DATA_DESC_ID IN [7,15,19,6,14,18]',gaintype='GSPLINE',calmode='a',splinetime=3000.,refant='1',phasewrap=260,gaincurve=False,opacity=False,bptable='l.3mmcont.bpoly',preavg=0.,gaintable='l.3mmcont.gcal')
calamp3mmtime=time.time()

## Correct the target source C34S (J=2-1)
##
##  note that only the 60 central channels will be calibrated
##   since the BPOLY solution is only defined for these
default('correct')
correct('l.ms',msselect='FIELD_ID==3 && DATA_DESC_ID==3',gaincurve=False,opacity=False,gaintable='l.3mmcont.gcal',bptable='l.3mmC34S.bpoly',pointtable='l.3mmcont.amp.gcal')
# Correct Target other 3mm line data (CH3OH)
correct('l.ms',msselect='FIELD_ID==3 && DATA_DESC_ID==11',gaincurve=False,opacity=False,gaintable='l.3mmcont.gcal',bptable='l.3mmch3oh.bpoly',pointtable='l.3mmcont.amp.gcal')
# Correct Target source 3mm continuum 
correct('l.ms',msselect='FIELD_ID==3 && DATA_DESC_ID IN [7,15,19,6,14,18]',gaincurve=False,opacity=False,gaintable='l.3mmcont.gcal',bptable='l.3mmcont.bpoly',pointtable='l.3mmcont.amp.gcal')
correct3mmtime=time.time()

# Split calibrated data
#print '--split calibrater--'
#default('split')
#split('l.ms','l.c34s.split.ms',fieldid=3,spwid=3,nchan=512,start=0,step=1,datacolumn='corrected')
#default('split')
#split('l.ms','l.ch3oh.split.ms',fieldid=3,spwid=11,nchan=512,start=0,step=1,datacolumn='corrected')
#default('split')
#split('l.ms','l.3mmcont.split.ms',fieldid=3,spwid=[6,7,14,15,18,19],nchan=1,start=10,step=100,datacolumn='corrected')
#splitsrctime=time.time()

# 1mm Calibration

###########################################################
## Get first cut phase solutions to improve S/N for BPass determination:
## 
default('gaincal')
gaincal('l.ms','l.1mmUSB.gcal0','channel',nchan=110,start=9,step=1,msselect='FIELD_ID==0 && DATA_DESC_ID IN [27,31]',gaintype='GSPLINE',calmode='p',splinetime=7000.,refant='1',phasewrap=260,preavg=0.,gaincurve=False,opacity=False)

## Derive bandpass calibration for 1mm 
##
bandpass('l.ms','l.1mmcont.bpoly','none',msselect='FIELD_ID==0 && DATA_DESC_ID IN [27,31]',bandtype='BPOLY',degamp=10,degphase=20,visnorm=False,bpnorm=False,maskcenter=2,maskedge=0,refant='1',gaintable='l.1mmUSB.gcal0',gaincurve=False,opacity=False)
bandpass('l.ms','l.1mmcont.bpoly','none',msselect='FIELD_ID==0 && DATA_DESC_ID IN [26,30]',bandtype='BPOLY',degamp=10,degphase=20,visnorm=False,bpnorm=False,maskcenter=2,maskedge=0,refant='1',gaintable='l.1mmUSB.gcal0',gaincurve=False,opacity=False,append=True)
bandpass('l.ms','l.1mmUSB.bpoly','none',msselect='FIELD_ID==0 && DATA_DESC_ID IN [23]',bandtype='BPOLY',degamp=10,degphase=20,visnorm=False,bpnorm=False,maskcenter=2,maskedge=0,refant='1',gaintable='l.1mmUSB.gcal0',gaincurve=False,opacity=False)

## Determine phase solutions for 1mm LSB & USB 
##
print '--Phase solutions--'
default('gaincal')
gaincal('l.ms','l.1mmcont.gcal','none',msselect='FIELD_ID IN [0,1,2] && DATA_DESC_ID IN [26,27,30,31]',gaintype='GSPLINE',calmode='p',splinetime=3000.,refant='1',phasewrap=260,gaincurve=False,opacity=False,bptable='l.1mmcont.bpoly',preavg=0.)

##  Apply all solutions derived so far, determine
##  calibrators' flux densities using a solve for T and
##  fluxscale
print '--Solve for solutions, fluxscale--'
default('gaincal')
gaincal('l.ms','l.1mmcont.temp','none',msselect='FIELD_ID IN [0,1,2] && DATA_DESC_ID IN [26,27,30,31]',solint=1200.,refant='1',gaintype='T',opacity=False,gaincurve=False,gaintable='l.1mmcont.gcal',bptable='l.1mmcont.bpoly')
#
default('fluxscale')
fluxscale('l.ms','l.1mmcont.temp','l.1mmcont.flux',reference=['3C273   '],transfer=['MWC349  ','2013+370'])
#MWC349   2.7	1.74
#2013+370 3.8	2.44
calphase1mmtime=time.time()

## Record flux values from logger window.  Manually insert
## fluxes with imgr.setjy:
print '--Setjy 1mm --'
setjy('l.ms',fieldid=1,spwid=27,fluxdensity=[1.7,0.,0.,0.])
setjy('l.ms',fieldid=1,spwid=31,fluxdensity=[1.7,0.,0.,0.])
setjy('l.ms',fieldid=2,spwid=27,fluxdensity=[1.8,0.,0.,0.])
setjy('l.ms',fieldid=2,spwid=31,fluxdensity=[1.8,0.,0.,0.])
setjy('l.ms',fieldid=1,spwid=27,fluxdensity=[1.7,0.,0.,0.])
setjy('l.ms',fieldid=1,spwid=31,fluxdensity=[1.7,0.,0.,0.])
setjy('l.ms',fieldid=2,spwid=27,fluxdensity=[1.8,0.,0.,0.])
setjy('l.ms',fieldid=2,spwid=31,fluxdensity=[1.8,0.,0.,0.])
setjy1mmtime=time.time()


## Amplitude calibration of 1mm LSB:
default('gaincal')
gaincal('l.ms','l.1mmcont.amp.gcal','none',msselect='FIELD_ID IN [0,1,2] && DATA_DESC_ID IN [26,27,30,31]',gaintype='GSPLINE',calmode='a',splinetime=3000.,refant='1',phasewrap=260,gaincurve=False,opacity=False,bptable='l.1mmcont.bpoly',preavg=0.,gaintable='l.1mmcont.gcal')
calamp1mmtime=time.time()

## Correct the target source 1mm LSB & USB continuum data:
##
##  Note that edge channels in won't be calibrated because
##   BPOLY solution is not defined for them
default('correct')
correct('l.ms',msselect='FIELD_ID==3 && DATA_DESC_ID IN [26,27,30,31]',gaincurve=False,opacity=False,gaintable='l.1mmcont.gcal',bptable='l.1mmcont.bpoly',pointtable='l.1mmcont.amp.gcal')
default('correct')
correct('l.ms',msselect='FIELD_ID==3 && DATA_DESC_ID IN [23]',gaincurve=False,opacity=False,gaintable='l.1mmcont.gcal',bptable='l.1mmcont.bpoly',pointtable='l.1mmcont.amp.gcal')
correct1mmtime=time.time()

## Split out calibrated target source 1 mm continuum data:
default('split')
split('l.ms','l.1mm.split.ms',fieldid=3,spwid=[26,27,30,31],nchan=110,start=8,step=1,datacolumn='corrected')
default('split')
split('l.ms','l.1mmc34s5_4.split.ms',fieldid=3,spwid=[23],nchan=512,start=0,step=1,datacolumn='corrected')
splitsrc1mmtime=time.time()

print '--Image 1mm continuum--'
## Get a first image the target source in 1 mm continuum emission:
default('clean')
clean('l.1mm.split.ms',imagename='l.1mm',alg='clark',niter=100,gain=0.1,nchan=1,start=0,width=1,spwid=[0,1,2,3],fieldid=0,stokes='I',weighting='natural',cell=[0.2,0.2],imsize=[128,128],mode='mfs')
image1mmtime=time.time()

print '--Resplit data--'
default('split')
split('l.ms','l1.c34s.split.ms',fieldid=3,spwid=3,nchan=512,start=0,step=1,datacolumn='corrected')
default('split')
split('l.ms','l1.ch3oh.split.ms',fieldid=3,spwid=[11],nchan=512,start=0,step=1,datacolumn='corrected')
default('split')
split('l.ms','l1.3mmcont.split.ms',fieldid=3,spwid=[6,7,14,15,18,19],nchan=1,start=9,step=100,datacolumn='corrected')
default('split')
split('l.ms','l1.1mmcont.split.ms',fieldid=3,spwid=[26,27,30,31],nchan=1,start=9,step=110,datacolumn='corrected')
default('split')
split('l.ms','l1.1mmc34s5_4.split.ms',fieldid=3,spwid=[23],nchan=512,start=0,step=1,datacolumn='corrected')
splitsrc1mmtime=time.time()



###############################################
#                                             #
# Regression/Benchmarking Script for L02D     #
#            29 JAN                (U Cam)    #
###############################################

print '--Flag data--'
default('flagdata')
flagdata('l2.ms',timerange=['2002/01/29/03:37:50.0', '2002/01/29/04:02:48.0'],fieldid=-1)
flagtime=time.time()

#
# l2 calibration 
#
print '--Calibration phase/bandpass (3mm)--'
#setjy
default('setjy')
setjy('l2.ms',fieldid=0,spwid=6,fluxdensity=[10.87,0.,0.,0.])
setjy('l2.ms',fieldid=0,spwid=14,fluxdensity=[10.87,0.,0.,0.])
setjy('l2.ms',fieldid=0,spwid=18,fluxdensity=[10.87,0.,0.,0.])
setjy('l2.ms',fieldid=0,spwid=7,fluxdensity=[10.87,0.,0.,0.])
setjy('l2.ms',fieldid=0,spwid=15,fluxdensity=[10.87,0.,0.,0.])
setjy('l2.ms',fieldid=0,spwid=19,fluxdensity=[10.87,0.,0.,0.])
setjy('l2.ms',fieldid=0,spwid=3,fluxdensity=[10.87,0.,0.,0.])
setjy('l2.ms',fieldid=0,spwid=11,fluxdensity=[10.87,0.,0.,0.])
setjy('l2.ms',fieldid=0,spwid=26,fluxdensity=[6.0,0,0,0])
setjy('l2.ms',fieldid=0,spwid=30,fluxdensity=[6.0,0,0,0])
setjy('l2.ms',fieldid=0,spwid=27,fluxdensity=[6.0,0,0,0])
setjy('l2.ms',fieldid=0,spwid=31,fluxdensity=[6.0,0,0,0])

# 3mm USB
#preliminary time-dependent phase solutions to improve coherent average for bandpass solution
default('gaincal')
gaincal('l2.ms','l2.3mmUSB.gcal0','channel',nchan=110,start=9,step=1,msselect='FIELD_ID==0 && DATA_DESC_ID IN [7,15,19]',gaintype='GSPLINE',calmode='p',splinetime=10000.,refant='1',phasewrap=260,gaincurve=False,opacity=False,preavg=0)

print '--bandpass (3mm)--'
#derive bandpass calibration for 3mm (C34S J=2-1) using USB continuum phase solns
default('bandpass')
bandpass('l2.ms','l2.3mmC34S.bpoly','none',msselect='FIELD_ID==0 && DATA_DESC_ID==3',bandtype='BPOLY',degamp=2,degphase=2,bpnorm=False,maskcenter=2,maskedge=0,refant='1',gaintable='l2.3mmUSB.gcal0',gaincurve=False,opacity=False)

# Band pass for CH3OH,3mm continuum USB/LSB
default('bandpass')
bandpass('l2.ms','l2.3mmch3oh.bpoly','none',msselect='FIELD_ID==0 && DATA_DESC_ID==11',bandtype='BPOLY',degamp=10,degphase=20,bpnorm=False,maskcenter=2,maskedge=0,refant='1',gaintable='l2.3mmUSB.gcal0',gaincurve=False,opacity=False)
default('bandpass')
bandpass('l2.ms','l2.3mmcont.bpoly','none',msselect='FIELD_ID==0 && DATA_DESC_ID IN [7,15,19]',bandtype='BPOLY',degamp=10,degphase=25,bpnorm=False,maskcenter=2,maskedge=0,refant='1',gaintable='l2.3mmUSB.gcal0',gaincurve=False,opacity=False)
default('bandpass')
bandpass('l2.ms','l2.3mmcont.bpoly','none',msselect='FIELD_ID==0 && DATA_DESC_ID IN [6,14,18]',bandtype='BPOLY',degamp=10,degphase=25,bpnorm=False,maskcenter=2,maskedge=0,refant='1',gaintable='l2.3mmUSB.gcal0',gaincurve=False,opacity=False,append=True)

print '--gaincal phase (3mm)--'
#derive new and better phase solutions for 3mm LSB
default('gaincal')
gaincal('l2.ms','l2.3mmcont.gcal','none',msselect='FIELD_ID IN [0,1,2] && DATA_DESC_ID IN [6,14,18,7,15,19]',gaintype='GSPLINE',calmode='p',splinetime=10000.,refant='1',phasewrap=260,gaincurve=False,opacity=False,bptable='l2.3mmcont.bpoly',preavg=0.)

#Apply all solutions derived so far, determine calibrator's flux densities by solving for T 
#and use fluxscale
default('gaincal')
gaincal('l2.ms','l2.3mmcont.temp','none',msselect='FIELD_ID IN [0,1,2] && DATA_DESC_ID IN [7,15,19,6,14,18]',solint=600.,refant='1',gaintype='T',opacity=False,gaincurve=False,gaintable='l2.3mmcont.gcal',bptable='l2.3mmcont.bpoly')

#fluxscale
default('fluxscale')
fluxscale('l2.ms','l2.3mmcont.temp','l2.3mmcont.flux',reference=['3C273   '],transfer=['MWC349  ','2013+370'])
calphase3mmtime=time.time()
#MWC249 = 0.82 Jy
#2013+370 = 2.60 Jy

print '--Set fluxscale (setjy)--'
default('setjy')
setjy('l2.ms',fieldid=1,spwid=6,fluxdensity=[1.01,0.,0.,0.])
setjy('l2.ms',fieldid=1,spwid=7,fluxdensity=[1.01,0.,0.,0.])
setjy('l2.ms',fieldid=1,spwid=14,fluxdensity=[1.01,0.,0.,0.])
setjy('l2.ms',fieldid=1,spwid=15,fluxdensity=[1.01,0.,0.,0.])
setjy('l2.ms',fieldid=1,spwid=18,fluxdensity=[1.01,0.,0.,0.])
setjy('l2.ms',fieldid=1,spwid=19,fluxdensity=[1.01,0.,0.,0.])
setjy('l2.ms',fieldid=2,spwid=6,fluxdensity=[2.79,0.,0.,0.])
setjy('l2.ms',fieldid=2,spwid=7,fluxdensity=[2.79,0.,0.,0.])
setjy('l2.ms',fieldid=2,spwid=14,fluxdensity=[2.79,0.,0.,0.])
setjy('l2.ms',fieldid=2,spwid=15,fluxdensity=[2.79,0.,0.,0.])
setjy('l2.ms',fieldid=2,spwid=18,fluxdensity=[2.79,0.,0.,0.])
setjy('l2.ms',fieldid=2,spwid=19,fluxdensity=[2.79,0.,0.,0.])
setjy3mmtime=time.time()

## Amplitude calibration of 3mm continuum:
##
##  phase solutions will be pre-applied as well as carried forward 
##   to the output solution table.
print '--gaincal amp (3mm)--'
default('gaincal')
gaincal('l2.ms','l2.3mmcont.amp.gcal','none',msselect='FIELD_ID IN [0,1,2] && DATA_DESC_ID IN [7,15,19,6,14,18]',gaintype='GSPLINE',calmode='a',splinetime=10000.,refant='1',phasewrap=260,gaincurve=False,opacity=False,bptable='l2.3mmcont.bpoly',preavg=0.,gaintable='l2.3mmcont.gcal')
calamp3mmtime=time.time()

## Correct the target source C34S (J=2-1)
##
##  note that only the 60 central channels will be calibrated
##   since the BPOLY solution is only defined for these
default('correct')
correct('l2.ms',msselect='FIELD_ID==3 && DATA_DESC_ID==3',gaincurve=False,opacity=False,gaintable='l2.3mmcont.gcal',bptable='l2.3mmC34S.bpoly',pointtable='l2.3mmcont.amp.gcal')
# Correct Target other 3mm line data (CH3OH)
default('correct')
correct('l2.ms',msselect='FIELD_ID==3 && DATA_DESC_ID==11',gaincurve=False,opacity=False,gaintable='l2.3mmcont.gcal',bptable='l2.3mmch3oh.bpoly',pointtable='l2.3mmcont.amp.gcal')
# Correct Target source 3mm continuum 
default('correct')
correct('l2.ms',msselect='FIELD_ID==3 && DATA_DESC_ID IN [7,15,19,6,14,18]',gaincurve=False,opacity=False,gaintable='l2.3mmcont.gcal',bptable='l2.3mmcont.bpoly',pointtable='l2.3mmcont.amp.gcal')
correct3mmtime=time.time()

# Split calibrated data
print '--split calibrater--'
default('split')
split('l2.ms','l2.3mmconta.split.ms',fieldid=3,spwid=[6,7,14,15,18,19],nchan=1,start=10,step=100,datacolumn='corrected')
splitsrctime=time.time()

#Image target source in 3mm continuum emission
default('clean')
clean('l2.3mmconta.split.ms',imagename='l2.3mmcont',alg='clark',niter=100,gain=0.1,nchan=1,start=0,width=1,spwid=range(0,6),fieldid=0,stokes='I',weighting='natural',cell=[0.4,0.4],imsize=[256,256],mode='mfs')


# 1mm Calibration

###########################################################
## Get first cut phase solutions to improve S/N for BPass determination:
## 
default('gaincal')
gaincal('l2.ms','l2.1mmUSB.gcal0','channel',nchan=110,start=9,step=1,msselect='FIELD_ID==0 && DATA_DESC_ID IN [27,31]',gaintype='GSPLINE',calmode='p',splinetime=10000.,refant='1',phasewrap=260,preavg=0.,gaincurve=False,opacity=False)

## Derive bandpass calibration for 1mm 
##
default('bandpass')
bandpass('l2.ms','l2.1mmcont.bpoly','none',msselect='FIELD_ID==0 && DATA_DESC_ID IN [27,31]',bandtype='BPOLY',degamp=10,degphase=20,visnorm=False,bpnorm=False,maskcenter=2,maskedge=0,refant='1',gaintable='l2.1mmUSB.gcal0',gaincurve=False,opacity=False)
default('bandpass')
bandpass('l2.ms','l2.1mmcont.bpoly','none',msselect='FIELD_ID==0 && DATA_DESC_ID IN [26,30]',bandtype='BPOLY',degamp=10,degphase=20,visnorm=False,bpnorm=False,maskcenter=2,maskedge=0,refant='1',gaintable='l2.1mmUSB.gcal0',gaincurve=False,opacity=False,append=True)
default('bandpass')
bandpass('l2.ms','l2.1mmUSB.bpoly','none',msselect='FIELD_ID==0 && DATA_DESC_ID IN [23]',bandtype='BPOLY',degamp=10,degphase=20,visnorm=False,bpnorm=False,maskcenter=2,maskedge=0,refant='1',gaintable='l2.1mmUSB.gcal0',gaincurve=False,opacity=False)

## Determine phase solutions for 1mm LSB & USB 
##
print '--Phase solutions--'
default('gaincal')
gaincal('l2.ms','l2.1mmcont.gcal','none',msselect='FIELD_ID IN [0,1,2] && DATA_DESC_ID IN [26,27,30,31]',gaintype='GSPLINE',calmode='p',splinetime=10000.,refant='1',phasewrap=260,gaincurve=False,opacity=False,bptable='l2.1mmcont.bpoly',preavg=0.)

##  Apply all solutions derived so far, determine
##  calibrators' flux densities using a solve for T and
##  fluxscale
print '--Solve for solutions, fluxscale--'
default('gaincal')
gaincal('l2.ms','l2.1mmcont.temp','none',msselect='FIELD_ID IN [0,1,2] && DATA_DESC_ID IN [26,27,30,31]',solint=1200.,refant='1',gaintype='T',opacity=False,gaincurve=False,gaintable='l2.1mmcont.gcal',bptable='l2.1mmcont.bpoly')
#
default('fluxscale')
fluxscale('l2.ms','l2.1mmcont.temp','l2.1mmcont.flux',reference=['3C273   '],transfer=['MWC349  ','2013+370'])
calphase1mmtime=time.time()
# MWC349:   0.8x (AIPS++ 0.75)
# 2013+370: 1.3  (AIPS++: 1.3)

## Record flux values from logger window.  Manually insert
## fluxes with imgr.setjy:
print '--Setjy 1mm --'
setjy('l2.ms',fieldid=1,spwid=27,fluxdensity=[1.75,0.,0.,0.])
setjy('l2.ms',fieldid=1,spwid=31,fluxdensity=[1.75,0.,0.,0.])
setjy('l2.ms',fieldid=2,spwid=27,fluxdensity=[1.8,0.,0.,0.])
setjy('l2.ms',fieldid=2,spwid=31,fluxdensity=[1.8,0.,0.,0.])
setjy('l2.ms',fieldid=1,spwid=26,fluxdensity=[1.75,0.,0.,0.])
setjy('l2.ms',fieldid=1,spwid=30,fluxdensity=[1.75,0.,0.,0.])
setjy('l2.ms',fieldid=2,spwid=26,fluxdensity=[1.8,0.,0.,0.])
setjy('l2.ms',fieldid=2,spwid=30,fluxdensity=[1.8,0.,0.,0.])
setjy1mmtime=time.time()


## Amplitude calibration of 1mm LSB:
default('gaincal')
gaincal('l2.ms','l2.1mmcont.amp.gcal','none',msselect='FIELD_ID IN [0,1,2] && DATA_DESC_ID IN [26,27,30,31]',gaintype='GSPLINE',calmode='a',splinetime=10000.,refant='1',phasewrap=260,gaincurve=False,opacity=False,bptable='l2.1mmcont.bpoly',preavg=0.,gaintable='l2.1mmcont.gcal')
calamp1mmtime=time.time()

## Correct the target source 1mm LSB & USB continuum data:
##
##  Note that edge channels in won't be calibrated because
##   BPOLY solution is not defined for them
default('correct')
correct('l2.ms',msselect='FIELD_ID==3 && DATA_DESC_ID IN [26,27,30,31]',gaincurve=False,opacity=False,gaintable='l2.1mmcont.gcal',bptable='l2.1mmcont.bpoly',pointtable='l2.1mmcont.amp.gcal')
default('correct')
correct('l2.ms',msselect='FIELD_ID==3 && DATA_DESC_ID IN [23]',gaincurve=False,opacity=False,gaintable='l2.1mmcont.gcal',bptable='l2.1mmcont.bpoly',pointtable='l2.1mmcont.amp.gcal')
correct1mmtime=time.time()

## Split out calibrated target source 1 mm continuum data:
default('split')
split('l2.ms','l2.1mma.split.ms',fieldid=3,spwid=[26,27,30,31],nchan=110,start=8,step=1,datacolumn='corrected')
default('split')
splitsrc1mmtime=time.time()

print '--Image 1mm continuum--'
## Get a first image the target source in 1 mm continuum emission:
default('clean')
clean('l2.1mma.split.ms',imagename='l2.1mm',alg='clark',niter=100,gain=0.1,nchan=1,start=0,width=1,spwid=[0,1,2,3],fieldid=0,stokes='I',weighting='natural',cell=[0.2,0.2],imsize=[128,128],mode='mfs')
image1mmtime=time.time()

print '--Resplit data--'
default('split')
split('l2.ms','l2.c34s.split.ms',fieldid=3,spwid=3,nchan=512,start=0,step=1,datacolumn='corrected')
default('split')
split('l2.ms','l2.ch3oh.split.ms',fieldid=3,spwid=[11],nchan=512,start=0,step=1,datacolumn='corrected')
default('split')
split('l2.ms','l2.3mmcont.split.ms',fieldid=3,spwid=[6,7,14,15,18,19],nchan=1,start=9,step=100,datacolumn='corrected')
default('split')
split('l2.ms','l2.1mmcont.split.ms',fieldid=3,spwid=[26,27,30,31],nchan=1,start=9,step=110,datacolumn='corrected')
default('split')
split('l2.ms','l2.1mmc34s5_4.split.ms',fieldid=3,spwid=[23],nchan=512,start=0,step=1,datacolumn='corrected')
splitsrc1mmtime=time.time()

### COMBINE DATA SETS

print '--Concatenate the data sets--'
os.system('cp -r l1.3mmcont.split.ms l02d.3mmcont.split.ms')
os.system('cp -r l1.c34s.split.ms l02d.c34s.split.ms')
os.system('cp -r l1.ch3oh.split.ms l02d.ch3oh.split.ms')
os.system('cp -r l1.1mmcont.split.ms l02d.1mmcont.split.ms')
default('concat')
concat('l02d.3mmcont.split.ms','l2.3mmcont.split.ms',freqtol='10MHz',dirtol='0.1arcsec')
concat('l02d.1mmcont.split.ms','l2.1mmcont.split.ms',freqtol='10MHz',dirtol='0.1arcsec')
concat('l02d.c34s.split.ms','l2.c34s.split.ms',freqtol='10MHz',dirtol='0.1arcsec')
concat('l02d.ch3oh.split.ms','l2.ch3oh.split.ms',freqtol='10MHz',dirtol='0.1arcsec')

#image target source in 3mm continuum emission
default('clean')
clean('l02d.3mmcont.split.ms',imagename='l02d.3mmcont',alg='clark',niter=1000,gain=0.1,nchan=1,start=0,width=1,spwid=range(0,6),fieldid=0,stokes='I',weighting='natural',cell=[0.6,0.6],imsize=[128,128],mode='mfs',threshold=0.7)

#image target source in 1mm continuum emission
default('clean')
clean('l02d.1mmcont.split.ms',imagename='l02d.1mmcont',alg='clark',niter=1000,gain=0.1,nchan=1,start=0,width=1,spwid=[0,1,2,3],fieldid=0,stokes='I',weighting='natural',cell=[0.2,0.2],imsize=[128,128],mode='mfs',threshold=4.)

#image target source in 3mm c34s (2-1) line
default('clean')
clean('l02d.c34s.split.ms',imagename='l02d.c34s',alg='clark',niter=1000,gain=0.1,nchan=160,start=19,width=3,spwid=0,fieldid=0,stokes='I',weighting='natural',cell=[0.6,0.6],imsize=[128,128],mode='channel',threshold=25.)

#image target source in 3mm ch3oh 
default('clean')
clean('l02d.ch3oh.split.ms',imagename='l02d.ch3oh',alg='clark',niter=1000,gain=0.1,nchan=160,start=19,width=3,spwid=0,fieldid=0,stokes='I',weighting='natural',cell=[0.6,0.6],imsize=[128,128],mode='channel',threshold=18.)

endProc=time.clock()
endTime=time.time()

# Regression

ms.open('l02d.1mmcont.split.ms')
thistest_1mm=max(ms.range(['amplitude']).get('amplitude'))
ms.close()
ms.open('l02d.3mmcont.split.ms')
thistest_3mm=max(ms.range(['amplitude']).get('amplitude'))
ms.close()
ms.open('l02d.ch3oh.split.ms')
thistest_ch3oh=max(ms.range(['amplitude']).get('amplitude'))
ms.close()
ms.open('l02d.c34s.split.ms')
thistest_c34s=max(ms.range(['amplitude']).get('amplitude'))
ms.close()
ia.open('l02d.3mmcont.image')
statistics=ia.statistics()
cont3mmmax=statistics['max'][0]
cont3mmrms=statistics['rms'][0]
ia.close()
ia.open('l02d.1mmcont.image')
statistics=ia.statistics()
cont1mmmax=statistics['max'][0]
cont1mmrms=statistics['rms'][0]
ia.close()
ia.open('l02d.c34s.image')
statistics=ia.statistics()
c34smax=statistics['max'][0]
c34srms=statistics['rms'][0]
ia.close()
ia.open('l02d.ch3oh.image')
statistics=ia.statistics()
ch3ohmax=statistics['max'][0]
ch3ohrms=statistics['rms'][0]
ia.close()

src3mm=0.471
src1mm=2.46
srcch3oh=419.210
srcc34s=318.025
im3mm=0.020
im1mm=0.136
imch3oh=0.464
imc34s=0.167

diff_3mm=abs((src3mm-thistest_3mm)/src3mm)
diff_1mm=abs((src1mm-thistest_1mm)/src1mm)
diff_ch3oh=abs((srcch3oh-thistest_ch3oh)/srcch3oh)
diff_c34s=abs((srcc34s-thistest_c34s)/srcc34s)
diff_im3mm=abs((im3mm-cont3mmmax)/im3mm)
diff_im1mm=abs((im1mm-cont1mmmax)/im1mm)
diff_imch3oh=abs((imch3oh-ch3ohmax)/imch3oh)
diff_imc34s=abs((imc34s-c34smax)/imc34s)

import datetime
datestring=datetime.datetime.isoformat(datetime.datetime.today())
outfile='l02d.'+datestring+'.log'
logfile=open(outfile,'w')

print >>logfile,'********** Regression ***********'
print >>logfile,'*                               *'
if (diff_3mm<0.05): print >>logfile,'* Passed 3mm max amplitude test '
print >>logfile,'--3mm max amp '+str(thistest_3mm)
if (diff_1mm<0.05): print >>logfile,'* Passed 1mm max amplitude test '
print >>logfile,'--1mm max amp '+str(thistest_1mm)
if (diff_ch3oh<0.05): print >>logfile,'* Passed CH3OH max amplitude test '
print >>logfile,'--CH3OH max amp '+str(thistest_ch3oh)
if (diff_c34s<0.05): print >>logfile,'* Passed C34S max amplitude test '
print >>logfile,'--C34S max amp '+str(thistest_c34s)
if (diff_im3mm<0.05): print >>logfile,'* Passed image (3mm) max test'
print >>logfile,'--Image max (3mm) '+str(cont3mmmax)
if (diff_im1mm<0.05): print >>logfile,'* Passed image (1mm) max test'
print >>logfile,'--Image max (1mm) '+str(cont1mmmax)
if (diff_imch3oh<0.05): print >>logfile,'* Passed image (ch3oh) max test'
print >>logfile,'--Image max (CH3OH) '+str(ch3ohmax)
if (diff_imc34s<0.05): print >>logfile,'* Passed image (c34s) max test'
print >>logfile,'--Image max (C34S) '+str(c34smax)


if ((diff_3mm<0.05) & (diff_1mm<0.05) & (diff_ch3oh<0.05) & (diff_ch3oh<0.05) & (diff_im1mm<0.05) & (diff_im3mm<0.05) & (diff_imch3oh < 0.05) & (diff_imc34s<0.05)):
	regstate=True
        print >>logfile,'---'
        print >>logfile,'Passed Regression test for L02D'
        print >>logfile,'---'
else:
        regstate=False
        print >>logfile,'----FAILED Regression test for L02D'
print >>logfile,'*********************************'

print >>logfile,''
print >>logfile,''
print >>logfile,'********* Benchmarking *****************'
print >>logfile,'*                                      *'
print >>logfile,'Total wall clock time was: '+str(endTime - startTime)
print >>logfile,'Total CPU        time was: '+str(endProc - startProc)
print >>logfile,'Processing rate MB/s  was: '+str(278./(endTime - startTime))
#print >>logfile,'* Breakdown:                           *'
#print >>logfile,'*   copy         time was: '+str(copytime-startTime)
#print >>logfile,'*   flag         time was: '+str(flagtime-copytime)
#print >>logfile,'*   cal ph 3mm   time was: '+str(calphase3mmtime-flagtime)
#print >>logfile,'*   setjy 3mm    time was: '+str(setjy3mmtime-calphase3mmtime)
#print >>logfile,'*   cal amp 3mm  time was: '+str(calamp3mmtime-setjy3mmtime)
#print >>logfile,'*   correct 3mm  time was: '+str(correct3mmtime-calamp3mmtime)
#print >>logfile,'*   split cal    time was: '+str(splitcaltime-correct3mmtime)
#print >>logfile,'*   split src    time was: '+str(splitsrctime-splitcaltime)
#print >>logfile,'*   image src    time was: '+str(image3mmtime-splitsrctime)
#print >>logfile,'*   cal ph 1mm   time was: '+str(calphase1mmtime-image3mmtime)
#print >>logfile,'*   setjy 1mm    time was: '+str(setjy1mmtime-calphase1mmtime)
#print >>logfile,'*   cal amp 1mm  time was: '+str(calamp1mmtime-setjy1mmtime)
#print >>logfile,'*   correct 1mm  time was: '+str(correct1mmtime-calamp1mmtime)
#print >>logfile,'*   splitsrc 1mm time was: '+str(splitsrc1mmtime-correct1mmtime)
#print >>logfile,'*   image 1mm    time was: '+str(image1mmtime-splitsrc1mmtime)
#print >>logfile,'*   HCO cal      time was: '+str(hcocaltime-image1mmtime)
#print >>logfile,'*   image HCO    time was: '+str(imagehcotime-hcocaltime)
#print >>logfile,'*   CO cal       time was: '+str(cocaltime-imagehcotime)
#print >>logfile,'*   image CO     time was: '+str(imagecotime-cocaltime)

logfile.close()