Changeset 1391

Timestamp:

07/30/07 11:59:36 (17 years ago)

Author:

Malte Marquarding

Message:

merge from alma branch to get alma/GBT support. Commented out fluxUnit changes as they are using a chnaged interface to PKSreader/writer. Also commented out progress meter related code.

Location:

trunk

Files:

• python/asapfitter.py (modified) (6 diffs)
• python/asapmath.py (modified) (2 diffs)
• python/asapplotter.py (modified) (1 diff)
• python/scantable.py (modified) (4 diffs)
• src/STAttr.cpp (modified) (4 diffs)
• src/STDefs.h (modified) (2 diffs)
• src/STFiller.cpp (modified) (8 diffs)
• src/STFiller.h (modified) (1 diff)
• src/STFitter.cpp (modified) (1 diff)
• src/STFitter.h (modified) (2 diffs)
• src/STMath.cpp (modified) (3 diffs)
• src/STMath.h (modified) (1 diff)
• src/STMathWrapper.h (modified) (1 diff)
• src/STSelector.cpp (modified) (1 diff)
• src/STTcal.cpp (modified) (2 diffs)
• src/STWriter.cpp (modified) (4 diffs)
• src/STWriter.h (modified) (1 diff)
• src/Scantable.cpp (modified) (1 diff)
• src/Scantable.h (modified) (1 diff)
• src/ScantableWrapper.h (modified) (1 diff)
• src/python_Fitter.cpp (modified) (1 diff)
• src/python_STMath.cpp (modified) (1 diff)
• src/python_Scantable.cpp (modified) (2 diffs)

trunk/python/asapfitter.py

@@ -25 +25 @@
         self._p = None
         self._selection = None
+        self.uselinear = False
 
     def set_data(self, xdat, ydat, mask=None):
@@ -74 +75 @@
         Set the function to be fit.
         Parameters:
-            poly:    use a polynomial of the order given
+            poly:    use a polynomial of the order given with nonlinear least squares fit 
+            lpoly:   use polynomial of the order given with linear least squares fit
             gauss:   fit the number of gaussian specified
         Example:
             fitter.set_function(gauss=2) # will fit two gaussians
-            fitter.set_function(poly=3)  # will fit a 3rd order polynomial
+            fitter.set_function(poly=3)  # will fit a 3rd order polynomial via nonlinear method
+            fitter.set_function(lpoly=3)  # will fit a 3rd order polynomial via linear method
         """
         #default poly order 0
@@ -86 +89 @@
             n = kwargs.get('poly')
             self.components = [n]
+            self.uselinear = False 
+        elif kwargs.has_key('lpoly'):
+            self.fitfunc = 'poly'
+            n = kwargs.get('lpoly')
+            self.components = [n]
+            self.uselinear = True
         elif kwargs.has_key('gauss'):
             n = kwargs.get('gauss')
@@ -91 +100 @@
             self.fitfuncs = [ 'gauss' for i in range(n) ]
             self.components = [ 3 for i in range(n) ]
+            self.uselinear = False 
         else:
             msg = "Invalid function type."
@@ -146 +156 @@
             if len(fxdpar) and fxdpar.count(0) == 0:
                  raise RuntimeError,"No point fitting, if all parameters are fixed."
-            converged = self.fitter.fit()
+            if self.uselinear:
+                converged = self.fitter.lfit()
+            else:
+                converged = self.fitter.fit()
             if not converged:
                 raise RuntimeError,"Fit didn't converge."
@@ -501 +514 @@
                              array(self.data._getmask(self._fittedrow),
                                    copy=False))
-            
+                             
             ylab = self.data._get_ordinate_label()
 

trunk/python/asapmath.py

@@ -2 +2 @@
 from asap import rcParams
 from asap import print_log
+from asap import selector
 
 def average_time(*args, **kwargs):
@@ -113 +114 @@
     return s
 
+def dototalpower(calon, caloff, tcalval=0.0):
+    """
+    Do calibration for CAL on,off signals.
+    Adopted from GBTIDL dototalpower
+    Parameters:
+        calon:         the 'cal on' subintegration
+        caloff:        the 'cal off' subintegration
+        tcalval:       user supplied Tcal value
+    """
+    varlist = vars()
+    from asap._asap import stmath
+    stm = stmath()
+    stm._setinsitu(False)
+    s = scantable(stm._dototalpower(calon, caloff, tcalval))
+    s._add_history("dototalpower",varlist)
+    print_log()
+    return s
+
+def dosigref(sig, ref, smooth, tsysval=0.0, tauval=0.0):
+    """
+    Calculate a quotient (sig-ref/ref * Tsys)
+    Adopted from GBTIDL dosigref
+    Parameters:
+        sig:         on source data
+        ref:         reference data
+        smooth:      width of box car smoothing for reference
+        tsysval:     user specified Tsys (scalar only)
+        tauval:      user specified Tau (required if tsysval is set)
+    """
+    varlist = vars()
+    from asap._asap import stmath
+    stm = stmath()
+    stm._setinsitu(False)
+    s = scantable(stm._dosigref(sig, ref, smooth, tsysval, tauval))
+    s._add_history("dosigref",varlist)
+    print_log()
+    return s
+
+def calps(scantab, scannos, smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0):
+    """
+    Calibrate GBT position switched data
+    Adopted from GBTIDL getps
+    Currently calps identify the scans as position switched data if they
+    contain '_ps' in the source name. The data must contains 'CAL' signal
+    on/off in each integration. To identify 'CAL' on state, the word, 'calon'
+    need to be present in the source name field.
+    (GBT MS data reading process to scantable automatically append these
+    id names to the source names)
+
+    Parameters:
+        scantab:       scantable
+        scannos:       list of scan numbers
+        smooth:        optional box smoothing order for the reference
+                       (default is 1 = no smoothing)
+        tsysval:       optional user specified Tsys (default is 0.0,
+                       use Tsys in the data)
+        tauval:        optional user specified Tau
+        tcalval:       optional user specified Tcal (default is 0.0,
+                       use Tcal value in the data)
+    """
+    varlist = vars()
+    # check for the appropriate data
+    s = scantab.get_scan('*_ps*')
+    if s is None:
+        msg = "The input data appear to contain no position-switch mode data."
+        if rcParams['verbose']:
+            print msg
+            return
+        else:
+            raise TypeError(msg)
+    ssub = s.get_scan(scannos)
+    if ssub is None:
+        msg = "No data was found with given scan numbers!"
+        if rcParams['verbose']:
+            print msg
+            return
+        else:
+            raise TypeError(msg)
+    ssubon = ssub.get_scan('*calon')
+    ssuboff = ssub.get_scan('*[^calon]')
+    if ssubon.nrow() != ssuboff.nrow():
+        msg = "mismatch in numbers of CAL on/off scans. Cannot calibrate. Check the scan numbers."
+        if rcParams['verbose']:
+            print msg
+            return
+        else:
+            raise TypeError(msg)
+    cals = dototalpower(ssubon, ssuboff, tcalval)
+    sig = cals.get_scan('*ps')
+    ref = cals.get_scan('*psr')
+    if sig.nscan() != ref.nscan():
+        msg = "mismatch in numbers of on/off scans. Cannot calibrate. Check the scan numbers."
+        if rcParams['verbose']:
+            print msg
+            return
+        else:
+            raise TypeError(msg)
+
+    #for user supplied Tsys
+    if tsysval>0.0:
+        if tauval<=0.0:
+            msg = "Need to supply a valid tau to use the supplied Tsys"
+            if rcParams['verbose']:
+                print msg
+                return
+            else:
+                raise TypeError(msg)
+        else:
+            sig.recalc_azel()
+            ref.recalc_azel()
+            #msg = "Use of user specified Tsys is not fully implemented yet."
+            #if rcParams['verbose']:
+            #    print msg
+            #    return
+            #else:
+            #    raise TypeError(msg)
+            # use get_elevation to get elevation and
+            # calculate a scaling factor using the formula
+            # -> tsys use to dosigref
+
+    #ress = dosigref(sig, ref, smooth, tsysval)
+    ress = dosigref(sig, ref, smooth, tsysval, tauval)
+    ress._add_history("calps", varlist)
+    print_log()
+    return ress
+
+def calnod(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0):
+    """
+    Do full (but a pair of scans at time) processing of GBT Nod data
+    calibration.
+    Adopted from  GBTIDL's getnod
+    Parameters:
+        scantab:     scantable
+        scannos:     a pair of scan numbers, or the first scan number of the pair
+        smooth:      box car smoothing order
+        tsysval:     optional user specified Tsys value
+        tauval:      optional user specified tau value (not implemented yet)
+        tcalval:     optional user specified Tcal value
+    """
+    varlist = vars()
+    from asap._asap import stmath
+    stm = stmath()
+    stm._setinsitu(False)
+
+    # check for the appropriate data
+    s = scantab.get_scan('*_nod*')
+    if s is None:
+        msg = "The input data appear to contain no Nod observing mode data."
+        if rcParams['verbose']:
+            print msg
+            return
+        else:
+            raise TypeError(msg)
+
+    # need check correspondance of each beam with sig-ref ...
+    # check for timestamps, scan numbers, subscan id (not available in
+    # ASAP data format...). Assume 1st scan of the pair (beam 0 - sig
+    # and beam 1 - ref...)
+    # First scan number of paired scans or list of pairs of
+    # scan numbers (has to have even number of pairs.)
+
+    #data splitting
+    scan1no = scan2no = 0
+
+    if len(scannos)==1:
+        scan1no = scannos[0]
+        scan2no = scannos[0]+1
+        pairScans = [scan1no, scan2no]
+    else:
+        #if len(scannos)>2:
+        #    msg = "calnod can only process a pair of nod scans at time."
+        #    if rcParams['verbose']:
+        #        print msg
+        #        return
+        #    else:
+        #        raise TypeError(msg)
+        #
+        #if len(scannos)==2:
+        #    scan1no = scannos[0]
+        #    scan2no = scannos[1]
+        pairScans = list(scannos)
+
+    if tsysval>0.0:
+        if tauval<=0.0:
+            msg = "Need to supply a valid tau to use the supplied Tsys"
+            if rcParams['verbose']:
+                print msg
+                return
+            else:
+                raise TypeError(msg)
+        else:
+            scantab.recalc_azel()
+    resspec = scantable(stm._donod(scantab, pairScans, smooth, tsysval,tauval,tcalval))
+    resspec._add_history("calnod",varlist)
+    print_log()
+    return resspec
+
+def calfs(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0):
+    """
+    Calibrate GBT frequency switched data.
+    Adopted from GBTIDL getfs.
+    Currently calfs identify the scans as frequency switched data if they
+    contain '_fs' in the source name. The data must contains 'CAL' signal
+    on/off in each integration. To identify 'CAL' on state, the word, 'calon'
+    need to be present in the source name field.
+    (GBT MS data reading via scantable automatically append these
+    id names to the source names)
+
+    Parameters:
+        scantab:       scantable
+        scannos:       list of scan numbers
+        smooth:        optional box smoothing order for the reference
+                       (default is 1 = no smoothing)
+        tsysval:       optional user specified Tsys (default is 0.0,
+                       use Tsys in the data)
+        tauval:        optional user specified Tau
+    """
+    varlist = vars()
+    from asap._asap import stmath
+    stm = stmath()
+    stm._setinsitu(False)
+
+#    check = scantab.get_scan('*_fs*')
+#    if check is None:
+#        msg = "The input data appear to contain no Nod observing mode data."
+#        if rcParams['verbose']:
+#            print msg
+#            return
+#        else:
+#            raise TypeError(msg)
+    s = scantab.get_scan(scannos)
+    del scantab
+
+    resspec = scantable(stm._dofs(s, scannos, smooth, tsysval,tauval,tcalval))
+    resspec._add_history("calfs",varlist)
+    print_log()
+    return resspec
+
 def simple_math(left, right, op='add', tsys=True):
     """

