# H121

print '--Copy MS--'
# Fill - really copy over pristine MS
msfile2='09-3mmmay-1998-h121.ms'
clearcal(msfile2)

## Set flux density scale for the flux density calibrater MWC349
setjy(vis=msfile2,spw='7',field='4',fluxdensity=[1.12,0.,0.,0.])

## Get preliminary time-dependent phase solutions to improve
## coherent average for bandpass solution and derive bandpass
## calibration for 3mm LSB
default('gaincal')
gaincal(vis=msfile2,caltable='h121b.3mm.ph.gcal0',
	field='1',spw='7:8~55',
	gaintype='GSPLINE',calmode='p',splinetime=10000.,
	refant='1',phasewrap=260,
	gaincurve=False,opacity=0.0,preavg=120.)

print '--bandpass (3mm)--'
#derive bandpass calibration 
# solint value
default('bandpass')
bandpass(vis=msfile2,caltable='h121b.3mm.bpoly',
	 field='1',spw='7:8~55',
	 bandtype='BPOLY',degamp=6,degphase=12,solnorm=False,
	 maskcenter=6,maskedge=0,refant='1',
	 gaintable='h121.3mm.ph.gcal0',gaincurve=False,opacity=0.0,solint=0.)

# 3mm CO Band
default('gaincal')
gaincal(vis=msfile2,caltable='2200b.3mm.ph.gcal0',
	field='0',spw='3:64~191',
	gaintype='GSPLINE',calmode='p',splinetime=10000.,
	refant='1',preavg=10.,phasewrap=260,
	gaincurve=False,opacity=0.0)
default('bandpass')
bandpass(vis=msfile2,caltable='2200b.3mm.bpoly',
	 field='0',spw='3:8~247',
	 bandtype='BPOLY',degamp=6,degphase=12,solnorm=False,
	 maskcenter=6,maskedge=0,refant='1',
	 gaintable='2200b.3mm.ph.gcal0',gaincurve=False,opacity=0.0)

# Determine new and better phase solutions for 3mm LSB for all calibrators
default('gaincal')
gaincal(vis=msfile2,caltable='h121b.3mm.ph.gcal',
	field='0,1,3,4',spw='7:8~55',
	gaintype='GSPLINE',calmode='p',splinetime=5000.,refant='1',
	phasewrap=260,
	gaincurve=False,opacity=0.0,gaintable='h121b.3mm.bpoly',preavg=0.)

#Apply all solutions derived so far, determine calibrator's flux densities by solving for T 
#and use fluxscale
default('gaincal')
gaincal(vis=msfile2,caltable='h121b.3mm.temp',
	field='0,1,3,4',spw='7:2~61',
	solint=600.,refant='1',gaintype='T',
	opacity=0.0,gaincurve=False,
	gaintable=['h121b.3mm.ph.gcal','h121b.3mm.bpoly'])

#fluxscale
default('fluxscale')
fluxscale(vis=msfile2,caltable='h121b.3mm.temp',fluxtable='h121b.3mm.flux',
	  reference='MWC349*',transfer='3C454*,2200+420*,0224+671*')
calphase3mmtime=time.time()
#2200+420: 3.29
#0224+671: 1.53
#C454.3 : 6.127

#use old values
print '--Set fluxscale (setjy)--'
default('setjy')
setjy(vis=msfile2,field='0',spw='7',fluxdensity=[3.235,0.,0.,0.])
setjy(vis=msfile2,field='1',spw='7',fluxdensity=[1.507,0.,0.,0.])
setjy(vis=msfile2,field='3',spw='7',fluxdensity=[5.982,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(vis=msfile2,caltable='h121b.3mm.amp.gcal',
	field='0,1,3,4',spw='7',
	gaintype='GSPLINE',calmode='a',splinetime=20000.,refant='1',
	phasewrap=260,
	gaincurve=False,opacity=0.0,
	preavg=2500.,
	gaintable=['h121b.3mm.ph.gcal','h121b.3mm.bpoly'])
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('applycal')
applycal(vis=msfile2,
	 spw='7',
	 gaincurve=False,opacity=0.0,
	 gaintable=['h121b.3mm.ph.gcal','h121b.3mm.amp.gcal','h121b.3mm.bpoly'])
# Correct Target CO (1-0)
default('applycal')
applycal(vis=msfile2,
	 spw='3',
	 gaincurve=False,opacity=0.0,
	 gaintable=['h121b.3mm.ph.gcal','h121b.3mm.amp.gcal','2200b.3mm.bpoly'])

# Split calibrated target source data
print '--split calibrater--'
default('split')
split(vis=msfile2,outputvis='h121b.3mm.split.ms',
#      field=2,spw=7,nchan=48,start=8,step=1,datacolumn='corrected')
	field='2',spw='7:8~56',datacolumn='corrected')
default('split')
split(vis=msfile2,outputvis='0224b.3mm.split.ms',
#      field=1,spw=7,nchan=48,start=8,step=1,datacolumn='corrected')
	field='1',spw='7:8~56',datacolumn='corrected')
default('split')
split(vis=msfile2,outputvis='h121b.co10.split.ms',
#      field=2,spw=3,nchan=240,start=8,step=1,datacolumn='corrected')
	field='2',spw='3:8~247',datacolumn='corrected')	
splitsrctime=time.time()

## Get a first image the target source in 3 mm continuum emission:
default('clean')
clean(vis='h121b.co10.split.ms',imagename='h121b.co10',
      alg='clark',niter=1000,gain=0.1,threshold=20.,
      nchan=80,start=0,width=3,
      spw=0,field='0',stokes='I',
      weighting='briggs',rmode='norm',robust=2.,
      cell=[1.,1.],imsize=[64,64],mode='channel')
image1time=time.time()

#Image target in 3mm continuum emission
default('clean')
clean(vis='h121.3mm.concat.ms',imagename='h121all.3mm',
      alg='clark',niter=2000,gain=0.1,nchan=1,start=0,width=48,
      spw=0,field='0',stokes='I',
      weighting='briggs',rmode='norm',robust=2.,
      cell=[1.,1.],imsize=[64,64],mode='mfs')