###############################################
#                                             #
# Regression/Benchmarking Script for GG TAU   #
#                                             #
###############################################
import time
import os
# Data must be pre-filled
# ggtau_07feb97.ms

pathname=os.environ.get('AIPSPATH')
datapath=pathname+'/data/regression/ATST1/GGTAU/07feb97-g067.ms'

restore()

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

print '--Copy MS--'
# Fill - really copy over pristine MS
os.system('rm -rf ggtau.3mm* ggtau.1mm* ggtau.hco.* ggtau.co.* ggtau_07feb97.ms')
#os.system('cp -r 07feb97-g067.ms ggtau_07feb97.ms')
copystring='cp -r '+datapath+' ggtau_07feb97.ms'
os.system(copystring)
clearcal('ggtau_07feb97.ms')
copytime=time.time()

print '--Flag data--'
# flag three end channels for spectral window 2, fieldid=2
#flagdata('ggtau_07feb97.ms',fieldid=[2],spwid=2,chans=[0,1,2,61,62,63])
# flag PdBI central data channels (32,33) that suffer ringing due to Gibbs phenomenon
default('flagdata')
flagdata('ggtau_07feb97.ms',fieldid=[0,1,3],spwid=2,chans=[31,32])
# flag channels 33,34 which are low for source 0415+379 (fieldid=0)
default('flagdata')
flagdata('ggtau_07feb97.ms',fieldid=[0],spwid=2,chans=[33,34])
# flag channels 34-37 for CRL 618 (fieldid=3) and baseline 1-3,2-4
default('flagdata')
flagdata('ggtau_07feb97.ms',fieldid=[3],baseline=[0,2],chans=[34,37])
default('flagdata')
flagdata('ggtau_07feb97.ms',fieldid=[3],baseline=[1,3],chans=[34,37])
# flag channels 250-256 for spectral window 7
#default('flagdata')
#flagdata('ggtau_07feb97.ms',spwid=[6])
flagtime=time.time()

#
# 3 mm Continuum calibration
#
print '--Continuum calibration phase/bandpass (3mm)--'
#setjy
default('setjy')
setjy('ggtau_07feb97.ms',fieldid=3,spwid=2,fluxdensity=[1.55,0.,0.,0.])

#preliminary time-dependent phase solutions to improve coherent average for bandpass solution
default('gaincal')
gaincal('ggtau_07feb97.ms','ggtau.3mm.ph.gcal0','channel',nchan=30,start=14,step=1,msselect='FIELD_ID==1 && DATA_DESC_ID==2',gaintype='GSPLINE',calmode='p',splinetime=10000.,refant='1',phasewrap=260,gaincurve=False,opacity=False,preavg=120)

print '--bandpass (3mm)--'
#derive bandpass calibration for 3mm LSB
default('bandpass')
bandpass('ggtau_07feb97.ms','ggtau.3mm.bpoly','channel',nchan=60,start=3,step=1,msselect='FIELD_ID==1 && DATA_DESC_ID==2',bandtype='BPOLY',degamp=6,degphase=12,bpnorm=False,maskcenter=4,maskedge=0,refant='1',gaintable='ggtau.3mm.ph.gcal0',gaincurve=False,opacity=False)

print '--gaincal phase (3mm)--'
#derive new and better phase solutions for 3mm LSB
default('gaincal')
gaincal('ggtau_07feb97.ms','ggtau.3mm.ph.gcal','channel',nchan=60,start=2,step=1,msselect='FIELD_ID IN [0,1,3] && DATA_DESC_ID==2',gaintype='GSPLINE',calmode='p',splinetime=10000.,refant='1',phasewrap=260,gaincurve=False,opacity=False,bptable='ggtau.3mm.bpoly',preavg=0.)

#Apply all solutions derived so far, determine calibrator's flux densities by solving for T and using fluxscale
default('gaincal')
gaincal('ggtau_07feb97.ms','ggtau.3mm.temp','channel',nchan=60,start=2,step=1,msselect='FIELD_ID IN [0,1,3] && DATA_DESC_ID==2',solint=600.,refant='1',gaintype='T',opacity=False,gaincurve=False,gaintable='ggtau.3mm.ph.gcal',bptable='ggtau.3mm.bpoly')