trunk/python/asapplotter.py

@@ -741 +741 @@
         return userlabel or d[mode]
 
+    def plotazel(self, scan=None):
+        """
+        plot azimuth and elevation  versus time of a scantable
+        """
+        import pylab as PL
+        from matplotlib.dates import DateFormatter, timezone, HourLocator, MinuteLocator, DayLocator
+        from matplotlib.ticker import MultipleLocator
+        from matplotlib.numerix import array, pi
+        self._data = scan
+        dates = self._data.get_time()
+        t = PL.date2num(dates)
+        tz = timezone('UTC')
+        PL.cla()
+        PL.ioff()
+        PL.clf()
+        tdel = max(t) - min(t)
+        ax = PL.subplot(2,1,1)
+        el = array(self._data.get_elevation())*180./pi
+        PL.ylabel('El [deg.]')
+        dstr = dates[0].strftime('%Y/%m/%d')
+        if tdel > 1.0:
+            dstr2 = dates[len(dates)-1].strftime('%Y/%m/%d')
+            dstr = dstr + " - " + dstr2
+            majloc = DayLocator()
+            minloc = HourLocator(range(0,23,12))
+            timefmt = DateFormatter("%b%d")
+        else:
+            timefmt = DateFormatter('%H')
+            majloc = HourLocator()
+            minloc = MinuteLocator(20)
+        PL.title(dstr)
+        PL.plot_date(t,el,'b,', tz=tz)
+        #ax.grid(True)
+        ax.yaxis.grid(True)
+        yloc = MultipleLocator(30)
+        ax.set_ylim(0,90)
+        ax.xaxis.set_major_formatter(timefmt)
+        ax.xaxis.set_major_locator(majloc)
+        ax.xaxis.set_minor_locator(minloc)
+        ax.yaxis.set_major_locator(yloc)
+        if tdel > 1.0:
+            labels = ax.get_xticklabels()
+        #    PL.setp(labels, fontsize=10, rotation=45)
+            PL.setp(labels, fontsize=10)
+        # Az plot
+        az = array(self._data.get_azimuth())*180./pi
+        if min(az) < 0:
+            for irow in range(len(az)):
+                if az[irow] < 0: az[irow] += 360.0
+
+        ax = PL.subplot(2,1,2)
+        PL.xlabel('Time (UT)')
+        PL.ylabel('Az [deg.]')
+        PL.plot_date(t,az,'b,', tz=tz)
+        ax.set_ylim(0,360)
+        #ax.grid(True)
+        ax.yaxis.grid(True)
+        #hfmt = DateFormatter('%H')
+        #hloc = HourLocator()
+        yloc = MultipleLocator(60)
+        ax.xaxis.set_major_formatter(timefmt)
+        ax.xaxis.set_major_locator(majloc)
+        ax.xaxis.set_minor_locator(minloc)
+        ax.yaxis.set_major_locator(yloc)
+        if tdel > 1.0:
+            labels = ax.get_xticklabels()
+            PL.setp(labels, fontsize=10)
+        PL.ion()
+        PL.draw()
+
+    def plotpointing(self, scan=None):
+        """
+        plot telescope pointings
+        """
+        import pylab as PL
+        from matplotlib.dates import DateFormatter, timezone
+        from matplotlib.ticker import MultipleLocator
+        from matplotlib.numerix import array, pi, zeros
+        self._data = scan
+        dir = array(self._data.get_directionval()).transpose()
+        ra = dir[0]*180./pi
+        dec = dir[1]*180./pi
+        PL.cla()
+        PL.ioff()
+        PL.clf()
+        ax = PL.axes([0.1,0.1,0.8,0.8])
+        ax = PL.axes([0.1,0.1,0.8,0.8])
+        ax.set_aspect('equal')
+        PL.plot(ra,dec, 'b,')
+        PL.xlabel('RA [deg.]')
+        PL.ylabel('Declination [deg.]')
+        PL.title('Telescope pointings')
+        [xmin,xmax,ymin,ymax] = PL.axis()
+        PL.axis([xmax,xmin,ymin,ymax])
+        PL.ion()
+        PL.draw()
+

