#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#!!  This script describes SMA data reduction using miriad to fill and
#!! apply Tsys correction, and then writing uvfits. CASA is then used
#!! to read uvfits, flag, calibrate, continuum subtract, and image the
#!! data. There are several steps that are currently VERY slow for 
#!! datasets with many spw such as SMA data (plotxy, flagdata, 
#!! uvcontsub, scratch columns). Please be patient, CASA is a work in
#!! progress. This script requires Beta patch 2.4 or later.
#
#!! The data themselves are a six-pointing mosaic of the infrared dark
#!! cloud G19.3+0.07 with 12CO(2-1) in the USB from Brogan et al. in
#!! prep. Please do not use these data for scientific purposes.
#  
#!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#######################################################################
#######################################################################
# DATA LOADING, TSYS CORRECTION, AND UVFITS WRITING IN MIRIAD
#######################################################################
#######################################################################
# 
# smalod in='../050819_04.23.28' out=g19_d2 sideband=2 rxif=0
# 
# smafix vis=g19_d2_rx0.usb out=g19_d2_rx0.usb.tsys device=/xs xaxis=time 
#                  \  yaxis=systemp nxy=3,3 options=tsyscor
# 
# fits in=g19_d2_rx0.usb.tsys out=g19_d2usb.spw1 op=uvout select=win\(1\)
# fits in=g19_d2_rx0.usb.tsys out=g19_d2usb.spw2 op=uvout select=win\(2\)
# fits in=g19_d2_rx0.usb.tsys out=g19_d2usb.spw3 op=uvout select=win\(3\)
#   ...
#   ...
#   ...
# fits in=g19_d2_rx0.usb.tsys out=g19_d2usb.spw24 op=uvout select=win\(24\)
# 
#### Note that currently these really do have to be written out
#### individually (of course you could write a little miriad script to
#### do it). If not the spw are concatenated together into a single
#### spw which is inappropriate before calibration.
#
# mkdir spw_fits_files
# 
# mv g19_d2usb.spw* spw_fits_files/
#
# the file spw_fits_files.tar.gz is what is distributed as the "data" 
#
#######################################################################
#######################################################################
#### INSIDE CASAPY 
#######################################################################
#######################################################################

os.system('rm -rf spw_fits.files')
os.system('gunzip spw_fits_files.tar.gz')
os.system('tar xvf spw_fits_files.tar')
os.system('rm -rf spw_ms_files')
os.system('rm -rf g19_d2usb*')
os.system('mkdir spw_ms_files')


#### Read in each uvfits file (need to put upper range to one more than 
#### have, to satisfy ipython peculiarity). 

for i in range(1,25):
	importuvfits(fitsfile="spw_fits_files/g19_d2usb.spw"+str(i),
		vis="spw_ms_files/g19_d2usb.spw"+str(i)+".ms")

#### Create "dictionary" of filenames to be concatenated

import glob
myfiles = glob.glob('spw_ms_files/g19_d2usb.spw*.ms')

concat(vis=myfiles,concatvis='g19_d2usb.ms',freqtol="",
       dirtol="0.1arcsec",async=False)

#### Get a listing of the properties of the concatenated files. 

listobs(vis="g19_d2usb.ms",verbose=False)