#fluxscale
default('fluxscale')
fluxscale('ggtau_07feb97.ms','ggtau.3mm.temp','ggtau.3mm.flux',reference=['CRL618  '],transfer=['0415+379','0528+134'])
### Flux density for 0415+379 (spw=3) is:     5.74898 +/-  0.0267088 Jy
### Flux density for 0528+134 (spw=3) is:     2.56862 +/-  0.011453  Jy
calphase3mmtime=time.time()

print '--Set fluxscale (setjy)--'
default('setjy')
setjy('ggtau_07feb97.ms',fieldid=0,spwid=2,fluxdensity=[5.849,0.,0.,0.])
setjy('ggtau_07feb97.ms',fieldid=1,spwid=2,fluxdensity=[2.634,0.,0.,0.])
setjy3mmtime=time.time()

## Amplitude calibration of 3mm LSB:
##
##  phase solutions will be pre-applied as well as carried forward 
##   to the output solution table.
print '--gaincal amp (3mm)--'
default('gaincal')
gaincal('ggtau_07feb97.ms','ggtau.3mm.amp.gcal','channel',nchan=60,start=2,step=1,msselect='FIELD_ID IN [0,1,3] && DATA_DESC_ID==2',gaintype='GSPLINE',calmode='a',splinetime=20000.,refant='1',phasewrap=260,gaincurve=False,opacity=False,bptable='ggtau.3mm.bpoly',preavg=2500.,pointtable='ggtau.3mm.ph.gcal')
calamp3mmtime=time.time()

## Correct the target source and all other 3mm LSB data: 
##
##  note that only the 60 central channels will be calibrated
##   since the BPOLY solution is only defined for these
default('correct')
correct('ggtau_07feb97.ms',msselect='DATA_DESC_ID==2',gaincurve=False,opacity=False,gaintable='ggtau.3mm.ph.gcal',bptable='ggtau.3mm.bpoly',pointtable='ggtau.3mm.amp.gcal')
correct3mmtime=time.time()

# Split calibrater data
print '--split calibrater--'
default('split')
split('ggtau_07feb97.ms','ggtau.cal.split.ms',fieldid=3,spwid=2,nchan=64,start=0,step=1,datacolumn='model')
splitcaltime=time.time()

#Split calibrated target data
print '--split source--'
default('split')
split('ggtau_07feb97.ms','ggtau.3mm.split.ms',fieldid=2,spwid=2,nchan=60,start=2,step=1,datacolumn='corrected')
splitsrctime=time.time()

#First image the target source in 3mm continuum emission
print '--Image continuum (3mm)--'
default('clean')
clean('ggtau.3mm.split.ms',imagename='ggtau.3mm',alg='clark',niter=100,gain=0.1,nchan=1,start=2,width=58,spwid=0,fieldid=0,stokes='I',weighting='briggs',rmode='norm',robust=0.5,cell=[0.2,0.2],imsize=[256,256])
image3mmtime=time.time()

# 1mm Calibration

print '--Continuum calibration phase/bandpass (1mm)--'
default('setjy')
setjy('ggtau_07feb97.ms',fieldid=3,spwid=10,fluxdensity=[2.2,0.,0.,0.])
setjy('ggtau_07feb97.ms',fieldid=3,spwid=14,fluxdensity=[2.2,0.,0.,0.])
setjy('ggtau_07feb97.ms',fieldid=3,spwid=18,fluxdensity=[2.2,0.,0.,0.])

   ### CRL618    spwid= 11  [I=2.2, Q=0, U=0, V=0] Jy, (user-specified)
   ### CRL618    spwid= 15  [I=2.2, Q=0, U=0, V=0] Jy, (user-specified)
   ### CRL618    spwid= 19  [I=2.2, Q=0, U=0, V=0] Jy, (user-specified)