trunk/python/scantable.py

@@ -516 +516 @@
         """
         return self._get_column(self._getdirection, row)
+
+    def get_directionval(self, row=-1):
+        """
+        Get a list of Positions on the sky (direction) for the observations.
+        Return a float for each integration in the scantable.
+        Parameters:
+            row:    row no of integration. Default -1 return all rows
+        Example:
+            none
+        """
+        return self._get_column(self._getdirectionvec, row)
+
 
     def set_unit(self, unit='channel'):
@@ -1215 +1227 @@
 
 
-    def poly_baseline(self, mask=None, order=0, plot=False, insitu=None):
+    def poly_baseline(self, mask=None, order=0, plot=False, uselin=False, insitu=None):
         """
         Return a scan which has been baselined (all rows) by a polynomial.
@@ -1224 +1236 @@
                         indivual fit has to be approved, by typing 'y'
                         or 'n'
+            uselin:     use linear polynomial fit
             insitu:     if False a new scantable is returned.
                         Otherwise, the scaling is done in-situ
@@ -1240 +1253 @@
             f = fitter()
             f.set_scan(self, mask)
-            f.set_function(poly=order)
+            #f.set_function(poly=order)
+            if uselin:
+                f.set_function(lpoly=order)
+            else:
+                f.set_function(poly=order)
             s = f.auto_fit(insitu, plot=plot)
             s._add_history("poly_baseline", varlist)

trunk/src/STAttr.cpp

@@ -73 +73 @@
    Float D = 1.0;
    switch (inst) {
+      case ALMA:
+        {
+           D = 12.0;
+        }
+        break;
+
       case ATMOPRA:
         {
@@ -92 +98 @@
         {
            D = 30.0;
+        }
+        break;
+      case GBT:
+        {
+           D = 104.9;
         }
         break;
@@ -338 +349 @@
   if (t==String("DSS-43")) {
     inst = TIDBINBILLA;
+  } else if (t==String("ALMA")) {
+    inst = ALMA;
   } else if (t==String("ATPKSMB")) {
     inst = ATPKSMB;
@@ -346 +359 @@
   } else if (t==String("CEDUNA")) {
     inst = CEDUNA;
+  } else if (t==String("GBT")) {
+    inst = GBT;
   } else if (t==String("HOBART")) {
     inst = HOBART;

trunk/src/STDefs.h

@@ -41 +41 @@
 
   enum Instrument {UNKNOWNINST=0,
+                   ALMA,
                    ATPKSMB,
                    ATPKSHOH,
@@ -46 +47 @@
                    TIDBINBILLA,
                    CEDUNA,
+                   GBT,
                    HOBART,
                    N_INSTRUMENTS};

trunk/src/STFiller.cpp

@@ -26 +26 @@
 
 #include <atnf/PKSIO/PKSreader.h>
+//#include <casa/System/ProgressMeter.h>
+
 
 #include "STDefs.h"
@@ -134 +136 @@
   }
   // Determine Telescope and set brightness unit
+
 
   Bool throwIt = False;
@@ -184 +187 @@
   reader_->select(beams, ifs, start, end, ref, True, haveXPol_[0], False);
   table_->setHeader(*header_);
+  //For MS, add the location of POINTING of the input MS so one get
+  //pointing data from there, if necessary.
+  //Also find nrow in MS 
+  nInDataRow = 0;
+  if (format == "MS2") {
+    Path datapath(inName); 
+    String ptTabPath = datapath.absoluteName();
+    Table inMS(ptTabPath);
+    nInDataRow = inMS.nrow();
+    ptTabPath.append("/POINTING");
+    table_->table().rwKeywordSet().define("POINTING", ptTabPath);
+    if ((header_->antennaname).matches("GBT")) {
+      String GOTabPath = datapath.absoluteName();
+      GOTabPath.append("/GBT_GO");
+      table_->table().rwKeywordSet().define("GBT_GO", GOTabPath);
+    }
+  }
+     
 }
 
@@ -209 +230 @@
   Complex         xCalFctr;
   Vector<Complex> xPol;
+  Double min = 0.0;
+  Double max = nInDataRow;
+  //ProgressMeter fillpm(min, max, "Data importing progress");
+  int n = 0;
   while ( status == 0 ) {
     status = reader_->read(scanNo, cycleNo, mjd, interval, fieldName,
@@ -220 +245 @@
                           spectra, flagtra, xCalFctr, xPol);
     if ( status != 0 ) break;
+    n += 1;
+    
     Regex filterrx(".*[SL|PA]$");
     Regex obsrx("^AT.+");
@@ -250 +277 @@
     RecordFieldPtr<String> srcnCol(rec, "SRCNAME");
     RecordFieldPtr<Int> srctCol(rec, "SRCTYPE");
+    RecordFieldPtr<String> fieldnCol(rec, "FIELDNAME");
+    *fieldnCol = fieldName;
     // try to auto-identify if it is on or off.
     Regex rx(".*[e|w|_R]$");
@@ -340 +369 @@
       row.put(table_->table().nrow()-1, rec);
     }
+    //fillpm._update(n);
   }
   if (status > 0) {
@@ -345 +375 @@
     throw(AipsError("Reading error occured, data possibly corrupted."));
   }
+  //fillpm.done();
   return status;
 }

trunk/src/STFiller.h

@@ -99 +99 @@
   casa::String filename_;
   casa::CountedPtr< Scantable > table_;