#      0         jupiter       13:01:25.48      -05.18.53.00  J2000
#      1         1743-038      17:43:58.85      -03.50.04.62  J2000
#      2         g19.3a        18:25:58.70      -12.03.57.80  J2000
#      3         g19.3b        18:25:57.40      -12.04.18.50  J2000
#      4         g19.3c        18:25:56.09      -12.04.28.20  J2000
#      5         g19.3d        18:25:54.60      -12.04.23.10  J2000
#      6         1908-201      19:11:09.65      -20.06.55.11  J2000
#      7         g19.3e        18:25:54.80      -12.04.43.80  J2000
#      8         g19.3f        18:25:53.30      -12.04.51.00  J2000
#      9         3c454.3       22:53:57.74      +16.08.53.56  J2000
#      10        0359+509      03:59:29.75      +50.57.50.15  J2000
#      11        hip19167      04:06:35.04      +50.21.04.54  J2000
#      12        polaris       02:31:48.70      +89.15.50.72  J2000
#      13        uranus        22:44:40.51      -08.50.23.77  J2000
# 
#     Spectral Windows:  (24 unique spectral windows and 1 unique polarization setups)
#      SpwID  #Chans Frame Ch1(MHz)    ChanWid(kHz)TotBW(kHz)  Ref(MHz)    Corrs
#      0         128 LSRK  229347.598  812.5       104000      229347.598  XX
#      1         128 LSRK  229429.601  812.5       104000      229429.601  XX
#      2         128 LSRK  229504.79   812.5       104000      229504.79   XX
#      3         128 LSRK  229586.792  812.5       104000      229586.792  XX
#      4         128 LSRK  229675.607  812.5       104000      229675.607  XX
#      5         128 LSRK  229757.609  812.5       104000      229757.609  XX
#      6         128 LSRK  229832.798  812.5       104000      229832.798  XX
#      7         128 LSRK  229914.8    812.5       104000      229914.8    XX
#      8         128 LSRK  230003.615  812.5       104000      230003.615  XX
#      9         128 LSRK  230085.617  812.5       104000      230085.617  XX
#      10        128 LSRK  230160.806  812.5       104000      230160.806  XX
#      11        128 LSRK  230242.808  812.5       104000      230242.808  XX
#      12        128 LSRK  230331.623  812.5       104000      230331.623  XX
#      13        128 LSRK  230413.625  812.5       104000      230413.625  XX
#      14        128 LSRK  230488.814  812.5       104000      230488.814  XX
#      15        128 LSRK  230570.817  812.5       104000      230570.817  XX
#      16        128 LSRK  230659.631  812.5       104000      230659.631  XX
#      17        128 LSRK  230741.634  812.5       104000      230741.634  XX
#      18        128 LSRK  230816.823  812.5       104000      230816.823  XX
#      19        128 LSRK  230898.825  812.5       104000      230898.825  XX
#      20        128 LSRK  230987.639  812.5       104000      230987.639  XX
#      21        128 LSRK  231069.641  812.5       104000      231069.641  XX
#      22        128 LSRK  231144.831  812.5       104000      231144.831  XX
#      23        128 LSRK  231226.833  812.5       104000      231226.833  XX


#### Brief description of the purposes of the sources in this project; 
#### () gives source id number. This is a six-pointing mosaic 
#### 
#### !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
#### Bandpass calibrator: 3c454.3 (9)
#### Gain calibrators: 1743-038 (1), 1908-201 (6)
#### Flux calibrator: 3c454.3 (9) and uranus (13)
#### Science target g19: (2,3,4,5,7,8)
#### other sources were not part of this program
#### !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

#### Initialize scratch columns 

clearcal(vis='g19_d2usb.ms')

#######################################################################
#### Plot sources as function of time
#######################################################################

plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",datacolumn="data",
       iteration="baseline",
       antenna="",spw="",field="1~9,13",
       averagemode="vector",width="allspw",timebin="0",
       crossscans=T,crossbls=False,stackspw=F,
       extendflag=F,extendchan="",extendspw="",
       extendant="",extendtime="",subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")

#### At the moment, this plot takes a long time because the scratch
#### columns must be initialized.
# Pause script if you are running scriptmode
# NOTE: if you are running interactively and cutting/pasting you do not
# need to include these lines that come after plotting
user_check=raw_input('Return to continue script\n')

#######################################################################
#### BANDPASS CALIBRATION
#######################################################################

#### The following can be used to plot the antenna positions, note
#### that antenna 8 was in the barn for repair.

default('plotxy')
plotxy(vis="g19_d2usb.ms",xaxis="x")

#### Make an inital amplitude vs. channel plot of the bandpass
#### calibrator 3c454.3 averaging all the data together in
#### time. crossscan=T needed because there are two scans on
#### 3c454.3. Note the rolloff of 24 individual spw. Try zooming to
#### see better.

plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp",datacolumn="data",
       iteration="baseline",
       antenna="",spw="",field="9",
       averagemode="vector",width="1",timebin="all",
       crossscans=T,crossbls=False,stackspw=F,
       extendflag=F,extendchan="",extendspw="",
       extendant="",extendtime="",subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

#### Next plot phase vs. channel. Note the large descrepancy in phase
####  for spw=4-7 on some baselines, this will be important later.

plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="phase",datacolumn="data",
       iteration="baseline",
       antenna="",spw="",field="9",
       averagemode="vector",width="1",timebin="all",
       crossscans=T,crossbls=False,stackspw=F,
       extendflag=F,extendchan="",extendspw="",
       extendant="",extendtime="",subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

####  SMA spw overlap by a fair amount, flag first 8 channels on either 
####  side of each spw to aid coherence of inital phase calibration 
####  below. 

default('flagdata')
flagdata(vis='g19_d2usb.ms', mode='manualflag',spw='0~23:0~7;120~127')