###########################################################
## Get first cut phase solutions to improve S/N for BPass determination:
## 
default('gaincal')
gaincal('ggtau_07feb97.ms','ggtau.1mm.ph.gcal0','channel',nchan=30,start=14,step=1,msselect='FIELD_ID==1 && DATA_DESC_ID IN [10,14,18]',gaintype='GSPLINE',calmode='p',splinetime=10000.,refant='1',phasewrap=260,preavg=120,gaincurve=False,opacity=False)

## Derive bandpass calibration for 1mm LSB:
##
##  Avoid edge channels in each spw using setdata, not maskedge
print '--Bandpass solution (1mm)--'
default('bandpass')
bandpass('ggtau_07feb97.ms','ggtau.1mm.bpoly','channel',nchan=60,start=2,step=1,msselect='FIELD_ID==1 && DATA_DESC_ID IN [10,14,18]',bandtype='BPOLY',degamp=6,degphase=6,visnorm=False,bpnorm=False,maskcenter=4,maskedge=0,refant='1',gaintable='ggtau.1mm.ph.gcal0',gaincurve=False,opacity=False)

## Determine phase solutions for 1mm LSB for all calibrators
##
##  Note, we are using the 3mm phase solution, transferred
##   to the 1mm bands, to reduce the phase rate.  The
##   net solution (transferred 3mm + 1mm) is written
##   to the output table
print '--Phase solutions--'
default('gaincal')
gaincal('ggtau_07feb97.ms','ggtau.1mm.ph.gcal','channel',nchan=60,start=2,step=1,msselect='FIELD_ID IN [0,1,3] && DATA_DESC_ID IN [10,14,18]',gaintype='GSPLINE',calmode='p',splinetime=10000.,refant='1',phasewrap=260,gaincurve=False,opacity=False,bptable='ggtau.1mm.bpoly',gaintable='ggtau.3mm.ph.gcal',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('ggtau_07feb97.ms','ggtau.1mm.temp','channel',nchan=60,start=2,step=1,msselect='FIELD_ID IN [0,1,3] && DATA_DESC_ID IN [10,14,18]',solint=600.,refant='1',gaintype='T',opacity=False,gaincurve=False,gaintable='ggtau.1mm.ph.gcal',bptable='ggtau.1mm.bpoly')
#
default('fluxscale')
fluxscale('ggtau_07feb97.ms','ggtau.1mm.temp','ggtau.1mm.flux',reference=['CRL618  '],transfer=['0415+379','0528+134'])
calphase1mmtime=time.time()

## Record flux values from logger window.  Manually insert
## fluxes with imgr.setjy:
print '--Setjy 1mm --'
setjy('ggtau_07feb97.ms',fieldid=0,spwid=10,fluxdensity=[4.310,0.,0.,0.])
setjy('ggtau_07feb97.ms',fieldid=0,spwid=14,fluxdensity=[4.310,0.,0.,0.])
setjy('ggtau_07feb97.ms',fieldid=0,spwid=18,fluxdensity=[4.310,0.,0.,0.])
setjy('ggtau_07feb97.ms',fieldid=1,spwid=10,fluxdensity=[1.842,0.,0.,0.])
setjy('ggtau_07feb97.ms',fieldid=1,spwid=14,fluxdensity=[1.842,0.,0.,0.])
setjy('ggtau_07feb97.ms',fieldid=1,spwid=18,fluxdensity=[1.842,0.,0.,0.])
setjy1mmtime=time.time()


## Amplitude calibration of 1mm LSB:
default('gaincal')
gaincal('ggtau_07feb97.ms','ggtau.1mm.amp.gcal','channel',nchan=60,start=2,step=1,msselect='FIELD_ID IN [0,1,3] && DATA_DESC_ID IN [10,14,18]',gaintype='GSPLINE',calmode='a',splinetime=20000.,refant='1',phasewrap=260,gaincurve=False,opacity=False,bptable='ggtau.1mm.bpoly',preavg=2500.,gaintable='ggtau.1mm.ph.gcal')
calamp1mmtime=time.time()

## Correct the target source and all other 1mm LSB continuum data:
##
##  Note that edge channels in won't be calibrated because
##   BPOLY solution is not defined for them
default('correct')
correct('ggtau_07feb97.ms',msselect='DATA_DESC_ID IN [10,14,18]',gaincurve=False,opacity=False,gaintable='ggtau.1mm.ph.gcal',bptable='ggtau.1mm.bpoly',pointtable='ggtau.1mm.amp.gcal')
correct1mmtime=time.time()

## Split out calibrated target source 1 mm continuum data:
split('ggtau_07feb97.ms','ggtau.1mm.split.ms',fieldid=2,spwid=[10,14,18],nchan=60,start=2,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('ggtau.1mm.split.ms',imagename='ggtau.1mm',alg='clark',niter=100,gain=0.1,nchan=1,start=3,width=58,spwid=[0,1,2],fieldid=0,stokes='I',weighting='briggs',rmode='norm',robust=0.5,cell=[0.1,0.1],imsize=[256,256])
image1mmtime=time.time()

# 3mm HCO+(1-0) line calibration
 ### Note: during the continuum reduction above, bandpass and time-dependent
 ### phase and amplitude solutions were derived using the wider bandwidth 3 mm LSB
 ### continuum data. These can be applied directly to the line data.
 ### Since the bandpass solution for the wide, low-resolution continuum band
 ### may not be ideal for the narrow, high-resolution line band.  If there
 ### is sufficient S/R on the bandpass calibrator in the narrow line band,
 ### a bandpass solution may be obtained for it as follows:

print '--Cal HCO and split--'
default('bandpass')
bandpass('ggtau_07feb97.ms','ggtau.hco.bpoly','channel',nchan=206,start=25,step=1,msselect='FIELD_ID==1 && DATA_DESC_ID==6',bandtype='BPOLY',degamp=1,degphase=1,visnorm=False,bpnorm=False,maskcenter=40,maskedge=0,refant='1',gaintable='ggtau.3mm.ph.gcal0',gaincurve=False,opacity=False)

default('correct')
correct('ggtau_07feb97.ms',msselect='DATA_DESC_ID==6',gaincurve=False,opacity=False,gaintable='ggtau.3mm.ph.gcal',bptable='ggtau.hco.bpoly',pointtable='ggtau.3mm.amp.gcal')
default('split')
split('ggtau_07feb97.ms','ggtau.hco.split.ms',fieldid=2,spwid=6,nchan=206,start=25,step=1,datacolumn='corrected')
hcocaltime=time.time()

print 'Image HCO--'
default('clean')
clean('ggtau.hco.split.ms',imagename='ggtau.hco',alg='clark',niter=100,gain=0.1,nchan=14,start=74,width=4,spwid=0,fieldid=0,stokes='I',weighting='briggs',rmode='norm',robust=0.5,cell=[0.2,0.2],imsize=[256,256],mode='channel')
   ### Emission in central channels is barely visible (too weak to image
   ###  with only 1 day of data but you get the point).  
imagehcotime=time.time()

# 1mm 13CO(2-1) line calibration
 ### Note: during the continuum reduction above, bandpass and time-dependent
 ### phase and amplitude solutions were derived using the wider bandwidth 1 mm LSB
 ### continuum data. These can be applied directly to the line data.
 ### The bandpass solution for the wide, low-resolution continuum band 
 ### may not be ideal for the narrow, high-resolution line band.  If there
 ### is sufficient S/R on the bandpass calibrator in the narrow line band,
 ### a bandpass solution may be obtained for it as follows:
print '--Cal CO and split--'
default('bandpass')
bandpass('ggtau_07feb97.ms','ggtau.co.bpoly','channel',nchan=240,start=8,step=1,msselect='FIELD_ID==1 && DATA_DESC_ID==22',bandtype='BPOLY',degamp=1,degphase=1,visnorm=False,bpnorm=False,maskcenter=4,maskedge=0,refant='1',gaintable='ggtau.1mm.ph.gcal0',gaincurve=False,opacity=False)
default('correct')
correct('ggtau_07feb97.ms',msselect='DATA_DESC_ID==22',gaincurve=False,opacity=False,gaintable='ggtau.1mm.ph.gcal',bptable='ggtau.co.bpoly',pointtable='ggtau.1mm.amp.gcal')
default('split')
split('ggtau_07feb97.ms','ggtau.co.split.ms',fieldid=2,spwid=22,nchan=240,start=0,step=1,datacolumn='corrected')
cocaltime=time.time()

print '--Image CO--'
default('clean')
clean('ggtau.co.split.ms',imagename='ggtau.co',alg='clark',niter=100,gain=0.1,nchan=14,start=91,width=4,spwid=0,fieldid=0,stokes='I',weighting='briggs',rmode='norm',robust=0.5,cell=[0.1,0.1],imsize=[256,256])
   ### Emission is too weak to image with only 1 day of data.
imagecotime=time.time()

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

#Regression

ms.open('ggtau.cal.split.ms')
thistest_cal=max(ms.range(["amplitude"]).get('amplitude'))
ms.close()
ms.open('ggtau.3mm.split.ms')
thistest_3mm=max(ms.range(["amplitude"]).get('amplitude'))
ms.close()
ms.open('ggtau.1mm.split.ms')
thistest_1mm=max(ms.range(["amplitude"]).get('amplitude'))
ms.close()
ms.open('ggtau.hco.split.ms')
thistest_hco=max(ms.range(["amplitude"]).get('amplitude'))
ms.close()
ms.open('ggtau.co.split.ms')
thistest_co=max(ms.range(["amplitude"]).get('amplitude'))
ms.close()
ia.open('ggtau.3mm.image')
statistics=ia.statistics()
thistest_immax=statistics['max'][0]
thistest_imrms=statistics['rms'][0]

calmax=1.55
#src3mmmax=0.89268
src3mmmax=0.96
src1mmmax=17.32818
srchcomax=6.853867
srcco21max=52.9783
imgmax=0.0055763
imgrms=0.000717

diff_cal=abs((calmax-thistest_cal)/calmax)
diff_3mm=abs((src3mmmax-thistest_3mm)/src3mmmax)
diff_1mm=abs((src1mmmax-thistest_1mm)/src1mmmax)
diff_hco=abs((srchcomax-thistest_hco)/srchcomax)
diff_co=abs((srcco21max-thistest_co)/srcco21max)
diff_immax=abs((imgmax-thistest_immax)/imgmax)
diff_imrms=abs((imgrms-thistest_imrms)/imgrms)

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

print >>logfile,'MeasurementSet Name: /home/rohir2/jmcmulli/ALMATST1/Regression/GGTAU/ggtau_07feb97.ms MS Version 2'
print >>logfile,''
print >>logfile,'   Observer:      Project: G067'
print >>logfile,'Observation: IRAM_PDB(5 antennas)'
print >>logfile,''
print >>logfile,'Data records: 86320       Total integration time = 33045 seconds'
print >>logfile,'   Observed from   14:50:27   to   00:01:12'
print >>logfile,''
print >>logfile,'Fields: 4'
print >>logfile,'  ID   Name          Right Ascension  Declination   Epoch'
print >>logfile,'  0    0415+379      04:18:21.28      +38.01.35.85  J2000'
print >>logfile,'  1    0528+134      05:30:56.41      +13.31.55.08  J2000'
print >>logfile,'  2    GG_TAU        04:32:30.34      +17.31.40.52  J2000'
print >>logfile,'  3    CRL618        04:42:53.56      +36.06.53.59  J2000'
print >>logfile,''
print >>logfile,'Spectral Windows:  (24 unique spectral windows and 1 unique polarization setups)'
print >>logfile,'  SpwID  #Chans Frame Ch1(MHz)    Resoln(kHz) TotBW(kHz)  Ref(MHz)    Corrs'
print >>logfile,'  0           1 LSRK  89188.5209  128000      128000      89188.5209  XX'
print >>logfile,'  1           1 LSRK  92264.5127  -128000     128000      92264.5127  XX'
print >>logfile,'  2          64 LSRK  89109.7709  2500        160000      89189.7709  XX'
print >>logfile,'  3          64 LSRK  92343.2627  -2500       160000      92263.2627  XX'
print >>logfile,'  4           1 LSRK  89188.5209  6000        6000        89188.5209  XX'
print >>logfile,'  5           1 LSRK  92264.5127  -6000       6000        92264.5127  XX'
print >>logfile,'  6         256 LSRK  89183.5404  39.0625     10000       89188.5404  XX'
print >>logfile,'  7         256 LSRK  92269.4932  -39.0625    10000       92264.4932  XX'
print >>logfile,'  8           1 LSRK  220253.658  128000      128000      220253.658  XX'
print >>logfile,'  9           1 LSRK  223632.358  -128000     128000      223632.358  XX'
print >>logfile,'  10         64 LSRK  220174.908  2500        160000      220254.908  XX'
print >>logfile,'  11         64 LSRK  223711.108  -2500       160000      223631.108  XX'
print >>logfile,'  12          1 LSRK  220398.681  128000      128000      220398.681  XX'
print >>logfile,'  13          1 LSRK  223487.336  -128000     128000      223487.336  XX'
print >>logfile,'  14         64 LSRK  220319.931  2500        160000      220399.931  XX'
print >>logfile,'  15         64 LSRK  223566.086  -2500       160000      223486.086  XX'
print >>logfile,'  16          1 LSRK  220543.703  128000      128000      220543.703  XX'
print >>logfile,'  17          1 LSRK  223342.314  -128000     128000      223342.314  XX'
print >>logfile,'  18         64 LSRK  220464.953  2500        160000      220544.953  XX'
print >>logfile,'  19         64 LSRK  223421.064  -2500       160000      223341.064  XX'
print >>logfile,'  20          1 LSRK  220398.681  16000       16000       220398.681  XX'
print >>logfile,'  21          1 LSRK  223487.336  -16000      16000       223487.336  XX'
print >>logfile,'  22        256 LSRK  220388.72   78.125      20000       220398.72   XX'
print >>logfile,'  23        256 LSRK  223497.297  -78.125     20000       223487.297  XX'
print >>logfile,''
print >>logfile,'Antennas: 5'
print >>logfile,'   ID=   1-4: 1=223, 2=323, 3=227, 4=316, 5=129'
print >>logfile,''
print >>logfile,'Tables(rows):   (-1 = table absent)'
print >>logfile,'   MAIN(86320)'
print >>logfile,'   ANTENNA(5)   DATA_DESCRIPTION(24)   DOPPLER(-1)   FEED(5)   FIELD(4)'
print >>logfile,'   FLAG_CMD(0)   FREQ_OFFSET(-1)   HISTORY(213)   OBSERVATION(1)   POINTING(0)'
print >>logfile,'   POLARIZATION(1)   PROCESSOR(0)   SOURCE(4)   SPECTRAL_WINDOW(24)   STATE(0)'
print >>logfile,'   SYSCAL(-1)   WEATHER(-1)'

print >>logfile,'********** Regression ***********'
print >>logfile,'*                               *'
if (diff_cal<0.05): print >>logfile,'* Passed cal max amplitude test '
print >>logfile,'--Cal max amp '+str(thistest_cal)
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_hco<0.05): print >>logfile,'* Passed HCO max amplitude test '
print >>logfile,'--HCO max amp '+str(thistest_hco)
if (diff_co<0.05): print >>logfile,'* Passed CO max amplitude test '
print >>logfile,'--CO max amp '+str(thistest_co)
if (diff_immax<0.05): print >>logfile,'* Passed image max test'
print >>logfile,'--Image max '+str(thistest_immax)
if (diff_imrms<0.05): print >>logfile,'* Passed image rms test'
print >>logfile,'--Image rms '+str(thistest_imrms)

if ((diff_cal<0.05) & (diff_3mm<0.05) & (diff_1mm<0.05) & (diff_hco<0.05) & (diff_co<0.05) & (diff_immax<0.05) & (diff_imrms<0.05)):
	regstate=True
        print >>logfile,'---'
        print >>logfile,'Passed Regression test for GG TAU'
        print >>logfile,'---'
else:
	regstate=False
        print >>logfile,'----FAILED Regression test for GG TAU'
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()