-  casa::Int nIF_, nBeam_, nPol_, nChan_;
+  casa::Int nIF_, nBeam_, nPol_, nChan_, nInDataRow;
   casa::uInt ifOffset_, beamOffset_;
   casa::Vector<casa::Bool> haveXPol_;

trunk/src/STFitter.cpp

@@ -329 +329 @@
 }
 
+bool Fitter::lfit() {
+  LinearFit<Float> fitter;
+  CompoundFunction<Float> func;
+
+  uInt n = funcs_.nelements();
+  for (uInt i=0; i<n; ++i) {
+    func.addFunction(*funcs_[i]);
+  }
+
+  fitter.setFunction(func);
+  //fitter.setMaxIter(50+n*10);
+  // Convergence criterium
+  //fitter.setCriteria(0.001);
+
+  // Fit
+  Vector<Float> sigma(x_.nelements());
+  sigma = 1.0;
+
+  parameters_.resize();
+  parameters_ = fitter.fit(x_, y_, sigma, &m_);
+  std::vector<float> ps;
+  parameters_.tovector(ps);
+  setParameters(ps);
+
+  error_.resize();
+  error_ = fitter.errors();
+
+  chisquared_ = fitter.getChi2();
+
+  residual_.resize();
+  residual_ =  y_;
+  fitter.residual(residual_,x_);
+  // use fitter.residual(model=True) to get the model
+  thefit_.resize(x_.nelements());
+  fitter.residual(thefit_,x_,True);
+  return true;
+}
 
 std::vector<float> Fitter::evaluate(int whichComp) const

trunk/src/STFitter.h

@@ -27 +27 @@
 //#                        AUSTRALIA
 //#
-//# $Id:$
+//# $Id$
 //#---------------------------------------------------------------------------
 #ifndef STFITTER_H
@@ -64 +64 @@
   void reset();
   bool fit();
+  // Fit via linear method
+  bool lfit();
   bool computeEstimate();
 

trunk/src/STMath.cpp

@@ -48 +48 @@
 #include "STAttr.h"
 #include "STMath.h"
+#include "STSelector.h"
 
 using namespace casa;
@@ -536 +537 @@
 }
 
