###
### Script for PdBI project L02D data:
### A Keplerian disk around the protostar IRAS 20126+4104
### Single field, 1 & 3mm continuum DSB
### many lines in each sub-band with strong continuum emission.
### Day 2 - 26 Jan 2002 data
###
### Steven T. Myers (NRAO) Version 14-Jan-2004
###

###########################################################
## Tested version stable_ss4 v1.9 Build 355
## Expected run time 7 minutes

###########################################################
## Include relevant tools if not already done:
##
include 'almati2ms.g';
include 'ms.g';             
include 'msplot.g';
include 'flagger.g'
include 'imager.g';
include 'calibrater.g';
include 'iramcalibrater.g';
include 'viewer.g';
include 'image.g';
include 'imagepol.g';      


###########################################################
## Fill Data: (unless already in ms)
## New in SS4: do PHCOR step to swap out bad baselines
almatifiller(msfile='26-jan-2002-l02d.ms',
             fitsdir='/home/cluster4/jmcmulli/almatest/L102D/data',
             pattern='26-jan-2002-l02d-*',
             append=F, compress=F, obsmode='CORR', chanzero='TIME_AVG',
             dophcor=T);

###########################################################
## Obtain a detailed summary (with default catalog tool):
##
dc.summary('26-jan-2002-l02d.ms');

###########################################################
## Flag shadowed data:
##
shadow(msname="26-jan-2002-l02d.ms", trial=F, minsep=15.0);

###########################################################
## Manually plot the data for visual inspection,
##  interactively flag bad data:
##
## see the data in x-y plot format: 
#plot := msplot(msfile='26-jan-2002-l02d.ms', edit=T);
# Allow editing as well as inspection. See cookbook for editing steps.

###########################################################
## Flag bad and irrelevant data on first day 
## - the real data for this file is on 28-Jan
## Flag spectral windows of same shape together 
## - note flagger cannot handle different shaped windows
## - could use autoflag instead
##
fg:=flagger('26-jan-2002-l02d.ms');
fg.settimerange(starttime='27-Jan-2002/05:00:00.0',
        endtime='27-Jan-2002/07:00:00.0');
# First the 1-channel windows:
fg.setids(spectralwindowid=[1,2,5,6,9,10,13,14,17,18,21,22,25,26,29,30]);
fg.flag()
#Should see: Flagging 2400 rows
# Next the 128-channel windows:
fg.setids(spectralwindowid=[7,8,15,16,19,20,27,28,31,32]);
fg.flag()
#Should see: Flagging 1500 rows
# Next the 512-channel windows:
fg.setids(spectralwindowid=[3,4,11,12,23,24]);
fg.flag() 
#Should see: Flagging 900 rows
# Now the Gibbs channels for 128-ch 1.25MHz wins
fg.reset();
fg.setids(spectralwindowid=[7,8,19,20]);
fg.setchan(chan=[1,2,127,128]);
fg.flag();
#Should see: Flagging 4440 rows
# Now the Gibbs channels for 128-ch 2.5MHz wins
fg.reset();
fg.setids(spectralwindowid=[15,16,27,28,31,32]);
fg.setchan(chan=[1,2,64,65,127,128]);
fg.flag();
#Should see: Flagging 6500 rows
# Now the Gibbs channels for 512-ch wins
fg.reset();
fg.setids(spectralwindowid=[3,4,11,12,23,24]);
fg.setchan(chan=[1,2,3,4]);
fg.flag();
#Should see: Flagging 6580 rows

fg.done();

###########################################################
##  Set flux density scale for the flux density calibrators
##
imgr:=imager(filename='26-jan-2002-l02d.ms');
# Set for the 3mm windows
for(spw in 1:20){
     imgr.setjy(fieldid=1,spwid=spw,fluxdensity=[10.87,0,0,0]);
     }

for(spw in 1:20){
     imgr.setjy(fieldid=3,spwid=spw,fluxdensity=[2.79,0,0,0]);
     }
# Now set for the 1mm windows
for(spw in 21:32){
     imgr.setjy(fieldid=1,spwid=spw,fluxdensity=[6.00,0,0,0]);
     }

for(spw in 21:32){
     imgr.setjy(fieldid=2,spwid=spw,fluxdensity=[1.75,0,0,0]);
     }

for(spw in 21:32){
     imgr.setjy(fieldid=3,spwid=spw,fluxdensity=[1.80,0,0,0]);
     }

###########################################################
##########      3 mm Calibration                 ##########
###########################################################