####  May wish to repeat plotxy above to check outcome of flagdata

####  Now plot amplitude and phase as a function of time.

plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",datacolumn="data",
       iteration="baseline", antenna="",spw="0~23",field="9",
       averagemode="vector",width="allspw",timebin="0",
       subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="phase",datacolumn="data",
       iteration="baseline", antenna="",spw="0~23",field="9",
       averagemode="vector",width="allspw",timebin="0",
       subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

#### The gap in 3C454.3 observation is due to a pause to determine
#### pointing solutions.  The flags below catch times just before and
#### after when the amplitudes and phase are off. This flagging could
#### also have been accomplished interactively.

default('flagdata')
flagdata(vis="g19_d2usb.ms",field="9",timerange='12:16:58~12:19:00')
flagdata(vis="g19_d2usb.ms",field="9",timerange='12:08:00~12:09:30')

####  Start with an initial phase only solution for 3c454.3 to take out 
####  major phase variations with time. For maximum sensitivity we 
####  average all the spw together with the exception of spw 4-7 which
####  show large phase offsets from the other spw in phase vs. channel.  

gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.pcal',
	field='9', spw='0~3,8~23', gaintype='G', minsnr=2.0,
	refant='ANT3', calmode='p',solint='int', combine='spw')

plotcal(caltable="g19_d2usb.ms.pcal",xaxis="time",yaxis="phase",
	field="9",antenna="",spw="",timerange="",subplot=331,
	overplot=False,clearpanel="Auto",iteration="antennas",
	plotsymbol="o",plotcolor="blue",showgui=True,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

####  A regular 'B' type solution can work well if you have lots of
####  S/N on your bandpass calibrator. You may wish to make this table
####  (uncomment commands below) and compare to 'BPOLY' solution
####  further down.

# bandpass(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.bpcal',
# 	field='9', spw='0~23', bandtype='B', solint='inf', combine='scan', 
# 	refant='ANT3', gaincurve=False, opacity=0.0,
# 	gaintable='g19_d2usb.ms.pcal',spwmap=[0]) 
#
#plotcal(caltable="g19_d2usb.ms.bpcal",xaxis="freq",yaxis="amp",
#	field="9",antenna="",spw="",timerange="",subplot=311,
#	overplot=False,clearpanel="Auto",iteration="antenna",
#	plotsymbol="o",plotcolor="blue",showgui=True,figfile="") 

#### A polynomial solution is more sensitive than simply fitting a 
#### polynomial to the channel based solutions 'B' after the fact. For  
#### the mm/submm where the S/N is often low on the bandpass calibrator 
#### 'BPOLY' is often the best choice, but every dataset is different. 
#### You may especially need to play with degamp and degphase. The 
#### combine='scan' is needed to average together the two 3c454.3 scans.
#### The spwmap=[0] is telling it to apply the phase solution that was 
#### found from the average of all spw to each spw before determining 
#### the bandpass solution. 

bandpass(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.bpoly7',
	field='9', spw='0~23', bandtype='BPOLY', 
	solint='inf', combine='scan', 
	refant='ANT3', degamp=7, degphase=7,
	gaintable='g19_d2usb.ms.pcal',spwmap=[0])


plotcal(caltable="g19_d2usb.ms.bpoly7",xaxis="freq",yaxis="amp",
	field="9",antenna="",spw="",timerange="",subplot=311,
	overplot=False,clearpanel="Auto",iteration="antenna",
	plotsymbol="o",plotcolor="blue",showgui=True,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

#### Unfortunately, this kind of plot with many spw and colors 
#### takes a long time. Be patient for now its being worked on. You 
#### can safely ignore where the solution gets large at edges of 
#### spw -- its showing where we flagged channels -- this is 
#### being worked on. Zoom in if necessary to see that inner 
#### portions are nice and flat.

#### The next step is to apply the bandpass solution while
#### re-calculating the phase only solution for 3c454.3 on a short
#### (integration time) time interval, but now including all spw
#### (including 4~7 since antenna based solutions will have taken out
#### phase vs.spw differences).

gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.pcal2',
	field='9', spw='0~23', gaintype='G', minsnr=2.0,
	refant='ANT3', calmode='p',solint='int', combine='spw',
	gaintable=['g19_d2usb.ms.bpoly7'])


plotcal(caltable="g19_d2usb.ms.pcal2",xaxis="time",yaxis="phase",
	field="9",antenna="",spw="",timerange="",subplot=331,
	overplot=False,clearpanel="Auto",iteration="antennas",
	plotsymbol="o",plotcolor="blue",showgui=True,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

#### Next we apply the antenna based bandpass, and phase only solution
#### to solve for amplitude solutions over a longer solint (to
#### increase S/N).The spwmap=[[0],[]] tells it that the pcal solution
#### is independent of spw while the bpoly7 solution is per spw,
#### spwmap order must be in the same order as the tables in
#### gaintable.

gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.apcal',
	field='9', spw='0~23', gaintype='G', minsnr=2.0,
	refant='ANT3', calmode='a',solint='300s', combine='spw',
	gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.bpoly7'],
	spwmap=[[0],[]])

plotcal(caltable="g19_d2usb.ms.apcal",xaxis="time",yaxis="amp",
	field="9",antenna="",spw="",timerange="",subplot=331,
	overplot=False,clearpanel="Auto",iteration="antennas",
	plotsymbol="o",plotcolor="blue",showgui=True,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

#### Although in previous steps solutions were applied "on-the-fly",
#### the data were not actually changed. Running applycal will fill
#### the corrected data column in the ms. Again be careful with
#### specification of gaintable, and keeping the same order in spwmap
#### and gainfield. NOTE this applycal is not necessary for later
#### steps -- it is being done in order to allow plotting of data
#### calibration to this point.

applycal(vis='g19_d2usb.ms',spw='0~23', field='9',
	gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.apcal',
	'g19_d2usb.ms.bpoly7'],
	spwmap=[[0],[0],[]],gainfield=['9','9','9'])


plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp",
       datacolumn="corrected",iteration="baseline",
       antenna='',spw="0~23",field="9",
       averagemode="vector",width="1",timebin="all",crossscans=T,
       crossbls=False,stackspw=False,extendflag=F,
       extendchan="",extendspw="",extendant="",extendtime="",
       plotcolor='darkcyan',subplot=331,
       multicolor="none",figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

#### As this plotxy shows, after antenna based bandpass, SMA data
#### often show residual baseline based problems as a function of spw
#### (i.e. a few abnormally low amplitudes for a few spw on a few
#### baselines). To fix this problem we need to do a baseline based
#### amplitude only solution per spw, but we'll need to remove an
#### overall normalization factor first by combining all spw (this
#### option will be incorporated into blcal one day so that only one
#### run is necessary).

blcal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.blcal1',
	field='9', spw='0~23',solint='inf', combine='spw,scan', 
	gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.apcal',
	'g19_d2usb.ms.bpoly7'], 
	calmode='a',spwmap=[[0],[0],[]])