+// dototalpower (migration of GBTIDL procedure dototalpower.pro)
+// calibrate the CAL on-off pair. It calculate Tsys and average CAL on-off subintegrations
+// do it for each cycles in a specific scan.
+CountedPtr< Scantable > STMath::dototalpower( const CountedPtr< Scantable >& calon,
+                                              const CountedPtr< Scantable >& caloff, Float tcal )
+{
+if ( ! calon->conformant(*caloff) ) {
+    throw(AipsError("'CAL on' and 'CAL off' scantables are not conformant."));
+  }
+  setInsitu(false);
+  CountedPtr< Scantable > out = getScantable(caloff, false);
+  Table& tout = out->table();
+  const Table& tcon = calon->table();
+  Vector<Float> tcalout;
+  Vector<Float> tcalout2;  //debug
+
+  if ( tout.nrow() != tcon.nrow() ) {
+    throw(AipsError("Mismatch in number of rows to form cal on - off pair."));
+  }
+  // iteration by scanno or cycle no.
+  TableIterator sit(tout, "SCANNO");
+  TableIterator s2it(tcon, "SCANNO");
+  while ( !sit.pastEnd() ) {
+    Table toff = sit.table();
+    TableRow row(toff);
+    Table t = s2it.table();
+    ScalarColumn<Double> outintCol(toff, "INTERVAL");
+    ArrayColumn<Float> outspecCol(toff, "SPECTRA");
+    ArrayColumn<Float> outtsysCol(toff, "TSYS");
+    ArrayColumn<uChar> outflagCol(toff, "FLAGTRA");
+    ROScalarColumn<uInt> outtcalIdCol(toff, "TCAL_ID");
+    ROScalarColumn<uInt> outpolCol(toff, "POLNO");
+    ROScalarColumn<Double> onintCol(t, "INTERVAL");
+    ROArrayColumn<Float> onspecCol(t, "SPECTRA");
+    ROArrayColumn<Float> ontsysCol(t, "TSYS");
+    ROArrayColumn<uChar> onflagCol(t, "FLAGTRA");
+    //ROScalarColumn<uInt> ontcalIdCol(t, "TCAL_ID");
+
+    for (uInt i=0; i < toff.nrow(); ++i) {
+      //skip these checks -> assumes the data order are the same between the cal on off pairs
+      //
+      Vector<Float> specCalon, specCaloff;
+      // to store scalar (mean) tsys
+      Vector<Float> tsysout(1);
+      uInt tcalId, polno;
+      Double offint, onint;
+      outpolCol.get(i, polno);
+      outspecCol.get(i, specCaloff);
+      onspecCol.get(i, specCalon);
+      Vector<uChar> flagCaloff, flagCalon;
+      outflagCol.get(i, flagCaloff);
+      onflagCol.get(i, flagCalon);
+      outtcalIdCol.get(i, tcalId);
+      outintCol.get(i, offint);
+      onintCol.get(i, onint);
+      // caluculate mean Tsys
+      uInt nchan = specCaloff.nelements();
+      // percentage of edge cut off
+      uInt pc = 10;
+      uInt bchan = nchan/pc;
+      uInt echan = nchan-bchan;
+
+      Slicer chansl(IPosition(1,bchan-1), IPosition(1,echan-1), IPosition(1,1),Slicer::endIsLast);
+      Vector<Float> testsubsp = specCaloff(chansl);
+      MaskedArray<Float> spoff = maskedArray( specCaloff(chansl),flagCaloff(chansl) );
+      MaskedArray<Float> spon = maskedArray( specCalon(chansl),flagCalon(chansl) );
+      MaskedArray<Float> spdiff = spon-spoff;
+      uInt noff = spoff.nelementsValid();
+      //uInt non = spon.nelementsValid();
+      uInt ndiff = spdiff.nelementsValid();
+      Float meantsys;
+
+/**
+      Double subspec, subdiff;
+      uInt usednchan;
+      subspec = 0;
+      subdiff = 0;
+      usednchan = 0;
+      for(uInt k=(bchan-1); k<echan; k++) {
+        subspec += specCaloff[k];
+        subdiff += static_cast<Double>(specCalon[k]-specCaloff[k]);
+        ++usednchan;
+      }
+**/
+      // get tcal if input tcal <= 0
+      String tcalt;
+      Float tcalUsed;
+      tcalUsed = tcal;
+      if ( tcal <= 0.0 ) {
+        caloff->tcal().getEntry(tcalt, tcalout, tcalId);
+        if (polno<=3) {
+          tcalUsed = tcalout[polno];
+        }
+        else {
+          tcalUsed = tcalout[0];
+        }
+      }
+
+      Float meanoff;
+      Float meandiff;
+      if (noff && ndiff) {
+         //Debug
+         //if(noff!=ndiff) cerr<<"noff and ndiff is not equal"<<endl;
+         meanoff = sum(spoff)/noff;
+         meandiff = sum(spdiff)/ndiff;
+         meantsys= (meanoff/meandiff )*tcalUsed + tcalUsed/2;
+      }
+      else {
+         meantsys=1;
+      }
+
+      tsysout[0] = Float(meantsys);
+      MaskedArray<Float> mcaloff = maskedArray(specCaloff, flagCaloff);
+      MaskedArray<Float> mcalon = maskedArray(specCalon, flagCalon);
+      MaskedArray<Float> sig =   Float(0.5) * (mcaloff + mcalon);
+      //uInt ncaloff = mcaloff.nelementsValid();
+      //uInt ncalon = mcalon.nelementsValid();
+
+      outintCol.put(i, offint+onint);
+      outspecCol.put(i, sig.getArray());
+      outflagCol.put(i, flagsFromMA(sig));
+      outtsysCol.put(i, tsysout);
+    }
+    ++sit;
+    ++s2it;
+  }
+  return out;
+}
+
+//dosigref - migrated from GBT IDL's dosigref.pro, do calibration of position switch
+// observatiions.
+// input: sig and ref scantables, and an optional boxcar smoothing width(default width=0,
+//        no smoothing).
+// output: resultant scantable [= (sig-ref/ref)*tsys]
+CountedPtr< Scantable > STMath::dosigref( const CountedPtr < Scantable >& sig,
+                                          const CountedPtr < Scantable >& ref,
+                                          int smoothref,
+                                          casa::Float tsysv,
+                                          casa::Float tau )
+{
+if ( ! ref->conformant(*sig) ) {
+    throw(AipsError("'sig' and 'ref' scantables are not conformant."));
+  }
+  setInsitu(false);
+  CountedPtr< Scantable > out = getScantable(sig, false);
+  CountedPtr< Scantable > smref;
+  if ( smoothref > 1 ) {
+    float fsmoothref = static_cast<float>(smoothref);
+    std::string inkernel = "boxcar";
+    smref = smooth(ref, inkernel, fsmoothref );
+    ostringstream oss;
+    oss<<"Applied smoothing of "<<fsmoothref<<" on the reference."<<endl;
+    pushLog(String(oss));
+  }
+  else {
+    smref = ref;
+  }
+  Table& tout = out->table();
+  const Table& tref = smref->table();
+  if ( tout.nrow() != tref.nrow() ) {
+    throw(AipsError("Mismatch in number of rows to form on-source and reference pair."));
+  }
+  // iteration by scanno? or cycle no.
+  TableIterator sit(tout, "SCANNO");
+  TableIterator s2it(tref, "SCANNO");
+  while ( !sit.pastEnd() ) {
+    Table ton = sit.table();
+    Table t = s2it.table();
+    ScalarColumn<Double> outintCol(ton, "INTERVAL");
+    ArrayColumn<Float> outspecCol(ton, "SPECTRA");
+    ArrayColumn<Float> outtsysCol(ton, "TSYS");
+    ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
+    ArrayColumn<Float> refspecCol(t, "SPECTRA");
+    ROScalarColumn<Double> refintCol(t, "INTERVAL");
+    ROArrayColumn<Float> reftsysCol(t, "TSYS");
+    ArrayColumn<uChar> refflagCol(t, "FLAGTRA");
+    ROScalarColumn<Float> refelevCol(t, "ELEVATION");
+    for (uInt i=0; i < ton.nrow(); ++i) {
+
+      Double onint, refint;
+      Vector<Float> specon, specref;
+      // to store scalar (mean) tsys
+      Vector<Float> tsysref;
+      outintCol.get(i, onint);
+      refintCol.get(i, refint);
+      outspecCol.get(i, specon);
+      refspecCol.get(i, specref);
+      Vector<uChar> flagref, flagon;
+      outflagCol.get(i, flagon);
+      refflagCol.get(i, flagref);
+      reftsysCol.get(i, tsysref);
+
+      Float tsysrefscalar;
+      if ( tsysv > 0.0 ) {
+        ostringstream oss;
+        Float elev;
+        refelevCol.get(i, elev);
+        oss << "user specified Tsys = " << tsysv;
+        // do recalc elevation if EL = 0
+        if ( elev == 0 ) {
+          throw(AipsError("EL=0, elevation data is missing."));
+        } else {
+          if ( tau <= 0.0 ) {
+            throw(AipsError("Valid tau is not supplied."));
+          } else {
+            tsysrefscalar = tsysv * exp(tau/elev);
+          }
+        }
+        oss << ", corrected (for El) tsys= "<<tsysrefscalar;
+        pushLog(String(oss));
+      }
+      else {
+        tsysrefscalar = tsysref[0];
+      }
+      //get quotient spectrum
+      MaskedArray<Float> mref = maskedArray(specref, flagref);
+      MaskedArray<Float> mon = maskedArray(specon, flagon);
+      MaskedArray<Float> specres =   tsysrefscalar*((mon - mref)/mref);
+      Double resint = onint*refint*smoothref/(onint+refint*smoothref);
+
+      //Debug
+      //cerr<<"Tsys used="<<tsysrefscalar<<endl;
+      // fill the result, replay signal tsys by reference tsys
+      outintCol.put(i, resint);
+      outspecCol.put(i, specres.getArray());
+      outflagCol.put(i, flagsFromMA(specres));
+      outtsysCol.put(i, tsysref);
+    }
+    ++sit;
+    ++s2it;
+  }
+  return out;
+}
+
+CountedPtr< Scantable > STMath::donod(const casa::CountedPtr<Scantable>& s,
+                                     const std::vector<int>& scans,
+                                     int smoothref,
+                                     casa::Float tsysv,
+                                     casa::Float tau,
+                                     casa::Float tcal )
+
+{
+  setInsitu(false);
+  STSelector sel;
+  std::vector<int> scan1, scan2, beams;
+  std::vector< vector<int> > scanpair;
+  std::vector<string> calstate;
+  String msg;
+
+  CountedPtr< Scantable > s1b1on, s1b1off, s1b2on, s1b2off;
+  CountedPtr< Scantable > s2b1on, s2b1off, s2b2on, s2b2off;
+
+  std::vector< CountedPtr< Scantable > > sctables;
+  sctables.push_back(s1b1on);
+  sctables.push_back(s1b1off);
+  sctables.push_back(s1b2on);
+  sctables.push_back(s1b2off);
+  sctables.push_back(s2b1on);
+  sctables.push_back(s2b1off);
+  sctables.push_back(s2b2on);
+  sctables.push_back(s2b2off);
+
+  //check scanlist
+  int n=s->checkScanInfo(scans);
+  if (n==1) {
+     throw(AipsError("Incorrect scan pairs. "));
+  }
+
+  // Assume scans contain only a pair of consecutive scan numbers.
+  // It is assumed that first beam, b1,  is on target.
+  // There is no check if the first beam is on or not.
+  if ( scans.size()==1 ) {
+    scan1.push_back(scans[0]);
+    scan2.push_back(scans[0]+1);
+  } else if ( scans.size()==2 ) {
+   scan1.push_back(scans[0]);
+   scan2.push_back(scans[1]);
+  } else {
+    if ( scans.size()%2 == 0 ) {
+      for (uInt i=0; i<scans.size(); i++) {
+        if (i%2 == 0) {
+          scan1.push_back(scans[i]);
+        }
+        else {
+          scan2.push_back(scans[i]);
+        }
+      }
+    } else {
+      throw(AipsError("Odd numbers of scans, cannot form pairs."));
+    }
+  }
+  scanpair.push_back(scan1);
+  scanpair.push_back(scan2);
+  calstate.push_back("*calon");
+  calstate.push_back("*[^calon]");
+  CountedPtr< Scantable > ws = getScantable(s, false);
+  uInt l=0;
+  while ( l < sctables.size() ) {
+    for (uInt i=0; i < 2; i++) {
+      for (uInt j=0; j < 2; j++) {
+        for (uInt k=0; k < 2; k++) {
+          sel.reset();
+          sel.setScans(scanpair[i]);
+          sel.setName(calstate[k]);
+          beams.clear();
+          beams.push_back(j);
+          sel.setBeams(beams);
+          ws->setSelection(sel);
+          sctables[l]= getScantable(ws, false);
+          l++;
+        }
+      }
+    }
+  }
+
+  // replace here by splitData or getData functionality
+  CountedPtr< Scantable > sig1;
+  CountedPtr< Scantable > ref1;
+  CountedPtr< Scantable > sig2;
+  CountedPtr< Scantable > ref2;
+  CountedPtr< Scantable > calb1;
+  CountedPtr< Scantable > calb2;
+
+  msg=String("Processing dototalpower for subset of the data");
+  ostringstream oss1;
+  oss1 << msg  << endl;
+  pushLog(String(oss1));
+  // Debug for IRC CS data
+  //float tcal1=7.0;
+  //float tcal2=4.0;
+  sig1 = dototalpower(sctables[0], sctables[1], tcal=tcal);
+  ref1 = dototalpower(sctables[2], sctables[3], tcal=tcal);
+  ref2 = dototalpower(sctables[4], sctables[5], tcal=tcal);
+  sig2 = dototalpower(sctables[6], sctables[7], tcal=tcal);
+
+  // correction of user-specified tsys for elevation here
+
+  // dosigref calibration
+  msg=String("Processing dosigref for subset of the data");
+  ostringstream oss2;
+  oss2 << msg  << endl;
+  pushLog(String(oss2));
+  calb1=dosigref(sig1,ref2,smoothref,tsysv,tau);
+  calb2=dosigref(sig2,ref1,smoothref,tsysv,tau);
+
+  // iteration by scanno or cycle no.
+  Table& tcalb1 = calb1->table();
+  Table& tcalb2 = calb2->table();
+  TableIterator sit(tcalb1, "SCANNO");
+  TableIterator s2it(tcalb2, "SCANNO");
+  while ( !sit.pastEnd() ) {
+    Table t1 = sit.table();
+    Table t2= s2it.table();
+    ArrayColumn<Float> outspecCol(t1, "SPECTRA");
+    ArrayColumn<Float> outtsysCol(t1, "TSYS");
+    ArrayColumn<uChar> outflagCol(t1, "FLAGTRA");
+    ScalarColumn<Double> outintCol(t1, "INTERVAL");
+    ArrayColumn<Float> t2specCol(t2, "SPECTRA");
+    ROArrayColumn<Float> t2tsysCol(t2, "TSYS");
+    ArrayColumn<uChar> t2flagCol(t2, "FLAGTRA");
+    ROScalarColumn<Double> t2intCol(t2, "INTERVAL");
+    for (uInt i=0; i < t1.nrow(); ++i) {
+      Vector<Float> spec1, spec2;
+      // to store scalar (mean) tsys
+      Vector<Float> tsys1, tsys2;
+      Vector<uChar> flag1, flag2;
+      Double tint1, tint2;
+      outspecCol.get(i, spec1);
+      t2specCol.get(i, spec2);
+      outflagCol.get(i, flag1);
+      t2flagCol.get(i, flag2);
+      outtsysCol.get(i, tsys1);
+      t2tsysCol.get(i, tsys2);
+      outintCol.get(i, tint1);
+      t2intCol.get(i, tint2);
+      // average
+      // assume scalar tsys for weights
+      Float wt1, wt2, tsyssq1, tsyssq2;
+      tsyssq1 = tsys1[0]*tsys1[0];
+      tsyssq2 = tsys2[0]*tsys2[0];
+      wt1 = Float(tint1)/tsyssq1;
+      wt2 = Float(tint2)/tsyssq2;
+      Float invsumwt=1/(wt1+wt2);
+      MaskedArray<Float> mspec1 = maskedArray(spec1, flag1);
+      MaskedArray<Float> mspec2 = maskedArray(spec2, flag2);
+      MaskedArray<Float> avspec =  invsumwt * (wt1*mspec1 + wt2*mspec2);
+      //Array<Float> avtsys =  Float(0.5) * (tsys1 + tsys2);
+      // cerr<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<endl;
+      tsys1[0] = sqrt(tsyssq1 + tsyssq2);
+      Array<Float> avtsys =  tsys1;
+
+      outspecCol.put(i, avspec.getArray());
+      outflagCol.put(i, flagsFromMA(avspec));
+      outtsysCol.put(i, avtsys);
+    }
+    ++sit;
+    ++s2it;
+  }
+  return calb1;
+}
+
+//GBTIDL version of frequency switched data calibration
+CountedPtr< Scantable > STMath::dofs( const CountedPtr< Scantable >& s,
+                                      const std::vector<int>& scans,
+                                      int smoothref,
+                                      casa::Float tsysv,
+                                      casa::Float tau,
+                                      casa::Float tcal )
+{
+
+
+  STSelector sel;
+  CountedPtr< Scantable > ws = getScantable(s, false);
+  CountedPtr< Scantable > sig, sigwcal, ref, refwcal;
+  CountedPtr< Scantable > calsig, calref, out;
+
+  //split the data
+  sel.setName("*_fs");
+  ws->setSelection(sel);
+  sig = getScantable(ws,false);
+  sel.reset();
+  sel.setName("*_fs_calon");
+  ws->setSelection(sel);
+  sigwcal = getScantable(ws,false);
+  sel.reset();
+  sel.setName("*_fsr");
+  ws->setSelection(sel);
+  ref = getScantable(ws,false);
+  sel.reset();
+  sel.setName("*_fsr_calon");
+  ws->setSelection(sel);
+  refwcal = getScantable(ws,false);
+
+  calsig = dototalpower(sigwcal, sig, tcal=tcal);
+  calref = dototalpower(refwcal, ref, tcal=tcal);
+
+  out=dosigref(calsig,calref,smoothref,tsysv,tau);
+
+  return out;
+}
+
 
 CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