###########################################################
## Get first cut phase solutions or BPass determination:
## 
cal:=calibrater(filename='26-jan-2002-l02d.ms');
## Do 128-channel windows [7,8,15,16,19,20]
## Leave out first and last 2 unflagged channels
## First the LSB [7,15,19]
cal.setdata(mode='channel', start=4, step=1, nchan=122, 
     msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID IN [7,15,19]');
cal.setsolvegainspline(table='26-jan-2002-l02d.3mmLSB.gcal0',
     mode='PHAS',                 
     preavg=0,                         #  no preavg too short,
     splinetime=10000,                 #  spline fit timescale (~3hrs),
     refant=2,                         #  refant = 2,
     npointaver=10,                    #  minimize phase wrap pblms
     phasewrap=260,
     append=F);                        #  create new table
cal.solve();
##
## Plot files will be:
## 26-jan-2002-l02d.3mmLSB.gcal0.PHASE.ps
## 26-jan-2002-l02d.3mmLSB.gcal0.PHASE.log
##
## Now the USB [8,16,20]
cal.reset();
cal.setdata(mode='channel', start=4, step=1, nchan=122, 
     msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID IN [8,16,20]');
cal.setsolvegainspline(table='26-jan-2002-l02d.3mmUSB.gcal0',
     mode='PHAS',                 
     preavg=0,                         #  no preavg too short,
     splinetime=10000,                 #  spline fit timescale (~3hrs),
     refant=2,                         #  refant = 2,
     npointaver=10,                    #  minimize phase wrap pblms
     phasewrap=260,
     append=F);                        #  create new table
cal.solve();
##
## Plot files will be:
## 26-jan-2002-l02d.3mmUSB.gcal0.PHASE.ps
## 26-jan-2002-l02d.3mmUSB.gcal0.PHASE.log
##
###########################################################
## Derive BPass calibration:
## 
## First for 3mm LSB continuum phase solns
cal.reset();                                  
cal.setdata(msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID IN [7,15,19]'); 
cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mmUSB.gcal0');
# Note we are applying the USB phases to the LSB to preserve LSB/USB offsets     
cal.setsolvebandpoly(table='26-jan-2002-l02d.3mm.bpoly', 
    degamp=12, degphase=24, visnorm=F, bpnorm=F, refant=2, append=F);  
cal.solve();                                  
##
## Plot files will be:
## 26-jan-2002-l02d.3mm.bpoly.ps
## Need to move these
shell('mv 26-jan-2002-l02d.3mm.bpoly.ps 26-jan-2002-l02d.3mmLSB.bpoly.ps')
##
## Now for 3mm USB  continuum phase solns...
cal.reset();                                  
cal.setdata(msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID IN [8,16,20]'); 
cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mmUSB.gcal0');
## NOTE append=T in caltable 
cal.setsolvebandpoly(table='26-jan-2002-l02d.3mm.bpoly', 
    degamp=18, degphase=24, visnorm=F, bpnorm=F, refant=2, append=T);  
cal.solve();                                  
##
## Plot files will be:
## 26-jan-2002-l02d.3mm.bpoly.ps
## Need to move these
shell('mv 26-jan-2002-l02d.3mm.bpoly.ps 26-jan-2002-l02d.3mmUSB.bpoly.ps')
##
###########################################################
## Phase calibrate 3mm LSB & USB 128-ch data:
## all cals
##
cal.reset(); 
cal.setdata(mode='channel', start=4, step=1, nchan=122, 
    msselect='FIELD_ID IN [1,2,3] && SPECTRAL_WINDOW_ID IN [7,8,15,16,19,20]');
cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mm.bpoly'); 
cal.setsolvegainspline(table='26-jan-2002-l02d.3mm.gcal',
    mode='PHAS', preavg=0, splinetime=10000,
    refant=2, npointaver=10, phasewrap=260, append=F); 
cal.solve(); 
##
## Plot files will be:
## 26-jan-2002-l02d.3mm.gcal.PHASE.ps
## 26-jan-2002-l02d.3mm.gcal.PHASE.log
##
###########################################################
## Do a residual T solution to allow fluxscale
## 
cal.reset();                                  # Reset calibrator.
cal.setdata(mode='channel', start=4, step=1, nchan=122, 
    msselect='FIELD_ID IN [1,2,3] && SPECTRAL_WINDOW_ID IN [7,8,15,16,19,20]'); 
cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mm.bpoly');
cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal');  
cal.setsolve(type='T',                        # Solve for any residual phase 
    table='26-jan-2002-l02d.3mm.tcal',        # errors with a timescale of 
    t=300,                                    # 5 min.  
    refant=2);                                # Use ant 2 as reference.
cal.solve();                                  # Solve for solutions. 

###########################################################
## Get the fluxscale for MWC349
## 
cal.fluxscale(tablein='26-jan-2002-l02d.3mm.tcal',
    tableout='26-jan-2002-l02d.3mm.fcal',
    reference="3C273 2013+370", transfer='MWC349');
##
## Should get:
## Flux density for MWC349 (spw=7) is:     1.048 +/-  0.006 Jy
## Flux density for MWC349 (spw=8) is:     1.080 +/-  0.007 Jy
## Flux density for MWC349 (spw=15) is:     1.057 +/-  0.008 Jy
## Flux density for MWC349 (spw=16) is:     1.091 +/-  0.009 Jy
## Flux density for MWC349 (spw=19) is:     1.028 +/-  0.008 Jy
## Flux density for MWC349 (spw=20) is:     1.094 +/-  0.004 Jy
## 
###########################################################
## Apply average fluxscale for MWC349
## 
for(spw in 1:20){
       imgr.setjy(fieldid=2,spwid=spw,fluxdensity=[1.06,0,0,0]);
       }

###########################################################
## Amplitude calibration of 3mm 128-ch data:
## all cals
## note that the phase solution in the table will be combined with the
## amplitude solution...doing effectively an incremental combination
##
cal.reset(); 
cal.setdata(mode='channel', start=4, step=1, nchan=122, 
    msselect='FIELD_ID IN [1,2,3] && SPECTRAL_WINDOW_ID IN [7,8,15,16,19,20]');
cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mm.bpoly');
cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal');  
cal.setsolvegainspline(table='26-jan-2002-l02d.3mm.gcal',
    mode='AMP', preavg=0.0, splinetime=10000, 
    refant=2, npointaver=10, phasewrap=260);
cal.solve();  
##
## Plot files will be:
## 26-jan-2002-l02d.3mm.gcal.AMP.ps
## 26-jan-2002-l02d.3mm.gcal.AMP.log
##
###########################################################
## Correct all of the 3mm 128ch data
## 
cal.reset(); 
cal.setdata(msselect='FIELD_ID IN [1,2,3,4] && SPECTRAL_WINDOW_ID IN [7,8,15,16,19,20]');
cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mm.bpoly');  
cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal');
cal.correct();

###########################################################
## Do separate bandpasses for the 512ch windows
## 

cal.reset();
cal.setdata(msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID==3'); 
cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mm.bpoly');  
cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal');
cal.setsolvebandpoly(table='26-jan-2002-l02d.3mmLSBa.bpoly', 
  degamp=3, degphase=3, visnorm=F, bpnorm=F,                         
  refant=2, append=F); 
cal.solve();                                  
##
## Plot files will be:
## 26-jan-2002-l02d.3mmLSBa.bpoly.ps
## 
cal.reset();
cal.setdata(msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID==11'); 
cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mm.bpoly');  
cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal');
cal.setsolvebandpoly(table='26-jan-2002-l02d.3mmLSBb.bpoly', 
  degamp=3, degphase=3, visnorm=F, bpnorm=F,                         
  refant=2, append=F); 
cal.solve();                                  
##
## Plot files will be:
## 26-jan-2002-l02d.3mmLSBb.bpoly.ps
## 
cal.reset();
cal.setdata(msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID==4'); 
cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mm.bpoly');  
cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal');
cal.setsolvebandpoly(table='26-jan-2002-l02d.3mmUSBa.bpoly', 
  degamp=3, degphase=3, visnorm=F, bpnorm=F,                         
  refant=2, append=F); 
cal.solve();                                  
##
## Plot files will be:
## 26-jan-2002-l02d.3mmUSBa.bpoly.ps
## 
cal.reset();
cal.setdata(msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID==12'); 
cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mm.bpoly');  
cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal');
cal.setsolvebandpoly(table='26-jan-2002-l02d.3mmUSBb.bpoly', 
  degamp=3, degphase=3, visnorm=F, bpnorm=F,                         
  refant=2, append=F); 
cal.solve();                                  
##
## Plot files will be:
## 26-jan-2002-l02d.3mmUSBb.bpoly.ps
## 
###########################################################
## Correct the 3mm 512ch data
## 
cal.reset(); 
cal.setdata(msselect='FIELD_ID IN [1,2,3,4] && SPECTRAL_WINDOW_ID IN [3]');
cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mmLSBa.bpoly');  
cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal');
cal.correct();

cal.reset(); 
cal.setdata(msselect='FIELD_ID IN [1,2,3,4] && SPECTRAL_WINDOW_ID IN [11]');
cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mmLSBb.bpoly');  
cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal');
cal.correct();

cal.reset(); 
cal.setdata(msselect='FIELD_ID IN [1,2,3,4] && SPECTRAL_WINDOW_ID IN [4]');
cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mmUSBa.bpoly');  
cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal');
cal.correct();

cal.reset(); 
cal.setdata(msselect='FIELD_ID IN [1,2,3,4] && SPECTRAL_WINDOW_ID IN [12]');
cal.setapply(type='BPOLY', table='26-jan-2002-l02d.3mmUSBb.bpoly');  
cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal');
cal.correct();

###########################################################
##########      1 mm Calibration                 ##########
###########################################################

###########################################################
## Initial phase solutions:
## Use the 128ch windows [27,28,31,32]
## First the LSB channels [27,31]
##
cal.reset();
cal.setdata(mode='channel', start=4, step=1, nchan=122,   
   msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID IN [27,31]');
cal.setsolvegainspline(table='26-jan-2002-l02d.1mmLSB.gcal0', 
   mode='PHAS', preavg=0, splinetime=10000, refant=2, npointaver=10,  
   phasewrap=260, append=F);   
cal.solve();
##
## Plot files will be:
## 26-jan-2002-l02d.1mmLSB.gcal0.PHASE.ps
## 26-jan-2002-l02d.1mmLSB.gcal0.PHASE.log
##
## Now the USB channels [28,32]
cal.reset();
cal.setdata(mode='channel', start=4, step=1, nchan=122,
   msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID IN [28,32]');
cal.setsolvegainspline(table='26-jan-2002-l02d.1mmUSB.gcal0', 
   mode='PHAS', preavg=0, splinetime=10000, refant=2, npointaver=10,  
   phasewrap=260, append=F);   
cal.solve();                        
##
## Plot files will be:
## 26-jan-2002-l02d.1mmUSB.gcal0.PHASE.ps
## 26-jan-2002-l02d.1mmUSB.gcal0.PHASE.log
##
###########################################################
## Bandpass calibration:
## First for 1mm LSB
## 
cal.reset();
cal.setdata(msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID IN [27,31]'); 
cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.1mmUSB.gcal0');         
cal.setsolvebandpoly(table='26-jan-2002-l02d.1mm.bpoly', 
    degamp=6, degphase=12, visnorm=F, bpnorm=F,                         
    refant=2, append=F);  
cal.solve();                                  
##
## Plot files will be:
## 26-jan-2002-l02d.1mm.bpoly.ps
## Need to move these
shell('mv 26-jan-2002-l02d.1mm.bpoly.ps 26-jan-2002-l02d.1mmLSB.bpoly.ps')
##
## Now for 1mm USB
cal.reset();                                  
cal.setdata(msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID IN [28,32]'); 
cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.1mmUSB.gcal0');         
cal.setsolvebandpoly(table='26-jan-2002-l02d.1mm.bpoly', 
    degamp=6, degphase=12, visnorm=F, bpnorm=F,                         
    refant=2, append=T);  
cal.solve();                                  
##
## Plot files will be:
## 26-jan-2002-l02d.1mm.bpoly.ps
## Need to move these
shell('mv 26-jan-2002-l02d.1mm.bpoly.ps 26-jan-2002-l02d.1mmUSB.bpoly.ps')
##
###########################################################
## Phase calibrate 1mm LSB & USB:
## 
cal.reset(); 
cal.setdata(mode='channel', start=4, step=1, nchan=122, 
    msselect='FIELD_ID IN [1,2,3] && SPECTRAL_WINDOW_ID IN [27,28,31,32]');
cal.setapply(type='BPOLY', table='26-jan-2002-l02d.1mm.bpoly');
## 
## Doing a curve transfer of 3mm solution to 1mm to reduce phase wrapping
## incremental solution with cumulative output table
##
cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.3mm.gcal');
cal.setsolvegainspline(table='26-jan-2002-l02d.1mm.gcal',
    mode='PHAS', preavg=0, splinetime=10000,
    refant=2, npointaver=10, phasewrap=260, append=F); 
cal.solve(); 
##
## Plot files will be:
## 26-jan-2002-l02d.1mm.gcal.PHASE.ps
## 26-jan-2002-l02d.1mm.gcal.PHASE.log
##
###########################################################
## Amplitude calibration of 1mm :
## Note, we had used setjy for all 1mm calibrators
## so no need for fluxscale here
##
cal.reset(); 
cal.setdata(mode='channel', start=4, step=1, nchan=122,
    msselect='FIELD_ID IN [1,2,3] && SPECTRAL_WINDOW_ID IN [27,28,31,32]');
cal.setapply(type='BPOLY', table='26-jan-2002-l02d.1mm.bpoly');  
cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.1mm.gcal');   
cal.setsolvegainspline(table='26-jan-2002-l02d.1mm.gcal',
    mode='AMP', preavg=0.0, splinetime=10000,
    refant=2, npointaver=10, phasewrap=260); 
cal.solve();  
##
## Plot files will be:
## 26-jan-2002-l02d.1mm.gcal.AMP.ps
## 26-jan-2002-l02d.1mm.gcal.AMP.log
##
###########################################################
## Correct all of the 1mm 128ch data
## 
cal.reset(); 
cal.setdata(msselect='FIELD_ID IN [1,2,3,4] && SPECTRAL_WINDOW_ID IN [27,28,31,32]');
cal.setapply(type='BPOLY', table='26-jan-2002-l02d.1mm.bpoly');  
cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.1mm.gcal');
cal.correct();

###########################################################
## Do separate bandpasses for the 512ch windows
## 

cal.reset();
cal.setdata(msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID==23'); 
cal.setapply(type='BPOLY', table='26-jan-2002-l02d.1mm.bpoly');  
cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.1mm.gcal');
cal.setsolvebandpoly(table='26-jan-2002-l02d.1mmLSBa.bpoly', 
  degamp=3, degphase=3, visnorm=F, bpnorm=F,                         
  refant=2, append=F); 
cal.solve();                                  
##
## Plot files will be:
## 26-jan-2002-l02d.1mmLSBa.bpoly.ps
## 
cal.reset();
cal.setdata(msselect='FIELD_ID==1 && SPECTRAL_WINDOW_ID==24'); 
cal.setapply(type='BPOLY', table='26-jan-2002-l02d.1mm.bpoly');  
cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.1mm.gcal');
cal.setsolvebandpoly(table='26-jan-2002-l02d.1mmUSBa.bpoly', 
  degamp=3, degphase=3, visnorm=F, bpnorm=F,                         
  refant=2, append=F); 
cal.solve();                                  
##
## Plot files will be:
## 26-jan-2002-l02d.1mmUSBa.bpoly.ps
## 
###########################################################
## Correct the 1mm 512ch data
## 
cal.reset(); 
cal.setdata(msselect='FIELD_ID IN [1,2,3,4] && SPECTRAL_WINDOW_ID IN [23]');
cal.setapply(type='BPOLY', table='26-jan-2002-l02d.1mmLSBa.bpoly');  
cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.1mm.gcal');
cal.correct();

cal.reset(); 
cal.setdata(msselect='FIELD_ID IN [1,2,3,4] && SPECTRAL_WINDOW_ID IN [24]');
cal.setapply(type='BPOLY', table='26-jan-2002-l02d.1mmUSBa.bpoly');  
cal.setapply(type='GSPLINE', table='26-jan-2002-l02d.1mm.gcal');
cal.correct();

###########################################################
## Done with calibration
## 
cal.done();
imgr.done();

###########################################################
## Split out continuum, 128ch line, and 512ch line data
## 
myms:=ms('26-jan-2002-l02d.ms', readonly=F);
## 3mm
## averaging the 122 middle channel of the 128ch continuum spw to 1 channel
myms.split(outputms='20126+4104.26-jan-2002-l02d.3mmCONT.ms', fieldids=[4], 
    spwids=[7,8,15,16,19,20], nchan=[1], start=[4], step=[122], 
    whichcol='CORRECTED_DATA');

## Now the 6 128ch bands...
myms.split(outputms='20126+4104.26-jan-2002-l02d.3mmLSB128a.ms', fieldids=[4], 
    spwids=[7], nchan=[128], start=[1], step=[1], 
    whichcol='CORRECTED_DATA');
myms.split(outputms='20126+4104.26-jan-2002-l02d.3mmLSB128b.ms', fieldids=[4], 
    spwids=[15], nchan=[128], start=[1], step=[1], 
    whichcol='CORRECTED_DATA');
myms.split(outputms='20126+4104.26-jan-2002-l02d.3mmLSB128c.ms', fieldids=[4], 
    spwids=[19], nchan=[128], start=[1], step=[1], 
    whichcol='CORRECTED_DATA');
myms.split(outputms='20126+4104.26-jan-2002-l02d.3mmUSB128a.ms', fieldids=[4], 
    spwids=[8], nchan=[128], start=[1], step=[1], 
    whichcol='CORRECTED_DATA');
myms.split(outputms='20126+4104.26-jan-2002-l02d.3mmUSB128b.ms', fieldids=[4], 
    spwids=[16], nchan=[128], start=[1], step=[1], 
    whichcol='CORRECTED_DATA');
myms.split(outputms='20126+4104.26-jan-2002-l02d.3mmUSB128c.ms', fieldids=[4], 
    spwids=[20], nchan=[128], start=[1], step=[1], 
    whichcol='CORRECTED_DATA');

## Now the 4 512ch line bands...
myms.split(outputms='20126+4104.26-jan-2002-l02d.3mmLSB512a.ms', fieldids=[4], 
    spwids=[3], nchan=[512], start=[1], step=[1], 
    whichcol='CORRECTED_DATA');
myms.split(outputms='20126+4104.26-jan-2002-l02d.3mmLSB512b.ms', fieldids=[4], 
    spwids=[11], nchan=[512], start=[1], step=[1], 
    whichcol='CORRECTED_DATA');
myms.split(outputms='20126+4104.26-jan-2002-l02d.3mmUSB512a.ms', fieldids=[4], 
    spwids=[4], nchan=[512], start=[1], step=[1], 
    whichcol='CORRECTED_DATA');
myms.split(outputms='20126+4104.26-jan-2002-l02d.3mmUSB512b.ms', fieldids=[4], 
    spwids=[12], nchan=[512], start=[1], step=[1], 
    whichcol='CORRECTED_DATA');

## 1mm
## averaging the 122 middle channel of the 128ch continuum spw to 1 channel
myms.split(outputms='20126+4104.26-jan-2002-l02d.1mmCONT.ms', fieldids=[4], 
    spwids=[27,28,31,32], nchan=[1], start=[4], step=[122], 
    whichcol='CORRECTED_DATA');
## Now the 4 128ch line bands...
myms.split(outputms='20126+4104.26-jan-2002-l02d.1mmLSB128a.ms', fieldids=[4], 
    spwids=[27], nchan=[128], start=[1], step=[1], 
    whichcol='CORRECTED_DATA');
myms.split(outputms='20126+4104.26-jan-2002-l02d.1mmLSB128b.ms', fieldids=[4], 
    spwids=[31], nchan=[128], start=[1], step=[1], 
    whichcol='CORRECTED_DATA');
myms.split(outputms='20126+4104.26-jan-2002-l02d.1mmUSB128a.ms', fieldids=[4], 
    spwids=[28], nchan=[128], start=[1], step=[1], 
    whichcol='CORRECTED_DATA');
myms.split(outputms='20126+4104.26-jan-2002-l02d.1mmUSB128b.ms', fieldids=[4], 
    spwids=[32], nchan=[128], start=[1], step=[1], 
    whichcol='CORRECTED_DATA');
## Now the 2 512ch line bands...
myms.split(outputms='20126+4104.26-jan-2002-l02d.1mmLSB512a.ms', fieldids=[4], 
    spwids=[23], nchan=[512], start=[1], step=[1], 
    whichcol='CORRECTED_DATA');
myms.split(outputms='20126+4104.26-jan-2002-l02d.1mmUSB512a.ms', fieldids=[4], 
    spwids=[24], nchan=[512], start=[1], step=[1], 
    whichcol='CORRECTED_DATA');

myms.done();

###########################################################
###########################################################
###########################################################


###########################################################
##########    3 mm Continuum Imaging             ##########
###########################################################

imgr:=imager(filename='20126+4104.26-jan-2002-l02d.3mmCONT.ms'); 
##
## The six spectral windows we tagged as continuum 
imgr.setdata(mode='channel', nchan=1, start=1, step=1, spwid=[1:6], fieldid=1);
##
## IRAM 15m FWHM should be about 50" at 96GHz
##
imgr.setimage(nx=256, ny=256,        # Set image plane parameters
              cellx='0.25arcsec',
              celly='0.25arcsec',
              stokes='I',            # Image total intensity
              mode='mfs',            # Doesnt matter here
              nchan=1, start=1,      # 
              step=1, spwid=[1:6],
              fieldid=1);  
imgr.weight(type='briggs',           # Set up Robust weighting
            rmode='norm',
            robust=0.5);  
##
## Make dirty map
##      
imgr.makeimage(type='corrected',image='20126+4104.3mmCONT.dirty');
##
## Find peak and box around it
##
im:=image('20126+4104.3mmCONT.dirty');
im.statistics(statsout=statrec);
im.done();
print '===================================';
print 'Processed 3mm continuum dirty image';
print 'Pixel position of maximum should be [129 129 1 1]:', statrec.maxpos;
print 'Maximum dirty flux should be 0.02107Jy:', statrec.max;
print 'RMS of dirty image should be 0.00238Jy:', statrec.rms;

maxblc:= [statrec.maxpos[1]-10, statrec.maxpos[2]-10, 1, 1];
maxtrc:= [statrec.maxpos[1]+10, statrec.maxpos[2]+10, 1, 1];
imgr.boxmask(mask='maxbox', blc=maxblc, trc=maxtrc);
##
## Clean 50 iterations around peak
##
imgr.clean(algorithm='clark',          # use the Clark CLEAN algorithm
           niter=50,                   # 50 iterations 
           gain=0.1,                 
           model='20126+4104.3mmCONT.model',     
           image='20126+4104.3mmCONT.clean',
           residual='20126+4104.3mmCONT.resid',
           mask='maxbox',              # Use our box
           interactive=F);
##
## Clean 50 more iterations in inner quarter
##
shell('cp -r 20126+4104.3mmCONT.model 20126+4104.3mmCONT.model2')
imgr.clean(algorithm='clark',        
           niter=50,                   # 50 more iterations 
           gain=0.1,                 
           model='20126+4104.3mmCONT.model2',     
           image='20126+4104.3mmCONT.clean2',
           residual='20126+4104.3mmCONT.resid2',
           mask='',                    # Inner quarter (no mask set)
           interactive=F);
imgr.done();
##
## Check results
##
im:=image('20126+4104.3mmCONT.clean2');
im.statistics(statsout=statrec);
im.done();
print '===================================';
print 'Processed 3mm continuum clean image';
print 'Pixel position of maximum should be [129 129 1 1]:',statrec.maxpos;
print 'Maximum clean flux should be 0.01815Jy:',statrec.max;
im:=image('20126+4104.3mmCONT.resid2');
im.statistics(statsout=statrec);
im.done();
print 'RMS residual should be 0.00104Jy:',statrec.rms;

###########################################################
##########    1 mm Continuum Imaging             ##########
###########################################################

imgr:=imager(filename='20126+4104.26-jan-2002-l02d.1mmCONT.ms'); 
##
## The four spectral windows we tagged as continuum 
imgr.setdata(mode='channel', nchan=1, start=1, step=1, spwid=[1:4], fieldid=1);
imgr.setimage(nx=256, ny=256,        # Set image plane parameters
              cellx='0.1arcsec',
              celly='0.1arcsec',
              stokes='I',            # Image total intensity
              mode='mfs',            # Doesnt matter here
              nchan=1, start=1,      # 
              step=1, spwid=[1:4],
              fieldid=1);  
imgr.weight(type='briggs',           # Set up Robust weighting
            rmode='norm',
            robust=0.5);      
##
## Dirty image, find peak and box for clean
##
imgr.makeimage(type='corrected',image='20126+4104.1mmCONT.dirty');
im:=image('20126+4104.1mmCONT.dirty');
im.statistics(statsout=statrec);
im.done();
print '===================================';
print 'Processed 1mm continuum dirty image';
print 'Pixel position of maximum should be [130 130 1 1]:',statrec.maxpos;
print 'Maximum dirty flux should be 0.1880Jy:',statrec.max;
print 'RMS of dirty image should be 0.0330Jy:',statrec.rms;

maxblc:= [statrec.maxpos[1]-10, statrec.maxpos[2]-10, 1, 1];
maxtrc:= [statrec.maxpos[1]+10, statrec.maxpos[2]+10, 1, 1];
imgr.boxmask(mask='maxbox', blc=maxblc, trc=maxtrc);
##
## Clean in box
##
imgr.clean(algorithm='clark',          # use the Clark CLEAN algorithm
           niter=50,                   # 50 iterations 
           gain=0.1,                 
           model='20126+4104.1mmCONT.model',     
           image='20126+4104.1mmCONT.clean',
           residual='20126+4104.1mmCONT.resid',
           mask='maxbox',              # Use our box
           interactive=F);
##
## Further clean in inner quarter
##
shell('cp -r 20126+4104.1mmCONT.model 20126+4104.1mmCONT.model2')
imgr.clean(algorithm='clark',        
           niter=50,                   # 50 more iterations 
           gain=0.1,                 
           model='20126+4104.1mmCONT.model2',     
           image='20126+4104.1mmCONT.clean2',
           residual='20126+4104.1mmCONT.resid2',
           mask='',                    # Inner quarter (no mask set)
           interactive=F);
imgr.done();
##
## Check results
##
im:=image('20126+4104.1mmCONT.clean2');
im.statistics(statsout=statrec);
im.done();
print '===================================';
print 'Processed 1mm continuum clean image';
print 'Pixel position of maximum should be [130 130 1 1]:',statrec.maxpos;
print 'Maximum clean flux should be 0.1487Jy:',statrec.max;
im:=image('20126+4104.1mmCONT.resid2');
im.statistics(statsout=statrec);
im.done();
print 'RMS residual should be 0.0146:',statrec.rms;

###########################################################
##########    1 mm Line cubes (512ch)            ##########
###########################################################
##
## LSB
##
imgr:=imager(filename='20126+4104.26-jan-2002-l02d.1mmLSB512a.ms'); 
imgr.setdata(mode='channel', nchan=512, start=1, step=1, spwid=1, fieldid=1);
## Average together pairs of channels
imgr.setimage(nx=128, ny=128,      # Set image plane parameters
            cellx='0.15arcsec',
            celly='0.15arcsec',
            stokes='I',            # Image total intensity
            mode='channel',        # Cube
            nchan=256,             # Average together 2 chans
            start=1, step=2, 
            spwid=1, fieldid=1);  
imgr.weight(type='briggs',         # Set up Robust weighting
          rmode='norm',
          robust=0.5);      
imgr.makeimage(type='corrected',image='20126+4104.1mmLSB512a.dirty');
imgr.done();
##
## USB
##
imgr:=imager(filename='20126+4104.26-jan-2002-l02d.1mmUSB512a.ms'); 
imgr.setdata(mode='channel', nchan=512, start=1, step=1, spwid=1, fieldid=1);
## Average together pairs of channels
imgr.setimage(nx=128, ny=128,      # Set image plane parameters
            cellx='0.15arcsec',
            celly='0.15arcsec',
            stokes='I',            # Image total intensity
            mode='channel',        # Cube
            nchan=256,             # Average together 2 chans
            start=1, step=2, 
            spwid=1, fieldid=1);  
imgr.weight(type='briggs',         # Set up Robust weighting
          rmode='norm',
          robust=0.5);      
imgr.makeimage(type='corrected',image='20126+4104.1mmUSB512a.dirty');
imgr.done();

###########################################################
##########    3 mm Line cubes (512ch)            ##########
###########################################################
##
## LSB
##
imgr:=imager(filename='20126+4104.26-jan-2002-l02d.3mmLSB512a.ms'); 
imgr.setdata(mode='channel', nchan=512, start=1, step=1, spwid=1, fieldid=1);
## Average together pairs of channels
imgr.setimage(nx=128, ny=128,      # Set image plane parameters
            cellx='0.4arcsec',
            celly='0.4arcsec',
            stokes='I',            # Image total intensity
            mode='channel',        # Cube
            nchan=256,             # Average together 2 chans
            start=1, step=2, 
            spwid=1, fieldid=1);  
imgr.weight(type='briggs',         # Set up Robust weighting
          rmode='norm',
          robust=0.5);      
imgr.makeimage(type='corrected',image='20126+4104.3mmLSB512a.dirty');
imgr.done();

imgr:=imager(filename='20126+4104.26-jan-2002-l02d.3mmLSB512b.ms'); 
imgr.setdata(mode='channel', nchan=512, start=1, step=1, spwid=1, fieldid=1);
## Average together pairs of channels
imgr.setimage(nx=128, ny=128,      # Set image plane parameters
            cellx='0.4arcsec',
            celly='0.4arcsec',
            stokes='I',            # Image total intensity
            mode='channel',        # Cube
            nchan=256,             # Average together 2 chans
            start=1, step=2, 
            spwid=1, fieldid=1);  
imgr.weight(type='briggs',         # Set up Robust weighting
          rmode='norm',
          robust=0.5);      
imgr.makeimage(type='corrected',image='20126+4104.3mmLSB512b.dirty');
imgr.done();
##
## USB
##
imgr:=imager(filename='20126+4104.26-jan-2002-l02d.3mmUSB512a.ms'); 
imgr.setdata(mode='channel', nchan=512, start=1, step=1, spwid=1, fieldid=1);
## Average together pairs of channels
imgr.setimage(nx=128, ny=128,      # Set image plane parameters
            cellx='0.4arcsec',
            celly='0.4arcsec',
            stokes='I',            # Image total intensity
            mode='channel',        # Cube
            nchan=256,             # Average together 2 chans
            start=1, step=2, 
            spwid=1, fieldid=1);  
imgr.weight(type='briggs',         # Set up Robust weighting
          rmode='norm',
          robust=0.5);      
imgr.makeimage(type='corrected',image='20126+4104.3mmUSB512a.dirty');
imgr.done();

imgr:=imager(filename='20126+4104.26-jan-2002-l02d.3mmUSB512b.ms'); 
imgr.setdata(mode='channel', nchan=512, start=1, step=1, spwid=1, fieldid=1);
## Average together pairs of channels
imgr.setimage(nx=128, ny=128,      # Set image plane parameters
            cellx='0.4arcsec',
            celly='0.4arcsec',
            stokes='I',            # Image total intensity
            mode='channel',        # Cube
            nchan=256,             # Average together 2 chans
            start=1, step=2, 
            spwid=1, fieldid=1);  
imgr.weight(type='briggs',         # Set up Robust weighting
          rmode='norm',
          robust=0.5);      
imgr.makeimage(type='corrected',image='20126+4104.3mmUSB512b.dirty');
imgr.done();
##
## These cubes can be examined in the viewer
#dv.gui();
##
## I did not see any sign of line emission in this single day of data
##

##
## Done
##