blcal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.blcal2',
	field='9', spw='0~23',solint='inf', combine='scan', 
	gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.apcal',
	'g19_d2usb.ms.blcal1','g19_d2usb.ms.bpoly7',], 
	calmode='a',spwmap=[[0],[0],[0],[]])

#### Now we can do a new applycal including ONLY the second blcal2
#### solution. Note that applycal overwrites the corrected column so
#### there is no need to clearcal. Like the antenna based bandpass 
#### the blcal solution is independent for each spw, so its spwmap 
#### entry is [].

applycal(vis='g19_d2usb.ms',spw='0~23', field='9',
	gaintable=['g19_d2usb.ms.pcal2','g19_d2usb.ms.apcal',
	'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],
	spwmap=[[0],[0],[],[]],gainfield=['9','9','9','9'])

plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp",
       datacolumn="corrected",iteration="baseline",
       antenna='',spw="0~23",field="9",
       averagemode="vector",width="1",timebin="all",crossscans=T,
       crossbls=False,stackspw=False,extendflag=F,
       extendchan="",extendspw="",extendant="",extendtime="",
       plotcolor='darkcyan',subplot=331,
       multicolor="none",figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

#### Note how flat the full bandpass is after all this calibration

#######################################################################
#### INSPECTION OF GAIN CALIBRATORS AND SET FLUX SCALE
#######################################################################

#### Apply antenna and baseline based bandpass solutions to the 
#### gain calibrators

applycal(vis='g19_d2usb.ms',spw='0~23', field='1,6',
	gaintable=['g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],
	spwmap=[[],[]],gainfield=['9','9'])


#### Have a look at both calibrators with channel vs. amp.

plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp",datacolumn="corrected",
       iteration="baseline", antenna="",spw="",field="1",
       averagemode="vector",width="1",timebin="all",
       subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

plotxy(vis="g19_d2usb.ms",xaxis="channel",yaxis="amp",datacolumn="corrected",
       iteration="baseline", antenna="",spw="",field="6",
       averagemode="vector",width="1",timebin="all",
       subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

#### Have a look at both calibrators in phase and amplitude
#### vs. time. For the moment, unfortunately multicolor='field'
#### doesn't work, so easiest to do seperately.

plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",datacolumn="corrected",
       iteration="baseline", antenna="",spw="0~23",field="1",
       averagemode="vector",width="allspw",timebin="0",
       subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="amp",datacolumn="corrected",
       iteration="baseline", antenna="",spw="0~23",field="6",
       averagemode="vector",width="allspw",timebin="0",
       subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="red",
       overplot=F,showflags=False,interactive=True,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

#### On many baselines the first integration is low. These can be flagged 
#### with flagdata. The integration time of these data is 30.9s. More recent 
#### SMA data rarely shows this problem.

flagdata(vis="g19_d2usb.ms",mode="quack",field='',spw="0~23",quackinterval=40.0)

plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="phase",datacolumn="corrected",
       iteration="baseline", antenna="",spw="0~23",field="1",
       averagemode="vector",width="allspw",timebin="0",
       subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

plotxy(vis="g19_d2usb.ms",xaxis="time",yaxis="phase",datacolumn="corrected",
       iteration="baseline", antenna="",spw="0~23",field="6",
       averagemode="vector",width="allspw",timebin="0",
       subplot=331,plotsymbol=".",
       multicolor="none",plotcolor="darkcyan",
       overplot=False,showflags=False,interactive=True,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

#### At the time of these observations, the calibrator J1911 had the 
#### most recent and relatively stable monitor data (SMA website), so 
#### we will use it for absolute flux calibration. In future tools for 
#### using resolved planets will be implemented.
#### 1911-201
#### 1mm   17 Jul 2005 11:38   SMA      274.5    1.88 +/-  0.13    mgurwell  
#### 1mm   05 Sep 2005 08:10   SMA      221.4    1.99 +/-  0.12    mgurwell  

setjy(vis='g19_d2usb.ms', field='6', spw='0~23', fluxdensity=[2.0,0.,0.,0.])

#######################################################################
#### FINAL GAIN AND FLUX CALIBRATION FOR ALL CALIBRATORS
#######################################################################

####  To avoid decorrelation of the amplitude, the phase solutions fed
####  to the amplitude calibration step should be at the shortest
####  solint your S/N can support, typically the integration time.

gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.allpcal',
	field='1,6,9', spw='0~23', gaintype='G', minsnr=2.0,
	refant='ANT3', calmode='p',solint='int', combine='spw',
	gaintable=['g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],
	spwmap=[[],[]])

#### Unless you chose to do some more sophisticated smoothing, the
#### phase solutions applied to your target might as well be on the
#### scan time of your phase calibrator(s). Without combine='scan',
#### the solution will not cross scan boundaries even though
#### solint='inf'. This is an easy way to "manually" interpolate.

gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.allpcalscan',
	field='1,6,9', spw='0~23', gaintype='G', minsnr=2.0,
	refant='ANT3', calmode='p',solint='inf', combine='spw',
	gaintable=['g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],
	spwmap=[[],[]])

#### Plot the phase solutions of the phase calibrators. Can repeat for
#### both short and long solution intervals.

plotcal(caltable="g19_d2usb.ms.allpcalscan",xaxis="time",yaxis="phase",
	field="1,6",antenna="",spw="",timerange="",subplot=331,
	overplot=False,clearpanel="Auto",iteration="antennas",
	plotsymbol="o",plotcolor="blue",showgui=True,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

#### Amplitude calibration with longer solint and applying 
#### short solint phase solution.

gaincal(vis='g19_d2usb.ms', caltable='g19_d2usb.ms.allapcal',
	field='1,6,9', spw='0~23', gaintype='G', minsnr=2.0,
	refant='ANT3', calmode='ap',solint='300s', combine='spw',
	gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.bpoly7',
	'g19_d2usb.ms.blcal2'],spwmap=[[0],[],[]])

#### Plot the amplitude solutions

plotcal(caltable="g19_d2usb.ms.allapcal",xaxis="time",yaxis="amp",
	field="1,6",antenna="",spw="",timerange="",subplot=331,
	overplot=False,clearpanel="Auto",iteration="antennas",
	plotsymbol="o",plotcolor="blue",showgui=True,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

#### Derive absolute flux calibration (based on setjy above)

fluxscale(vis='g19_d2usb.ms',caltable='g19_d2usb.ms.allapcal',
	fluxtable='g19_d2usb.ms.fluxcal',reference='6')

#######################################################################
#### APPLY FINAL CALIBRATIONS
#######################################################################

#### Final application of calibration for each source is done
#### independently below. Flux densities from SMA flux monitoring are
#### provided as checks of the absolute flux calibration.

applycal(vis='g19_d2usb.ms',spw='0~23', field='9',
	gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.fluxcal',
	'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]],
	gainfield=['9','9','9','9'])

plotxy(vis="g19_d2usb.ms",xaxis="uvdist",yaxis="amp",datacolumn="corrected",
       iteration="", antenna="",spw="0~23",field="9",
       averagemode="vector",width="allspw",timebin="0",
       subplot=111,plotsymbol=".",interactive=T,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

####  3C454.3
####  1mm   19 Aug 2005 16:01   SMA      225.3   27.36 +/-  1.37    mgurwell 

applycal(vis='g19_d2usb.ms',spw='0~23', field='1',
	gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.fluxcal',
	'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]],
	gainfield=['1','1','9','9'])

plotxy(vis="g19_d2usb.ms",xaxis="uvdist",yaxis="amp",datacolumn="corrected",
       iteration="", antenna="",spw="0~23",field="1",
       averagemode="vector",width="allspw",timebin="0",
       subplot=111,plotsymbol=".",interactive=T,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')


#### J1743-038
#### 1mm   26 Aug 2005 06:05   SMA      225.6    1.71 +/-  0.09    mgurwell 


applycal(vis='g19_d2usb.ms',spw='0~23', field='6',
	gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.fluxcal',
	'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]],
	gainfield=['6','6','9','9'])

plotxy(vis="g19_d2usb.ms",xaxis="uvdist",yaxis="amp",datacolumn="corrected",
       iteration="", antenna="",spw="0~23",field="6",
       averagemode="vector",width="allspw",timebin="0",
       subplot=111,plotsymbol=".",interactive=T,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

#### J1911-201
#### 1mm   17 Jul 2005 11:38   SMA      274.5    1.88 +/-  0.13    mgurwell  
#### 1mm   05 Sep 2005 08:10   SMA      221.4    1.99 +/-  0.12    mgurwell  

#### Since Uranus is resolved we apply the solutions from 3C454.3 to
#### achieve at lest some reduction in instrumental gain variations.

applycal(vis='g19_d2usb.ms',spw='0~23', field='13',
	gaintable=['g19_d2usb.ms.allpcal','g19_d2usb.ms.fluxcal',
	'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]],
	gainfield=['9','9','9','9'])

plotxy(vis="g19_d2usb.ms",xaxis="uvdist",yaxis="amp",datacolumn="corrected",
       iteration="", antenna="",spw="0~23",field="13",
       averagemode="vector",width="allspw",timebin="0",
       subplot=111,plotsymbol=".",figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

#### A UV-distance plot of the expected flux of Uranus can be obtained
#### from the SMA webpage. In future we will implement tools to use 
#### planets for absolute flux calibration.

applycal(vis='g19_d2usb.ms',spw='0~23', field='2,3,4,5,7,8',
	gaintable=['g19_d2usb.ms.allpcalscan','g19_d2usb.ms.fluxcal',
	'g19_d2usb.ms.bpoly7','g19_d2usb.ms.blcal2'],spwmap=[[0],[0],[],[]],
	gainfield=['1,6','1,6','9','9'])

#### Field 2 (pointing a) has the strongest CO emission near the phase 
#### center, and the 12CO(2-1) line is clearly apparent in the plot below. 
#### Future enhancement will include ability to offset the phase center 
#### to a different position for plotting.

plotxy(vis="g19_d2usb.ms",xaxis="freq",yaxis="amp",
	datacolumn="corrected",iteration="",selectdata=True,
	antenna='',spw="0~23",field="2",timerange="",averagemode="vector", 
	width="1",timebin="all",crossscans=T,crossbls=T,
	stackspw=False,subplot=111,interactive=T,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

#### Note if you want to display velocity you will need to also set the
#### restfreq to the line of interest.

###################################################################
#### SPLIT OFF TARGET DATA FOR FURTHER INSPECTION
###################################################################

split(vis="g19_d2usb.ms",outputvis='g19_d2usb_targets.ms',
      field="2,3,4,5,7,8")

listobs(vis="g19_d2usb_targets.ms",verbose=F)

#### ID   Code Name          Right Ascension  Declination   Epoch
#### 0         g19.3a        18:25:58.70      -12.03.57.80  J2000
#### 1         g19.3b        18:25:57.40      -12.04.18.50  J2000
#### 2         g19.3c        18:25:56.09      -12.04.28.20  J2000
#### 3         g19.3d        18:25:54.60      -12.04.23.10  J2000
#### 4         g19.3e        18:25:54.80      -12.04.43.80  J2000
#### 5         g19.3f        18:25:53.30      -12.04.51.00  J2000

#### Note that after splitting, the "corrected" data column in the
#### original vis file will become the "data" and "corrected" columns
#### and the field IDs are renumbered.

plotxy(vis="g19_d2usb_targets.ms",xaxis="uvdist",yaxis="amp",datacolumn="data",
       iteration="", antenna="",spw="0~23",field="",
       averagemode="vector",width="allspw",timebin="0",
       subplot=111,plotsymbol=".",interactive=T,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

#### This plotxy will take more time than usual, because scratch 
#### columns will be made in the new ms.

plotxy(vis="g19_d2usb_targets.ms",xaxis="time",yaxis="amp",datacolumn="data",
       iteration="", antenna="",spw="0~23",field="",
       averagemode="vector",width="allspw",timebin="0",
       subplot=111,plotsymbol=".",interactive=T,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

#### Data at times > 11:04:00 show lots of scatter: the elevation of 
#### these data are very low and are worth flagging. One can plotxy 
#### with xaxis="time",yaxis="elevation" to see this.

flagdata(vis="g19_d2usb_targets.ms",field="",timerange='>11:04:00')

#### Note that because the Doppler tracked frequency was 13CO (+ 10GHz
#### in the USB), the velocity diplayed in in plotxy will not be
#### correct unless you set the restfreq to the line of interest.

#### CO(2-1) restfreq=230.53797 GHz

plotxy(vis="g19_d2usb_targets_line.ms",xaxis="channel",yaxis="amp",
	datacolumn="data",iteration="",selectdata=True,
	antenna='',spw="0~23",field="0",timerange="",averagemode="vector",
        restfreq='230.53797GHz',
	width="1",timebin="all",crossscans=T,crossbls=T,
	stackspw=False,subplot=111,interactive=T,figfile="")

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

#### The other weak line detected in this sideband is not visible in 
#### a vector average plot.
#### A-CH3OH (8 -1 8 0 -- 7 +0 7 0) 229.75876 GHz   

###################################################################
#### CONTINUUM IMAGING AND SUBSTRACTION
###################################################################
###################################################################

#### In the current release 2.3.1, you cannot use the uvcontsub task
#### to remove the continuum UNLESS you have sufficient line free
#### channels in the spw containing your line(s) of interest,
#### i.e. continuum isn't removed from spw with no line-free
#### channels. An alternative method is described for the 12CO case
#### which spans 3 spw.

#### copy the split dataset to a new name to preserve a "continuum" dataset
os.system('rm -rf g19_d2usb_targets_line.ms')

os.system('cp -r g19_d2usb_targets.ms g19_d2usb_targets_line.ms')

#### after steps below, g19_d2usb_targets_line.ms will have the continuum
#### subtracted.  If you want to re-image the continuum after these steps
#### use g19_d2usb_targets.ms

#### Set visibility file to use for all steps below. This is just to
#### demonstrate this alternative way of setting things, one could
#### also set it explicitly in each task. If you are cutting and pasting 
#### be careful to have set this in the current session, or re-set it if you 
#### change to a different vis in the meantime.

vis='g19_d2usb_targets_line.ms'

#### refresh model_data and corrected_data
clearcal(vis);

#### Image the continuum (interactively) -- it is important to only
#### include clean components you are certain represent real
#### emission. This means that you should make an interactive clean
#### mask. You will need to stop clean (interactively) when you are happy
#### with the residuals. 

#### Note that for sparse mosaics like this one, the subparameter
#### ftmachine='ft' will typically give better results than the
#### default ftmachine='mosaic'. Also note that in a mosaic, the rms
#### noise changes as a function of position. Therefore, the rms noise
#### that you calculate from the radiometer equation is only valid
#### where the combined primary beam response of the mosaic is near
#### one. Clean weights the decision of what to clean next by the
#### primary beam response function (PBRF), so emission near PBRF=1
#### will be cleaned faster than elsewhere. Thus, if you have emission
#### away from the "sweet spot(s)" you will need to clean to a
#### threshold lower than you might expect from the anticipated rms
#### noise. Development is on-going to provide other options which may 
#### work better for the case of strong emission near image boundaries.

contimag='g19_d2usb_cont'
clean(vis=vis,imagename=contimag,
      field='',spw='0~12,16~23',
      mode='mfs',
      niter=200,gain=0.1,threshold=0.0,
      psfmode='clark',imagermode='mosaic',scaletype='SAULT',
      ftmachine='ft',
      interactive=T,imsize=500,cell="0.5arcsec",
      phasecenter="J2000 18h25m56.09 -12d04m28.20",
      pbcor=T,minpb=0.15)

#### Here you can start the viewer by typing viewer and selecting this
#### image from file menu or start with

viewer(contimag'+.image')

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

#### establish the continuum model_data for all spws via mosaic clean
#### (niter=0) (can't use the task ft because it doesn't know about
#### mosaic primary beam).

imag='g19_d2usb_ft'  # (throwaway)
clean(vis=vis,imagename=imag,
      field='',spw='',
      mode='mfs',
      niter=0,gain=0.1,threshold=0.0,
      psfmode='clark',imagermode='mosaic',scaletype='SAULT',
      ftmachine='ft',
      interactive=F,imsize=500,cell="0.5arcsec",
      phasecenter="J2000 18h25m56.09 -12d04m28.20",
      pbcor=T,minpb=0.15,
      modelimage=contimag+'.model',
      cyclefactor=0.25,cyclespeedup = 50)


# subtract model from corrected to leave only lines for imaging
uvsub(vis=vis);

# verify that the continuum is gone...
contsubimag='g19_d2usb_cont_gone'
clean(vis=vis,imagename=contsubimag,
      field='',spw='0~12,16~23',
      mode='mfs',
      niter=0,gain=0.1,threshold=0.0,
      psfmode='clark',imagermode='mosaic',scaletype='SAULT',
      ftmachine='ft',
      interactive=F,imsize=500,cell="0.5arcsec",
      phasecenter="J2000 18h25m56.09 -12d04m28.20",
      pbcor=T,minpb=0.15,
      cyclefactor=0.25,cyclespeedup = 50)

viewer(contsubimag+'.image')

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')


###################################################################
#### Spectral line imaging
###################################################################

#### This image selects a channel range suitable to image the
#### 12CO(2-1) line.  Notice use of restfreq parameter to get velocity
#### scale correct.

vis='g19_d2usb_targets_line.ms'
lineimag='g19_d2usb_coline'
clean(vis=vis,imagename=lineimag,
      field='',spw='',
      mode='channel',start=730,nchan=60,width=2,
      niter=10000,gain=0.1,threshold=0.2,
      psfmode='clark',imagermode='mosaic',scaletype='SAULT',
      ftmachine='ft',
      restfreq='230.53797GHz',
      interactive=T,imsize=500,cell="0.5arcsec",
      phasecenter="J2000 18h25m56.09 -12d04m28.20",
      pbcor=T,minpb=0.15)

#### It is essential to run with interactive=T and make a clean mask.

viewer('g19_d2usb_coline.image')

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

#### The variations in apparent minpb as you go from channel to
#### channel are due to the channel weights being incorporated into
#### the mosaic primary beam response function. i.e. the channel
#### weights are better where the spw overlap, causing a change in the 
#### minpb. This behavior is being changed in a future release so a 
#### uniform cutoff can be applied.

#### We can fix this problem using the following immath to
#### mask the cube based on a channel with a reasonable mask

immath(outfile='g19_d2usb_coline.fluxch_25',mode='evalexpr',
	expr='"g19_d2usb_coline.flux"',chans='24')

#### Now make 0th and 1st moment maps.

immoments(imagename='g19_d2usb_coline.image',moments=[0],axis=3,
	  planes='15~48',outfile='g19_d2usb_coline.mom0')

immath(outfile='g19_d2usb_coline.mom0masked',mode='evalexpr',
       expr='"g19_d2usb_coline.mom0"["g19_d2usb_coline.fluxch_25">0.2]')

immoments(imagename='g19_d2usb_coline.image',moments=[1],axis=3,
	  planes='15~48',outfile='g19_d2usb_coline.mom1',
	  excludepix=[-100,2])

immath(outfile='g19_d2usb_coline.mom1masked',mode='evalexpr',
       expr='"g19_d2usb_coline.mom1"["g19_d2usb_coline.mom0masked">15]') 

viewer('g19_d2usb_coline.mom1masked')

# Pause script if you are running in scriptmode
user_check=raw_input('Return to continue script\n')

#### By clicking on "folder" icon in viewer window you can bring up
#### the moment0 as a contour image. The strange appearence of the the
#### mom1 map toward the NE source is probably due to more than one
#### outflow being superposed.