@@ -1250 +1692 @@
   vec = 0;
   pols->table_.rwKeywordSet().define("nPol", Int(1));
-  pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
+  //pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
+  pols->table_.rwKeywordSet().define("POLTYPE", in->getPolType());
   std::vector<CountedPtr<Scantable> > vpols;
   vpols.push_back(pols);

trunk/src/STMath.h

@@ -132 +132 @@
                                         bool preserve = true );
 
+  /**
+    * Calibrate total power scans (translated from GBTIDL)
+    * @param calon uncalibrated Scantable with CAL noise signal
+    * @param caloff uncalibrated Scantable with no CAL signal
+    * @param tcal optional scalar Tcal, CAL temperature (K)
+    * @return casa::CountedPtr<Scantable> which holds a calibrated Scantable
+    * (spectrum - average of the two CAL on and off spectra;
+    * tsys - mean Tsys = <caloff>*Tcal/<calon-caloff> + Tcal/2)
+    */
+  casa::CountedPtr<Scantable> dototalpower( const casa::CountedPtr<Scantable>& calon,
+                                            const casa::CountedPtr<Scantable>& caloff,
+                                            casa::Float tcal=1.0 );
+
+  /**
+    * Combine signal and reference scans (translated from GBTIDL)
+    * @param sig Scantable which contains signal scans
+    * @param ref Scantable which contains reference scans
+    * @param smoothref optional Boxcar smooth width of the reference scans
+    * default: no smoothing (=1)
+    * @param tsysv optional scalar Tsys value at the zenith, required to
+    * set tau, as well
+    * @param tau optional scalar Tau value
+    * @return casa::CountedPtr<Scantable> which holds combined scans
+    * (spectrum = (sig-ref)/ref * Tsys )
+    */
+  casa::CountedPtr<Scantable> dosigref( const casa::CountedPtr<Scantable>& sig,
+                                        const casa::CountedPtr<Scantable>& ref,
+                                        int smoothref=1,
+                                        casa::Float tsysv=0.0,
+                                        casa::Float tau=0.0 );
+
+ /**
+    * Calibrate GBT Nod scan pairs (translated from GBTIDL)
+    * @param s Scantable which contains Nod scans
+    * @param scans Vector of scan numbers
+    * @param smoothref optional Boxcar smooth width of the reference scans
+    * @param tsysv optional scalar Tsys value at the zenith, required to
+    * set tau, as well
+    * @param tau optional scalar Tau value
+    * @param tcal optional scalar Tcal, CAL temperature (K)
+    * @return casa::CountedPtr<Scantable> which holds calibrated scans
+    */
+  casa::CountedPtr<Scantable> donod( const casa::CountedPtr<Scantable>& s,
+                                     const std::vector<int>& scans,
+                                     int smoothref=1,
+                                     casa::Float tsysv=0.0,
+                                     casa::Float tau=0.0,
+                                     casa::Float tcal=0.0 );
+
+  /**
+    * Calibrate frequency switched scans (translated from GBTIDL)
+    * @param s Scantable which contains frequency switched  scans
+    * @param scans Vector of scan numbers
+    * @param smoothref optional Boxcar smooth width of the reference scans
+    * @param tsysv optional scalar Tsys value at the zenith, required to
+    * set tau, as well
+    * @param tau optional scalar Tau value
+    * @param tcal optional scalar Tcal, CAL temperature (K)
+    * @return casa::CountedPtr<Scantable> which holds calibrated scans
+    */
+  casa::CountedPtr<Scantable> dofs( const casa::CountedPtr<Scantable>& s,
+                                    const std::vector<int>& scans,
+                                    int smoothref=1,
+                                    casa::Float tsysv=0.0,
+                                    casa::Float tau=0.0,
+                                    casa::Float tcal=0.0 );
+
+
   casa::CountedPtr<Scantable>
     freqSwitch( const casa::CountedPtr<Scantable>& in );

trunk/src/STMathWrapper.h

@@ -91 +91 @@
                                                preserve ) ); }
 
+  ScantableWrapper dototalpower( const ScantableWrapper& calon,
+                             const ScantableWrapper& caloff, casa::Float tcal= 0 )
+  { return ScantableWrapper( STMath::dototalpower( calon.getCP(), caloff.getCP(), tcal ) ); }
+
+  ScantableWrapper dosigref( const ScantableWrapper& sig,
+                             const ScantableWrapper& ref,
+                             int smoothref = 0, casa::Float tsysv=0.0, casa::Float tau=0.0)
+  { return ScantableWrapper( STMath::dosigref( sig.getCP(), ref.getCP(), smoothref, tsysv, tau ) ); }
+
+  ScantableWrapper donod( const ScantableWrapper& s,
+                          const std::vector<int>& scans,
+                          int smoothref = 0,
+                          casa::Float tsysv=0.0, casa::Float tau=0.0, casa::Float tcal=0.0 )
+  { return ScantableWrapper( STMath::donod( s.getCP(), scans, smoothref, tsysv, tau, tcal ) ); }
+
+  ScantableWrapper dofs( const ScantableWrapper& s,
+                         const std::vector<int>& scans,
+                         int smoothref = 0,
+                         casa::Float tsysv=0.0, casa::Float tau=0.0, casa::Float tcal=0.0 )
+  { return ScantableWrapper( STMath::dofs( s.getCP(), scans, smoothref, tsysv, tau, tcal ) ); }
+
   ScantableWrapper
     freqSwitch( const ScantableWrapper& in )

trunk/src/STSelector.cpp

@@ -83 +83 @@
 void asap::STSelector::setName( const std::string& sname )
 {
-  std::string sql = "SELECT from $1 WHERE SRCNAME == pattern('"+sname+"')";
+  std::string sql = "SELECT FROM $1 WHERE SRCNAME == pattern('"+sname+"')";
   setTaQL(sql);
 }

trunk/src/STTcal.cpp

@@ -68 +68 @@
 }
 
+/*** rewrite this for handling of GBT data 
 uInt STTcal::addEntry( const String& time, const Vector<Float>& cal)
 {
@@ -90 +91 @@
   return resultid;
 }
+***/
+
+uInt STTcal::addEntry( const String& time, const Vector<Float>& cal)
+{
+  // test if this already exists
+  // TT - different Tcal values for each polarization, feed, and
+  // data description. So there may be multiple entries for the same
+  // time stamp.
+  uInt resultid;
+  uInt rno = table_.nrow();
+  //table_.addRow();
+  // get last assigned tcal_id and increment
+  if ( rno == 0 ) {
+    resultid = 0;
+  }
+  else {
+    idCol_.get(rno-1, resultid);
+    resultid++;
+  }
+  table_.addRow();
+  tcalCol_.put(rno, cal);
+  timeCol_.put(rno, time);
+  idCol_.put(rno, resultid);
+  return resultid;
+}
+
 
 void STTcal::getEntry( String& time, Vector<Float>& tcal, uInt id )

trunk/src/STWriter.cpp

@@ -116 +116 @@
 
   // Extract the header from the table.
+  // this is a little different from what I have done
+  // before. Need to check with the Offline User Test data
   STHeader hdr = in->getHeader();
   //const Int nPol  = hdr.npol;
@@ -123 +125 @@
   Vector<uInt> nPol(nIF),nChan(nIF);
   Vector<Bool> havexpol(nIF);
+  String fluxUnit = hdr.fluxunit;
+
   nPol = 0;nChan = 0; havexpol = False;
   for (uint i=0;i<ifs.size();++i) {
@@ -264 +268 @@
   pushLog(String(oss));
   writer_->close();
-
+  //if MS2 delete POINTING table exists and copy the one in the keyword
+  if ( format_ == "MS2" ) {
+    replacePtTab(table, filename);
+  }
   return 0;
 }
@@ -314 +321 @@
 }
 
-
-}
+// For writing MS data, if there is the reference to
+// original pointing table it replace it by it.
+void STWriter::replacePtTab (const Table& tab, const std::string& fname)
+{
+  String oldPtTabName = fname;
+  oldPtTabName.append("/POINTING");
+  if ( tab.keywordSet().isDefined("POINTING") ) {
+    String PointingTab = tab.keywordSet().asString("POINTING");
+    if ( Table::isReadable(PointingTab) ) {
+      Table newPtTab(PointingTab, Table::Old);
+      newPtTab.copy(oldPtTabName, Table::New);
+      ostringstream oss;
+      oss << "STWriter: copied  " <<PointingTab  << " to " << fname;
+      pushLog(String(oss));
+    }
+  }
+}
+
+}

trunk/src/STWriter.h

@@ -81 +81 @@
                       casa::Vector<casa::Complex>& xpol,
                       const casa::Table& tab);
+
+  void replacePtTab(const casa::Table& tab, const std::string& fname);
+
   std::string     format_;
   PKSwriter* writer_;

trunk/src/Scantable.cpp

@@ -1028 +1028 @@
 }
 
+std::string asap::Scantable::getAntennaName() const
+{
+  String out;
+  table_.keywordSet().get("AntennaName", out);
+  return out;
+}
+
+int asap::Scantable::checkScanInfo(const std::vector<int>& scanlist) const
+{
+  String tbpath;
+  int ret = 0;
+  if ( table_.keywordSet().isDefined("GBT_GO") ) {
+    table_.keywordSet().get("GBT_GO", tbpath);
+    Table t(tbpath,Table::Old);
+    // check each scan if other scan of the pair exist
+    int nscan = scanlist.size();
+    for (int i = 0; i < nscan; i++) {
+      Table subt = t( t.col("SCAN") == scanlist[i]+1 );
+      if (subt.nrow()==0) {
+        cerr <<"Scan "<<scanlist[i]<<" cannot be found in the scantable."<<endl;
+        ret = 1;
+        break;
+      }
+      ROTableRow row(subt);
+      const TableRecord& rec = row.get(0);
+      int scan1seqn = rec.asuInt("PROCSEQN");
+      int laston1 = rec.asuInt("LASTON");
+      if ( rec.asuInt("PROCSIZE")==2 ) {
+        if ( i < nscan-1 ) {
+          Table subt2 = t( t.col("SCAN") == scanlist[i+1]+1 );
+          if ( subt2.nrow() == 0) {
+            cerr<<"Scan "<<scanlist[i+1]<<" cannot be found in the scantable."<<endl;
+            ret = 1;
+            break;
+          }
+          ROTableRow row2(subt2);
+          const TableRecord& rec2 = row2.get(0);
+          int scan2seqn = rec2.asuInt("PROCSEQN");
+          int laston2 = rec2.asuInt("LASTON");
+          if (scan1seqn == 1 && scan2seqn == 2) {
+            if (laston1 == laston2) {
+              cerr<<"A valid scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
+              i +=1;
+            }
+            else {
+              cerr<<"Incorrect scan pair ["<<scanlist[i]<<","<<scanlist[i+1]<<"]"<<endl;
+            }
+          }
+          else if (scan1seqn==2 && scan2seqn == 1) {
+            if (laston1 == laston2) {
+              cerr<<"["<<scanlist[i]<<","<<scanlist[i+1]<<"] is a valid scan pair but in incorrect order."<<endl;
+              ret = 1;
+              break;
+            }
+          }
+          else {
+            cerr<<"The other scan for  "<<scanlist[i]<<" appears to be missing. Check the input scan numbers."<<endl;
+            ret = 1;
+            break;
+          }
+        }
+      }
+      else {
+        cerr<<"The scan does not appear to be standard obsevation."<<endl;
+      }
+    //if ( i >= nscan ) break;
+    }
+  }
+  else {
+    cerr<<"No reference to GBT_GO table."<<endl;
+    ret = 1;
+  }
+  return ret;
+}
+
+std::vector<double>  asap::Scantable::getDirectionVector(int whichrow) const
+{
+  Vector<Double> Dir = dirCol_(whichrow).getAngle("rad").getValue();
+  std::vector<double> dir;
+  Dir.tovector(dir);
+  return dir;
+}
+
 }
  //namespace asap

trunk/src/Scantable.h

@@ -381 +381 @@
     { STFitEntry fe; fitTable_.getEntry(fe, mfitidCol_(row)); return fe; }
 
+  //Added by TT
+  /**
+   * Get the antenna name
+   * @return antenna name string
+   */
+  std::string getAntennaName() const;
+
+  /**
+   * For GBT MS data only. check a scan list
+   * against the information found in GBT_GO table for
+   * scan number orders to get correct pairs.
+   * @param[in] scan list
+   * @return status
+   */
+  int checkScanInfo(const std::vector<int>& scanlist) const;
+
+  /**
+   * Get the direction as a vector, for a specific row
+   * @param[in] whichrow the row numbyyer
+   * @return the direction in a vector
+   */
+  std::vector<double> getDirectionVector(int whichrow) const;
+
 private:
 

trunk/src/ScantableWrapper.h

@@ -199 +199 @@
     { return table_->columnNames(); }
 
+  std::string getAntennaName() const
+    { return table_->getAntennaName(); }
+
+  int checkScanInfo(const vector<int>& scanlist) const
+    { return table_->checkScanInfo(scanlist); }
+
+  std::vector<double>  getDirectionVector(int whichrow) const
+    { return table_->getDirectionVector(whichrow); }
+
 private:
   casa::CountedPtr<Scantable> table_;

trunk/src/python_Fitter.cpp

@@ -55 +55 @@
         .def("reset", &Fitter::reset)
         .def("fit", &Fitter::fit)
+        .def("lfit", &Fitter::lfit)
         .def("evaluate", &Fitter::evaluate)
       ;

trunk/src/python_STMath.cpp

@@ -52 +52 @@
         .def("_auto_quotient", &STMathWrapper::autoQuotient)
         .def("_quotient", &STMathWrapper::quotient)
+        .def("_dototalpower", &STMathWrapper::dototalpower)
+        .def("_dosigref", &STMathWrapper::dosigref)
+        .def("_donod", &STMathWrapper::donod)
+        .def("_dofs", &STMathWrapper::dofs)
         .def("_stats", &STMathWrapper::statistic)
         .def("_freqswitch", &STMathWrapper::freqSwitch)

trunk/src/python_Scantable.cpp

@@ -102 +102 @@
     .def("_getdirection", &ScantableWrapper::getDirectionString,
          (boost::python::arg("whichrow")=0) )
+    .def("get_antennaname", &ScantableWrapper::getAntennaName)
     .def("_flag", &ScantableWrapper::flag)
     .def("_save",  &ScantableWrapper::makePersistent)
@@ -121 +122 @@
     .def("_recalcazel", &ScantableWrapper::calculateAZEL)
     .def("_setsourcetype", &ScantableWrapper::setSourceType)
+    .def("_getdirectionvec", &ScantableWrapper::getDirectionVector)
   ;
 };