source: branches/hpc33/src/STMath.cpp @ 2460

Last change on this file since 2460 was 2460, checked in by KohjiNakamura, 12 years ago

probably this change get average() faster.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 175.8 KB
Line 
1//
2// C++ Implementation: STMath
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12
13#include <sstream>
14
15#include <casa/iomanip.h>
16#include <casa/Arrays/MaskArrLogi.h>
17#include <casa/Arrays/MaskArrMath.h>
18#include <casa/Arrays/ArrayLogical.h>
19#include <casa/Arrays/ArrayMath.h>
20#include <casa/Arrays/Slice.h>
21#include <casa/Arrays/Slicer.h>
22#include <casa/BasicSL/String.h>
23#include <casa/Containers/Block.h>
24#include <casa/Containers/RecordField.h>
25#include <casa/Exceptions/Error.h>
26#include <casa/Logging/LogIO.h>
27
28#include <coordinates/Coordinates/CoordinateSystem.h>
29#include <coordinates/Coordinates/CoordinateUtil.h>
30#include <coordinates/Coordinates/FrequencyAligner.h>
31#include <coordinates/Coordinates/SpectralCoordinate.h>
32
33#include <lattices/Lattices/LatticeUtilities.h>
34
35#include <scimath/Functionals/Polynomial.h>
36#include <scimath/Mathematics/Convolver.h>
37#include <scimath/Mathematics/VectorKernel.h>
38
39#include <tables/Tables/ExprNode.h>
40#include <tables/Tables/ReadAsciiTable.h>
41#include <tables/Tables/TableCopy.h>
42#include <tables/Tables/TableIter.h>
43#include <tables/Tables/TableParse.h>
44#include <tables/Tables/TableRecord.h>
45#include <tables/Tables/TableRow.h>
46#include <tables/Tables/TableVector.h>
47#include <tables/Tables/TabVecMath.h>
48
49#include <atnf/PKSIO/SrcType.h>
50
51#include "RowAccumulator.h"
52#include "STAttr.h"
53#include "STMath.h"
54#include "STSelector.h"
55#include "Accelerator.h"
56
57using namespace casa;
58using namespace asap;
59
60
61// 2012/02/17 TN
62// Since STGrid is implemented, average doesn't consider direction
63// when accumulating
64// tolerance for direction comparison (rad)
65// #define TOL_OTF    1.0e-15
66// #define TOL_POINT  2.9088821e-4  // 1 arcmin
67
68STMath::STMath(bool insitu) :
69  insitu_(insitu)
70{
71}
72
73
74STMath::~STMath()
75{
76}
77
78CountedPtr<Scantable>
79STMath::average( const std::vector<CountedPtr<Scantable> >& in,
80                 const std::vector<bool>& mask,
81                 const std::string& weight,
82                 const std::string& avmode)
83{
84  LogIO os( LogOrigin( "STMath", "average()", WHERE ) ) ;
85  if ( avmode == "SCAN" && in.size() != 1 )
86    throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
87                    "Use merge first."));
88  WeightType wtype = stringToWeight(weight);
89
90  // 2012/02/17 TN
91  // Since STGrid is implemented, average doesn't consider direction
92  // when accumulating
93  // check if OTF observation
94//   String obstype = in[0]->getHeader().obstype ;
95//   Double tol = 0.0 ;
96//   if ( (obstype.find( "OTF" ) != String::npos) || (obstype.find( "OBSERVE_TARGET" ) != String::npos) ) {
97//     tol = TOL_OTF ;
98//   }
99//   else {
100//     tol = TOL_POINT ;
101//   }
102
103  // output
104  // clone as this is non insitu
105  bool insitu = insitu_;
106  setInsitu(false);
107  CountedPtr< Scantable > out = getScantable(in[0], true);
108  setInsitu(insitu);
109  std::vector<CountedPtr<Scantable> >::const_iterator stit = in.begin();
110  ++stit;
111  while ( stit != in.end() ) {
112    out->appendToHistoryTable((*stit)->history());
113    ++stit;
114  }
115
116  Table& tout = out->table();
117
118  /// @todo check if all scantables are conformant
119
120  ArrayColumn<Float> specColOut(tout,"SPECTRA");
121  ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
122  ArrayColumn<Float> tsysColOut(tout,"TSYS");
123  ScalarColumn<Double> mjdColOut(tout,"TIME");
124  ScalarColumn<Double> intColOut(tout,"INTERVAL");
125  ScalarColumn<uInt> cycColOut(tout,"CYCLENO");
126  ScalarColumn<uInt> scanColOut(tout,"SCANNO");
127
128  // set up the output table rows. These are based on the structure of the
129  // FIRST scantable in the vector
130  const Table& baset = in[0]->table();
131
132  Block<String> cols(3);
133  cols[0] = String("BEAMNO");
134  cols[1] = String("IFNO");
135  cols[2] = String("POLNO");
136  if ( avmode == "SOURCE" ) {
137    cols.resize(4);
138    cols[3] = String("SRCNAME");
139  }
140  if ( avmode == "SCAN"  && in.size() == 1) {
141    //cols.resize(4);
142    //cols[3] = String("SCANNO");
143    cols.resize(5);
144    cols[3] = String("SRCNAME");
145    cols[4] = String("SCANNO");
146  }
147  uInt outrowCount = 0;
148  TableIterator iter(baset, cols);
149//   int count = 0 ;
150  while (!iter.pastEnd()) {
151    Table subt = iter.table();
152    // copy the first row of this selection into the new table
153    tout.addRow();
154    TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
155    // re-index to 0
156    if ( avmode != "SCAN" && avmode != "SOURCE" ) {
157      scanColOut.put(outrowCount, uInt(0));
158    }
159    ++outrowCount;
160    // 2012/02/17 TN
161    // Since STGrid is implemented, average doesn't consider direction
162    // when accumulating
163//     MDirection::ScalarColumn dircol ;
164//     dircol.attach( subt, "DIRECTION" ) ;
165//     Int length = subt.nrow() ;
166//     vector< Vector<Double> > dirs ;
167//     vector<int> indexes ;
168//     for ( Int i = 0 ; i < length ; i++ ) {
169//       Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
170//       //os << << count++ << ": " ;
171//       //os << "[" << t[0] << "," << t[1] << "]" << LogIO::POST ;
172//       bool adddir = true ;
173//       for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
174//         //if ( allTrue( t == dirs[j] ) ) {
175//         Double dx = t[0] - dirs[j][0] ;
176//         Double dy = t[1] - dirs[j][1] ;
177//         Double dd = sqrt( dx * dx + dy * dy ) ;
178//         //if ( allNearAbs( t, dirs[j], tol ) ) {
179//         if ( dd <= tol ) {
180//           adddir = false ;
181//           break ;
182//         }
183//       }
184//       if ( adddir ) {
185//         dirs.push_back( t ) ;
186//         indexes.push_back( i ) ;
187//       }
188//     }
189//     uInt rowNum = dirs.size() ;
190//     tout.addRow( rowNum ) ;
191//     for ( uInt i = 0 ; i < rowNum ; i++ ) {
192//       TableCopy::copyRows( tout, subt, outrowCount+i, indexes[i], 1 ) ;
193//       // re-index to 0
194//       if ( avmode != "SCAN" && avmode != "SOURCE" ) {
195//         scanColOut.put(outrowCount+i, uInt(0));
196//       }       
197//     }
198//     outrowCount += rowNum ;
199    ++iter;
200  }
201  RowAccumulator acc(wtype);
202  Vector<Bool> cmask(mask);
203  acc.setUserMask(cmask);
204  ROTableRow row(tout);
205  ROArrayColumn<Float> specCol, tsysCol;
206  ROArrayColumn<uChar> flagCol;
207  ROScalarColumn<Double> mjdCol, intCol;
208  ROScalarColumn<Int> scanIDCol;
209
210  Vector<uInt> rowstodelete;
211
212  for (uInt i=0; i < tout.nrow(); ++i) {
213    for ( int j=0; j < int(in.size()); ++j ) {
214      const Table& tin = in[j]->table();
215      const TableRecord& rec = row.get(i);
216      ROScalarColumn<Double> tmp(tin, "TIME");
217      Double td;tmp.get(0,td);
218
219#if 1
220      static char const*const colNames1[] = { "IFNO", "BEAMNO", "POLNO" };
221      uInt const values1[] = { rec.asuInt("IFNO"), rec.asuInt("BEAMNO"), rec.asuInt("POLNO") };
222      SingleTypeEqPredicate<uInt, 3> myPred(tin, colNames1, values1);
223      CustomTableExprNode myExpr(tin, myPred);
224      Table basesubt = tin(myExpr);
225#else
226      Table basesubt = tin( tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
227                         && tin.col("IFNO") == Int(rec.asuInt("IFNO"))
228                         && tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
229#endif
230      Table subt;
231      if ( avmode == "SOURCE") {
232        subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME"));
233      } else if (avmode == "SCAN") {
234        subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME")
235                      && basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
236      } else {
237        subt = basesubt;
238      }
239
240      // 2012/02/17 TN
241      // Since STGrid is implemented, average doesn't consider direction
242      // when accumulating
243//       vector<uInt> removeRows ;
244//       uInt nrsubt = subt.nrow() ;
245//       for ( uInt irow = 0 ; irow < nrsubt ; irow++ ) {
246//         //if ( !allTrue((subt.col("DIRECTION").getArrayDouble(TableExprId(irow)))==rec.asArrayDouble("DIRECTION")) ) {
247//         Vector<Double> x0 = (subt.col("DIRECTION").getArrayDouble(TableExprId(irow))) ;
248//         Vector<Double> x1 = rec.asArrayDouble("DIRECTION") ;
249//         double dx = x0[0] - x1[0];
250//         double dy = x0[1] - x1[1];
251//         Double dd = sqrt( dx * dx + dy * dy ) ;
252//         //if ( !allNearAbs((subt.col("DIRECTION").getArrayDouble(TableExprId(irow))), rec.asArrayDouble("DIRECTION"), tol ) ) {
253//         if ( dd > tol ) {
254//           removeRows.push_back( irow ) ;
255//         }
256//       }
257//       if ( removeRows.size() != 0 ) {
258//         subt.removeRow( removeRows ) ;
259//       }
260     
261//       if ( nrsubt == removeRows.size() )
262//         throw(AipsError("Averaging data is empty.")) ;
263
264      specCol.attach(subt,"SPECTRA");
265      flagCol.attach(subt,"FLAGTRA");
266      tsysCol.attach(subt,"TSYS");
267      intCol.attach(subt,"INTERVAL");
268      mjdCol.attach(subt,"TIME");
269      Vector<Float> spec,tsys;
270      Vector<uChar> flag;
271      Double inter,time;
272      for (uInt k = 0; k < subt.nrow(); ++k ) {
273        flagCol.get(k, flag);
274        Vector<Bool> bflag(flag.shape());
275        convertArray(bflag, flag);
276        /*
277        if ( allEQ(bflag, True) ) {
278        continue;//don't accumulate
279        }
280        */
281        specCol.get(k, spec);
282        tsysCol.get(k, tsys);
283        intCol.get(k, inter);
284        mjdCol.get(k, time);
285        // spectrum has to be added last to enable weighting by the other values
286        acc.add(spec, !bflag, tsys, inter, time);
287      }
288
289      // If there exists a channel at which all the input spectra are masked,
290      // spec has 'nan' values for that channel and it may affect the following
291      // processes. To avoid this, replacing 'nan' values in spec with
292      // weighted-mean of all spectra in the following line.
293      // (done for CAS-2776, 2011/04/07 by Wataru Kawasaki)
294      acc.replaceNaN();
295    }
296    const Vector<Bool>& msk = acc.getMask();
297    if ( allEQ(msk, False) ) {
298      uint n = rowstodelete.nelements();
299      rowstodelete.resize(n+1, True);
300      rowstodelete[n] = i;
301      continue;
302    }
303    //write out
304    if (acc.state()) {
305      Vector<uChar> flg(msk.shape());
306      convertArray(flg, !msk);
307      for (uInt k = 0; k < flg.nelements(); ++k) {
308        uChar userFlag = 1 << 7;
309        if (msk[k]==True) userFlag = 0 << 7;
310        flg(k) = userFlag;
311      }
312
313      flagColOut.put(i, flg);
314      specColOut.put(i, acc.getSpectrum());
315      tsysColOut.put(i, acc.getTsys());
316      intColOut.put(i, acc.getInterval());
317      mjdColOut.put(i, acc.getTime());
318      // we should only have one cycle now -> reset it to be 0
319      // frequency switched data has different CYCLENO for different IFNO
320      // which requires resetting this value
321      cycColOut.put(i, uInt(0));
322    } else {
323      ostringstream oss;
324      oss << "For output row="<<i<<", all input rows of data are flagged. no averaging" << endl;
325      pushLog(String(oss));
326    }
327    acc.reset();
328  }
329
330  if (rowstodelete.nelements() > 0) {
331    os << rowstodelete << LogIO::POST ;
332    tout.removeRow(rowstodelete);
333    if (tout.nrow() == 0) {
334      throw(AipsError("Can't average fully flagged data."));
335    }
336  }
337  return out;
338}
339
340CountedPtr< Scantable >
341STMath::averageChannel( const CountedPtr < Scantable > & in,
342                          const std::string & mode,
343                          const std::string& avmode )
344{
345  (void) mode; // currently unused
346  // 2012/02/17 TN
347  // Since STGrid is implemented, average doesn't consider direction
348  // when accumulating
349  // check if OTF observation
350//   String obstype = in->getHeader().obstype ;
351//   Double tol = 0.0 ;
352//   if ( obstype.find( "OTF" ) != String::npos ) {
353//     tol = TOL_OTF ;
354//   }
355//   else {
356//     tol = TOL_POINT ;
357//   }
358
359  // clone as this is non insitu
360  bool insitu = insitu_;
361  setInsitu(false);
362  CountedPtr< Scantable > out = getScantable(in, true);
363  setInsitu(insitu);
364  Table& tout = out->table();
365  ArrayColumn<Float> specColOut(tout,"SPECTRA");
366  ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
367  ArrayColumn<Float> tsysColOut(tout,"TSYS");
368  ScalarColumn<uInt> scanColOut(tout,"SCANNO");
369  ScalarColumn<Double> intColOut(tout, "INTERVAL");
370  Table tmp = in->table().sort("BEAMNO");
371  Block<String> cols(3);
372  cols[0] = String("BEAMNO");
373  cols[1] = String("IFNO");
374  cols[2] = String("POLNO");
375  if ( avmode == "SCAN") {
376    cols.resize(4);
377    cols[3] = String("SCANNO");
378  }
379  uInt outrowCount = 0;
380  uChar userflag = 1 << 7;
381  TableIterator iter(tmp, cols);
382  while (!iter.pastEnd()) {
383    Table subt = iter.table();
384    ROArrayColumn<Float> specCol, tsysCol;
385    ROArrayColumn<uChar> flagCol;
386    ROScalarColumn<Double> intCol(subt, "INTERVAL");
387    specCol.attach(subt,"SPECTRA");
388    flagCol.attach(subt,"FLAGTRA");
389    tsysCol.attach(subt,"TSYS");
390
391    tout.addRow();
392    TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
393    if ( avmode != "SCAN") {
394      scanColOut.put(outrowCount, uInt(0));
395    }
396    Vector<Float> tmp;
397    specCol.get(0, tmp);
398    uInt nchan = tmp.nelements();
399    // have to do channel by channel here as MaskedArrMath
400    // doesn't have partialMedians
401    Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
402    Vector<Float> outspec(nchan);
403    Vector<uChar> outflag(nchan,0);
404    Vector<Float> outtsys(1);/// @fixme when tsys is channel based
405    for (uInt i=0; i<nchan; ++i) {
406      Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
407      MaskedArray<Float> ma = maskedArray(specs,flags);
408      outspec[i] = median(ma);
409      if ( allEQ(ma.getMask(), False) )
410        outflag[i] = userflag;// flag data
411    }
412    outtsys[0] = median(tsysCol.getColumn());
413    specColOut.put(outrowCount, outspec);
414    flagColOut.put(outrowCount, outflag);
415    tsysColOut.put(outrowCount, outtsys);
416    Double intsum = sum(intCol.getColumn());
417    intColOut.put(outrowCount, intsum);
418    ++outrowCount;
419    ++iter;
420
421    // 2012/02/17 TN
422    // Since STGrid is implemented, average doesn't consider direction
423    // when accumulating
424//     MDirection::ScalarColumn dircol ;
425//     dircol.attach( subt, "DIRECTION" ) ;
426//     Int length = subt.nrow() ;
427//     vector< Vector<Double> > dirs ;
428//     vector<int> indexes ;
429//     // Handle MX mode averaging
430//     if (in->nbeam() > 1 ) {     
431//       length = 1;
432//     }
433//     for ( Int i = 0 ; i < length ; i++ ) {
434//       Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
435//       bool adddir = true ;
436//       for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
437//         //if ( allTrue( t == dirs[j] ) ) {
438//         Double dx = t[0] - dirs[j][0] ;
439//         Double dy = t[1] - dirs[j][1] ;
440//         Double dd = sqrt( dx * dx + dy * dy ) ;
441//         //if ( allNearAbs( t, dirs[j], tol ) ) {
442//         if ( dd <= tol ) {
443//           adddir = false ;
444//           break ;
445//         }
446//       }
447//       if ( adddir ) {
448//         dirs.push_back( t ) ;
449//         indexes.push_back( i ) ;
450//       }
451//     }
452//     uInt rowNum = dirs.size() ;
453//     tout.addRow( rowNum );
454//     for ( uInt i = 0 ; i < rowNum ; i++ ) {
455//       TableCopy::copyRows(tout, subt, outrowCount+i, indexes[i], 1) ;
456//       // Handle MX mode averaging
457//       if ( avmode != "SCAN") {
458//         scanColOut.put(outrowCount+i, uInt(0));
459//       }
460//     }
461//     MDirection::ScalarColumn dircolOut ;
462//     dircolOut.attach( tout, "DIRECTION" ) ;
463//     for ( uInt irow = 0 ; irow < rowNum ; irow++ ) {
464//       Vector<Double> t = \
465//      dircolOut(outrowCount+irow).getAngle(Unit(String("rad"))).getValue() ;
466//       Vector<Float> tmp;
467//       specCol.get(0, tmp);
468//       uInt nchan = tmp.nelements();
469//       // have to do channel by channel here as MaskedArrMath
470//       // doesn't have partialMedians
471//       Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
472//       // mask spectra for different DIRECTION
473//       for ( uInt jrow = 0 ; jrow < subt.nrow() ; jrow++ ) {
474//         Vector<Double> direction = \
475//        dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
476//         //if ( t[0] != direction[0] || t[1] != direction[1] ) {
477//         Double dx = t[0] - direction[0];
478//         Double dy = t[1] - direction[1];
479//         Double dd = sqrt(dx*dx + dy*dy);
480//         //if ( !allNearAbs( t, direction, tol ) ) {
481//         if ( dd > tol &&  in->nbeam() < 2 ) {
482//           flags[jrow] = userflag ;
483//         }
484//       }
485//       Vector<Float> outspec(nchan);
486//       Vector<uChar> outflag(nchan,0);
487//       Vector<Float> outtsys(1);/// @fixme when tsys is channel based
488//       for (uInt i=0; i<nchan; ++i) {
489//         Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
490//         MaskedArray<Float> ma = maskedArray(specs,flags);
491//         outspec[i] = median(ma);
492//         if ( allEQ(ma.getMask(), False) )
493//           outflag[i] = userflag;// flag data
494//       }
495//       outtsys[0] = median(tsysCol.getColumn());
496//       specColOut.put(outrowCount+irow, outspec);
497//       flagColOut.put(outrowCount+irow, outflag);
498//       tsysColOut.put(outrowCount+irow, outtsys);
499//       Vector<Double> integ = intCol.getColumn() ;
500//       MaskedArray<Double> mi = maskedArray( integ, flags ) ;
501//       Double intsum = sum(mi);
502//       intColOut.put(outrowCount+irow, intsum);
503//     }
504//     outrowCount += rowNum ;
505//     ++iter;
506  }
507  return out;
508}
509
510CountedPtr< Scantable > STMath::getScantable(const CountedPtr< Scantable >& in,
511                                             bool droprows)
512{
513  if (insitu_) {
514    return in;
515  }
516  else {
517    // clone
518    return CountedPtr<Scantable>(new Scantable(*in, Bool(droprows)));
519  }
520}
521
522CountedPtr< Scantable > STMath::unaryOperate( const CountedPtr< Scantable >& in,
523                                              float val,
524                                              const std::string& mode,
525                                              bool tsys )
526{
527  CountedPtr< Scantable > out = getScantable(in, false);
528  Table& tab = out->table();
529  ArrayColumn<Float> specCol(tab,"SPECTRA");
530  ArrayColumn<Float> tsysCol(tab,"TSYS");
531  if (mode=="DIV") val = 1.0/val ;
532  else if (mode=="SUB") val *= -1.0 ;
533  for (uInt i=0; i<tab.nrow(); ++i) {
534    Vector<Float> spec;
535    Vector<Float> ts;
536    specCol.get(i, spec);
537    tsysCol.get(i, ts);
538    if (mode == "MUL" || mode == "DIV") {
539      //if (mode == "DIV") val = 1.0/val;
540      spec *= val;
541      specCol.put(i, spec);
542      if ( tsys ) {
543        ts *= val;
544        tsysCol.put(i, ts);
545      }
546    } else if ( mode == "ADD"  || mode == "SUB") {
547      //if (mode == "SUB") val *= -1.0;
548      spec += val;
549      specCol.put(i, spec);
550      if ( tsys ) {
551        ts += val;
552        tsysCol.put(i, ts);
553      }
554    }
555  }
556  return out;
557}
558
559CountedPtr< Scantable > STMath::arrayOperate( const CountedPtr< Scantable >& in,
560                                              const std::vector<float> val,
561                                              const std::string& mode,
562                                              const std::string& opmode,
563                                              bool tsys )
564{
565  CountedPtr< Scantable > out ;
566  if ( opmode == "channel" ) {
567    out = arrayOperateChannel( in, val, mode, tsys ) ;
568  }
569  else if ( opmode == "row" ) {
570    out = arrayOperateRow( in, val, mode, tsys ) ;
571  }
572  else {
573    throw( AipsError( "Unknown array operation mode." ) ) ;
574  }
575  return out ;
576}
577
578CountedPtr< Scantable > STMath::arrayOperateChannel( const CountedPtr< Scantable >& in,
579                                                     const std::vector<float> val,
580                                                     const std::string& mode,
581                                                     bool tsys )
582{
583  if ( val.size() == 1 ){
584    return unaryOperate( in, val[0], mode, tsys ) ;
585  }
586
587  // conformity of SPECTRA and TSYS
588  if ( tsys ) {
589    TableIterator titer(in->table(), "IFNO");
590    while ( !titer.pastEnd() ) {
591      ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ;
592      ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ;
593      Array<Float> spec = specCol.getColumn() ;
594      Array<Float> ts = tsysCol.getColumn() ;
595      if ( !spec.conform( ts ) ) {
596        throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ;
597      }
598      titer.next() ;
599    }
600  }
601
602  // check if all spectra in the scantable have the same number of channel
603  vector<uInt> nchans;
604  vector<uInt> ifnos = in->getIFNos() ;
605  for ( uInt i = 0 ; i < ifnos.size() ; i++ ) {
606    nchans.push_back( in->nchan( ifnos[i] ) ) ;
607  }
608  Vector<uInt> mchans( nchans ) ;
609  if ( anyNE( mchans, mchans[0] ) ) {
610    throw( AipsError("All spectra in the input scantable must have the same number of channel for vector operation." ) ) ;
611  }
612
613  // check if vector size is equal to nchan
614  Vector<Float> fact( val ) ;
615  if ( fact.nelements() != mchans[0] ) {
616    throw( AipsError("Vector size must be 1 or be same as number of channel.") ) ;
617  }
618
619  // check divided by zero
620  if ( ( mode == "DIV" ) && anyEQ( fact, (float)0.0 ) ) {
621    throw( AipsError("Divided by zero is not recommended." ) ) ;
622  }
623
624  CountedPtr< Scantable > out = getScantable(in, false);
625  Table& tab = out->table();
626  ArrayColumn<Float> specCol(tab,"SPECTRA");
627  ArrayColumn<Float> tsysCol(tab,"TSYS");
628  if (mode == "DIV") fact = (float)1.0 / fact;
629  else if (mode == "SUB") fact *= (float)-1.0 ;
630  for (uInt i=0; i<tab.nrow(); ++i) {
631    Vector<Float> spec;
632    Vector<Float> ts;
633    specCol.get(i, spec);
634    tsysCol.get(i, ts);
635    if (mode == "MUL" || mode == "DIV") {
636      //if (mode == "DIV") fact = (float)1.0 / fact;
637      spec *= fact;
638      specCol.put(i, spec);
639      if ( tsys ) {
640        ts *= fact;
641        tsysCol.put(i, ts);
642      }
643    } else if ( mode == "ADD"  || mode == "SUB") {
644      //if (mode == "SUB") fact *= (float)-1.0 ;
645      spec += fact;
646      specCol.put(i, spec);
647      if ( tsys ) {
648        ts += fact;
649        tsysCol.put(i, ts);
650      }
651    }
652  }
653  return out;
654}
655
656CountedPtr< Scantable > STMath::arrayOperateRow( const CountedPtr< Scantable >& in,
657                                                 const std::vector<float> val,
658                                                 const std::string& mode,
659                                                 bool tsys )
660{
661  if ( val.size() == 1 ) {
662    return unaryOperate( in, val[0], mode, tsys ) ;
663  }
664
665  // conformity of SPECTRA and TSYS
666  if ( tsys ) {
667    TableIterator titer(in->table(), "IFNO");
668    while ( !titer.pastEnd() ) {
669      ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ;
670      ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ;
671      Array<Float> spec = specCol.getColumn() ;
672      Array<Float> ts = tsysCol.getColumn() ;
673      if ( !spec.conform( ts ) ) {
674        throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ;
675      }
676      titer.next() ;
677    }
678  }
679
680  // check if vector size is equal to nrow
681  Vector<Float> fact( val ) ;
682  if (fact.nelements() != uInt(in->nrow())) {
683    throw( AipsError("Vector size must be 1 or be same as number of row.") ) ;
684  }
685
686  // check divided by zero
687  if ( ( mode == "DIV" ) && anyEQ( fact, (float)0.0 ) ) {
688    throw( AipsError("Divided by zero is not recommended." ) ) ;
689  }
690
691  CountedPtr< Scantable > out = getScantable(in, false);
692  Table& tab = out->table();
693  ArrayColumn<Float> specCol(tab,"SPECTRA");
694  ArrayColumn<Float> tsysCol(tab,"TSYS");
695  if (mode == "DIV") fact = (float)1.0 / fact;
696  if (mode == "SUB") fact *= (float)-1.0 ;
697  for (uInt i=0; i<tab.nrow(); ++i) {
698    Vector<Float> spec;
699    Vector<Float> ts;
700    specCol.get(i, spec);
701    tsysCol.get(i, ts);
702    if (mode == "MUL" || mode == "DIV") {
703      spec *= fact[i];
704      specCol.put(i, spec);
705      if ( tsys ) {
706        ts *= fact[i];
707        tsysCol.put(i, ts);
708      }
709    } else if ( mode == "ADD"  || mode == "SUB") {
710      spec += fact[i];
711      specCol.put(i, spec);
712      if ( tsys ) {
713        ts += fact[i];
714        tsysCol.put(i, ts);
715      }
716    }
717  }
718  return out;
719}
720
721CountedPtr< Scantable > STMath::array2dOperate( const CountedPtr< Scantable >& in,
722                                                const std::vector< std::vector<float> > val,
723                                                const std::string& mode,
724                                                bool tsys )
725{
726  // conformity of SPECTRA and TSYS
727  if ( tsys ) {
728    TableIterator titer(in->table(), "IFNO");
729    while ( !titer.pastEnd() ) {
730      ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ;
731      ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ;
732      Array<Float> spec = specCol.getColumn() ;
733      Array<Float> ts = tsysCol.getColumn() ;
734      if ( !spec.conform( ts ) ) {
735        throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ;
736      }
737      titer.next() ;
738    }
739  }
740
741  // some checks
742  vector<uInt> nchans;
743  for (Int i = 0 ; i < in->nrow() ; i++) {
744    nchans.push_back((in->getSpectrum(i)).size());
745  }
746  //Vector<uInt> mchans( nchans ) ;
747  vector< Vector<Float> > facts ;
748  for ( uInt i = 0 ; i < nchans.size() ; i++ ) {
749    Vector<Float> tmp( val[i] ) ;
750    // check divided by zero
751    if ( ( mode == "DIV" ) && anyEQ( tmp, (float)0.0 ) ) {
752      throw( AipsError("Divided by zero is not recommended." ) ) ;
753    }
754    // conformity check
755    if ( tmp.nelements() != nchans[i] ) {
756      stringstream ss ;
757      ss << "Row " << i << ": Vector size must be same as number of channel." ;
758      throw( AipsError( ss.str() ) ) ;
759    }
760    facts.push_back( tmp ) ;
761  }
762
763
764  CountedPtr< Scantable > out = getScantable(in, false);
765  Table& tab = out->table();
766  ArrayColumn<Float> specCol(tab,"SPECTRA");
767  ArrayColumn<Float> tsysCol(tab,"TSYS");
768  for (uInt i=0; i<tab.nrow(); ++i) {
769    Vector<Float> fact = facts[i] ;
770    Vector<Float> spec;
771    Vector<Float> ts;
772    specCol.get(i, spec);
773    tsysCol.get(i, ts);
774    if (mode == "MUL" || mode == "DIV") {
775      if (mode == "DIV") fact = (float)1.0 / fact;
776      spec *= fact;
777      specCol.put(i, spec);
778      if ( tsys ) {
779        ts *= fact;
780        tsysCol.put(i, ts);
781      }
782    } else if ( mode == "ADD"  || mode == "SUB") {
783      if (mode == "SUB") fact *= (float)-1.0 ;
784      spec += fact;
785      specCol.put(i, spec);
786      if ( tsys ) {
787        ts += fact;
788        tsysCol.put(i, ts);
789      }
790    }
791  }
792  return out;
793}
794
795CountedPtr<Scantable> STMath::binaryOperate(const CountedPtr<Scantable>& left,
796                                            const CountedPtr<Scantable>& right,
797                                            const std::string& mode)
798{
799  bool insitu = insitu_;
800  if ( ! left->conformant(*right) ) {
801    throw(AipsError("'left' and 'right' scantables are not conformant."));
802  }
803  setInsitu(false);
804  CountedPtr< Scantable > out = getScantable(left, false);
805  setInsitu(insitu);
806  Table& tout = out->table();
807  Block<String> coln(5);
808  coln[0] = "SCANNO";  coln[1] = "CYCLENO";  coln[2] = "BEAMNO";
809  coln[3] = "IFNO";  coln[4] = "POLNO";
810  Table tmpl = tout.sort(coln);
811  Table tmpr = right->table().sort(coln);
812  ArrayColumn<Float> lspecCol(tmpl,"SPECTRA");
813  ROArrayColumn<Float> rspecCol(tmpr,"SPECTRA");
814  ArrayColumn<uChar> lflagCol(tmpl,"FLAGTRA");
815  ROArrayColumn<uChar> rflagCol(tmpr,"FLAGTRA");
816
817  for (uInt i=0; i<tout.nrow(); ++i) {
818    Vector<Float> lspecvec, rspecvec;
819    Vector<uChar> lflagvec, rflagvec;
820    lspecvec = lspecCol(i);    rspecvec = rspecCol(i);
821    lflagvec = lflagCol(i);    rflagvec = rflagCol(i);
822    MaskedArray<Float> mleft = maskedArray(lspecvec, lflagvec);
823    MaskedArray<Float> mright = maskedArray(rspecvec, rflagvec);
824    if (mode == "ADD") {
825      mleft += mright;
826    } else if ( mode == "SUB") {
827      mleft -= mright;
828    } else if ( mode == "MUL") {
829      mleft *= mright;
830    } else if ( mode == "DIV") {
831      mleft /= mright;
832    } else {
833      throw(AipsError("Illegal binary operator"));
834    }
835    lspecCol.put(i, mleft.getArray());
836  }
837  return out;
838}
839
840
841
842MaskedArray<Float> STMath::maskedArray( const Vector<Float>& s,
843                                        const Vector<uChar>& f)
844{
845  Vector<Bool> mask;
846  mask.resize(f.shape());
847  convertArray(mask, f);
848  return MaskedArray<Float>(s,!mask);
849}
850
851MaskedArray<Double> STMath::maskedArray( const Vector<Double>& s,
852                                         const Vector<uChar>& f)
853{
854  Vector<Bool> mask;
855  mask.resize(f.shape());
856  convertArray(mask, f);
857  return MaskedArray<Double>(s,!mask);
858}
859
860Vector<uChar> STMath::flagsFromMA(const MaskedArray<Float>& ma)
861{
862  const Vector<Bool>& m = ma.getMask();
863  Vector<uChar> flags(m.shape());
864  convertArray(flags, !m);
865  return flags;
866}
867
868CountedPtr< Scantable > STMath::autoQuotient( const CountedPtr< Scantable >& in,
869                                              const std::string & mode,
870                                              bool preserve )
871{
872  /// @todo make other modes available
873  /// modes should be "nearest", "pair"
874  // make this operation non insitu
875  (void) mode; //currently unused
876  const Table& tin = in->table();
877  Table ons = tin(tin.col("SRCTYPE") == Int(SrcType::PSON));
878  Table offs = tin(tin.col("SRCTYPE") == Int(SrcType::PSOFF));
879  if ( offs.nrow() == 0 )
880    throw(AipsError("No 'off' scans present."));
881  // put all "on" scans into output table
882
883  bool insitu = insitu_;
884  setInsitu(false);
885  CountedPtr< Scantable > out = getScantable(in, true);
886  setInsitu(insitu);
887  Table& tout = out->table();
888
889  TableCopy::copyRows(tout, ons);
890  TableRow row(tout);
891  ROScalarColumn<Double> offtimeCol(offs, "TIME");
892  ArrayColumn<Float> outspecCol(tout, "SPECTRA");
893  ROArrayColumn<Float> outtsysCol(tout, "TSYS");
894  ArrayColumn<uChar> outflagCol(tout, "FLAGTRA");
895  for (uInt i=0; i < tout.nrow(); ++i) {
896    const TableRecord& rec = row.get(i);
897    Double ontime = rec.asDouble("TIME");
898    Table presel = offs(offs.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
899                        && offs.col("IFNO") == Int(rec.asuInt("IFNO"))
900                        && offs.col("POLNO") == Int(rec.asuInt("POLNO")) );
901    ROScalarColumn<Double> offtimeCol(presel, "TIME");
902
903    Double mindeltat = min(abs(offtimeCol.getColumn() - ontime));
904    // Timestamp may vary within a cycle ???!!!
905    // increase this by 0.01 sec in case of rounding errors...
906    // There might be a better way to do this.
907    // fix to this fix. TIME is MJD, so 1.0d not 1.0s
908    mindeltat += 0.01/24./60./60.;
909    Table sel = presel( abs(presel.col("TIME")-ontime) <= mindeltat);
910
911    if ( sel.nrow() < 1 )  {
912      throw(AipsError("No closest in time found... This could be a rounding "
913                      "issue. Try quotient instead."));
914    }
915    TableRow offrow(sel);
916    const TableRecord& offrec = offrow.get(0);//should only be one row
917    RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
918    RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
919    RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
920    /// @fixme this assumes tsys is a scalar not vector
921    Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
922    Vector<Float> specon, tsyson;
923    outtsysCol.get(i, tsyson);
924    outspecCol.get(i, specon);
925    Vector<uChar> flagon;
926    outflagCol.get(i, flagon);
927    MaskedArray<Float> mon = maskedArray(specon, flagon);
928    MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
929    MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
930    if (preserve) {
931      quot -= tsysoffscalar;
932    } else {
933      quot -= tsyson[0];
934    }
935    outspecCol.put(i, quot.getArray());
936    outflagCol.put(i, flagsFromMA(quot));
937  }
938  // renumber scanno
939  TableIterator it(tout, "SCANNO");
940  uInt i = 0;
941  while ( !it.pastEnd() ) {
942    Table t = it.table();
943    TableVector<uInt> vec(t, "SCANNO");
944    vec = i;
945    ++i;
946    ++it;
947  }
948  return out;
949}
950
951
952CountedPtr< Scantable > STMath::quotient( const CountedPtr< Scantable > & on,
953                                          const CountedPtr< Scantable > & off,
954                                          bool preserve )
955{
956  bool insitu = insitu_;
957  if ( ! on->conformant(*off) ) {
958    throw(AipsError("'on' and 'off' scantables are not conformant."));
959  }
960  setInsitu(false);
961  CountedPtr< Scantable > out = getScantable(on, false);
962  setInsitu(insitu);
963  Table& tout = out->table();
964  const Table& toff = off->table();
965  TableIterator sit(tout, "SCANNO");
966  TableIterator s2it(toff, "SCANNO");
967  while ( !sit.pastEnd() ) {
968    Table ton = sit.table();
969    TableRow row(ton);
970    Table t = s2it.table();
971    ArrayColumn<Float> outspecCol(ton, "SPECTRA");
972    ROArrayColumn<Float> outtsysCol(ton, "TSYS");
973    ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
974    for (uInt i=0; i < ton.nrow(); ++i) {
975      const TableRecord& rec = row.get(i);
976      Table offsel = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
977                          && t.col("IFNO") == Int(rec.asuInt("IFNO"))
978                          && t.col("POLNO") == Int(rec.asuInt("POLNO")) );
979      if ( offsel.nrow() == 0 )
980        throw AipsError("STMath::quotient: no matching off");
981      TableRow offrow(offsel);
982      const TableRecord& offrec = offrow.get(0);//should be ncycles - take first
983      RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
984      RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
985      RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
986      Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
987      Vector<Float> specon, tsyson;
988      outtsysCol.get(i, tsyson);
989      outspecCol.get(i, specon);
990      Vector<uChar> flagon;
991      outflagCol.get(i, flagon);
992      MaskedArray<Float> mon = maskedArray(specon, flagon);
993      MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
994      MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
995      if (preserve) {
996        quot -= tsysoffscalar;
997      } else {
998        quot -= tsyson[0];
999      }
1000      outspecCol.put(i, quot.getArray());
1001      outflagCol.put(i, flagsFromMA(quot));
1002    }
1003    ++sit;
1004    ++s2it;
1005    // take the first off for each on scan which doesn't have a
1006    // matching off scan
1007    // non <= noff:  matching pairs, non > noff matching pairs then first off
1008    if ( s2it.pastEnd() ) s2it.reset();
1009  }
1010  return out;
1011}
1012
1013// dototalpower (migration of GBTIDL procedure dototalpower.pro)
1014// calibrate the CAL on-off pair. It calculate Tsys and average CAL on-off subintegrations
1015// do it for each cycles in a specific scan.
1016CountedPtr< Scantable > STMath::dototalpower( const CountedPtr< Scantable >& calon,
1017                                              const CountedPtr< Scantable >& caloff, Float tcal )
1018{
1019  if ( ! calon->conformant(*caloff) ) {
1020    throw(AipsError("'CAL on' and 'CAL off' scantables are not conformant."));
1021  }
1022  setInsitu(false);
1023  CountedPtr< Scantable > out = getScantable(caloff, false);
1024  Table& tout = out->table();
1025  const Table& tcon = calon->table();
1026  Vector<Float> tcalout;
1027
1028  std::map<uInt,uInt> tcalIdToRecNoMap;
1029  const Table& calOffTcalTable = caloff->tcal().table();
1030  {
1031    ROScalarColumn<uInt> calOffTcalTable_IDcol(calOffTcalTable, "ID");
1032    const Vector<uInt> tcalIds(calOffTcalTable_IDcol.getColumn());
1033    size_t tcalIdsEnd = tcalIds.nelements();
1034    for (uInt i = 0; i < tcalIdsEnd; i++) {
1035      tcalIdToRecNoMap[tcalIds[i]] = i;
1036    }
1037  }
1038  ROArrayColumn<Float> calOffTcalTable_TCALcol(calOffTcalTable, "TCAL");
1039
1040  if ( tout.nrow() != tcon.nrow() ) {
1041    throw(AipsError("Mismatch in number of rows to form cal on - off pair."));
1042  }
1043  // iteration by scanno or cycle no.
1044  TableIterator sit(tout, "SCANNO");
1045  TableIterator s2it(tcon, "SCANNO");
1046  while ( !sit.pastEnd() ) {
1047    Table toff = sit.table();
1048    TableRow row(toff);
1049    Table t = s2it.table();
1050    ScalarColumn<Double> outintCol(toff, "INTERVAL");
1051    ArrayColumn<Float> outspecCol(toff, "SPECTRA");
1052    ArrayColumn<Float> outtsysCol(toff, "TSYS");
1053    ArrayColumn<uChar> outflagCol(toff, "FLAGTRA");
1054    ROScalarColumn<uInt> outtcalIdCol(toff, "TCAL_ID");
1055    ROScalarColumn<uInt> outpolCol(toff, "POLNO");
1056    ROScalarColumn<Double> onintCol(t, "INTERVAL");
1057    ROArrayColumn<Float> onspecCol(t, "SPECTRA");
1058    ROArrayColumn<Float> ontsysCol(t, "TSYS");
1059    ROArrayColumn<uChar> onflagCol(t, "FLAGTRA");
1060    //ROScalarColumn<uInt> ontcalIdCol(t, "TCAL_ID");
1061
1062    for (uInt i=0; i < toff.nrow(); ++i) {
1063      //skip these checks -> assumes the data order are the same between the cal on off pairs
1064      //
1065      Vector<Float> specCalon, specCaloff;
1066      // to store scalar (mean) tsys
1067      Vector<Float> tsysout(1);
1068      uInt tcalId, polno;
1069      Double offint, onint;
1070      outpolCol.get(i, polno);
1071      outspecCol.get(i, specCaloff);
1072      onspecCol.get(i, specCalon);
1073      Vector<uChar> flagCaloff, flagCalon;
1074      outflagCol.get(i, flagCaloff);
1075      onflagCol.get(i, flagCalon);
1076      outtcalIdCol.get(i, tcalId);
1077      outintCol.get(i, offint);
1078      onintCol.get(i, onint);
1079      // caluculate mean Tsys
1080      uInt nchan = specCaloff.nelements();
1081      // percentage of edge cut off
1082      uInt pc = 10;
1083      uInt bchan = nchan/pc;
1084      uInt echan = nchan-bchan;
1085
1086      Slicer chansl(IPosition(1,bchan-1), IPosition(1,echan-1), IPosition(1,1),Slicer::endIsLast);
1087      Vector<Float> testsubsp = specCaloff(chansl);
1088      MaskedArray<Float> spoff = maskedArray( specCaloff(chansl),flagCaloff(chansl) );
1089      MaskedArray<Float> spon = maskedArray( specCalon(chansl),flagCalon(chansl) );
1090      MaskedArray<Float> spdiff = spon-spoff;
1091      uInt noff = spoff.nelementsValid();
1092      //uInt non = spon.nelementsValid();
1093      uInt ndiff = spdiff.nelementsValid();
1094      Float meantsys;
1095
1096/**
1097      Double subspec, subdiff;
1098      uInt usednchan;
1099      subspec = 0;
1100      subdiff = 0;
1101      usednchan = 0;
1102      for(uInt k=(bchan-1); k<echan; k++) {
1103        subspec += specCaloff[k];
1104        subdiff += static_cast<Double>(specCalon[k]-specCaloff[k]);
1105        ++usednchan;
1106      }
1107**/
1108      // get tcal if input tcal <= 0
1109      Float tcalUsed;
1110      tcalUsed = tcal;
1111      if ( tcal <= 0.0 ) {
1112        uInt tcalRecNo = tcalIdToRecNoMap[tcalId];
1113        calOffTcalTable_TCALcol.get(tcalRecNo, tcalout);
1114//         if (polno<=3) {
1115//           tcalUsed = tcalout[polno];
1116//         }
1117//         else {
1118//           tcalUsed = tcalout[0];
1119//         }
1120        if ( tcalout.size() == 1 )
1121          tcalUsed = tcalout[0] ;
1122        else if ( tcalout.size() == nchan )
1123          tcalUsed = mean(tcalout) ;
1124        else {
1125          uInt ipol = polno ;
1126          if ( ipol > 3 ) ipol = 0 ;
1127          tcalUsed = tcalout[ipol] ;
1128        }
1129      }
1130
1131      Float meanoff;
1132      Float meandiff;
1133      if (noff && ndiff) {
1134         //Debug
1135         //if(noff!=ndiff) cerr<<"noff and ndiff is not equal"<<endl;
1136         //LogIO os( LogOrigin( "STMath", "dototalpower()", WHERE ) ) ;
1137         //if(noff!=ndiff) os<<"noff and ndiff is not equal"<<LogIO::POST;
1138         meanoff = sum(spoff)/noff;
1139         meandiff = sum(spdiff)/ndiff;
1140         meantsys= (meanoff/meandiff )*tcalUsed + tcalUsed/2;
1141      }
1142      else {
1143         meantsys=1;
1144      }
1145
1146      tsysout[0] = Float(meantsys);
1147      MaskedArray<Float> mcaloff = maskedArray(specCaloff, flagCaloff);
1148      MaskedArray<Float> mcalon = maskedArray(specCalon, flagCalon);
1149      MaskedArray<Float> sig =   Float(0.5) * (mcaloff + mcalon);
1150      //uInt ncaloff = mcaloff.nelementsValid();
1151      //uInt ncalon = mcalon.nelementsValid();
1152
1153      outintCol.put(i, offint+onint);
1154      outspecCol.put(i, sig.getArray());
1155      outflagCol.put(i, flagsFromMA(sig));
1156      outtsysCol.put(i, tsysout);
1157    }
1158    ++sit;
1159    ++s2it;
1160  }
1161  return out;
1162}
1163
1164//dosigref - migrated from GBT IDL's dosigref.pro, do calibration of position switch
1165// observatiions.
1166// input: sig and ref scantables, and an optional boxcar smoothing width(default width=0,
1167//        no smoothing).
1168// output: resultant scantable [= (sig-ref/ref)*tsys]
1169CountedPtr< Scantable > STMath::dosigref( const CountedPtr < Scantable >& sig,
1170                                          const CountedPtr < Scantable >& ref,
1171                                          int smoothref,
1172                                          casa::Float tsysv,
1173                                          casa::Float tau )
1174{
1175if ( ! ref->conformant(*sig) ) {
1176    throw(AipsError("'sig' and 'ref' scantables are not conformant."));
1177  }
1178  setInsitu(false);
1179  CountedPtr< Scantable > out = getScantable(sig, false);
1180  CountedPtr< Scantable > smref;
1181  if ( smoothref > 1 ) {
1182    float fsmoothref = static_cast<float>(smoothref);
1183    std::string inkernel = "boxcar";
1184    smref = smooth(ref, inkernel, fsmoothref );
1185    ostringstream oss;
1186    oss<<"Applied smoothing of "<<fsmoothref<<" on the reference."<<endl;
1187    pushLog(String(oss));
1188  }
1189  else {
1190    smref = ref;
1191  }
1192  Table& tout = out->table();
1193  const Table& tref = smref->table();
1194  if ( tout.nrow() != tref.nrow() ) {
1195    throw(AipsError("Mismatch in number of rows to form on-source and reference pair."));
1196  }
1197  // iteration by scanno? or cycle no.
1198  TableIterator sit(tout, "SCANNO");
1199  TableIterator s2it(tref, "SCANNO");
1200  while ( !sit.pastEnd() ) {
1201    Table ton = sit.table();
1202    Table t = s2it.table();
1203    ScalarColumn<Double> outintCol(ton, "INTERVAL");
1204    ArrayColumn<Float> outspecCol(ton, "SPECTRA");
1205    ArrayColumn<Float> outtsysCol(ton, "TSYS");
1206    ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
1207    ArrayColumn<Float> refspecCol(t, "SPECTRA");
1208    ROScalarColumn<Double> refintCol(t, "INTERVAL");
1209    ROArrayColumn<Float> reftsysCol(t, "TSYS");
1210    ArrayColumn<uChar> refflagCol(t, "FLAGTRA");
1211    ROScalarColumn<Float> refelevCol(t, "ELEVATION");
1212    for (uInt i=0; i < ton.nrow(); ++i) {
1213
1214      Double onint, refint;
1215      Vector<Float> specon, specref;
1216      // to store scalar (mean) tsys
1217      Vector<Float> tsysref;
1218      outintCol.get(i, onint);
1219      refintCol.get(i, refint);
1220      outspecCol.get(i, specon);
1221      refspecCol.get(i, specref);
1222      Vector<uChar> flagref, flagon;
1223      outflagCol.get(i, flagon);
1224      refflagCol.get(i, flagref);
1225      reftsysCol.get(i, tsysref);
1226
1227      Float tsysrefscalar;
1228      if ( tsysv > 0.0 ) {
1229        ostringstream oss;
1230        Float elev;
1231        refelevCol.get(i, elev);
1232        oss << "user specified Tsys = " << tsysv;
1233        // do recalc elevation if EL = 0
1234        if ( elev == 0 ) {
1235          throw(AipsError("EL=0, elevation data is missing."));
1236        } else {
1237          if ( tau <= 0.0 ) {
1238            throw(AipsError("Valid tau is not supplied."));
1239          } else {
1240            tsysrefscalar = tsysv * exp(tau/elev);
1241          }
1242        }
1243        oss << ", corrected (for El) tsys= "<<tsysrefscalar;
1244        pushLog(String(oss));
1245      }
1246      else {
1247        tsysrefscalar = tsysref[0];
1248      }
1249      //get quotient spectrum
1250      MaskedArray<Float> mref = maskedArray(specref, flagref);
1251      MaskedArray<Float> mon = maskedArray(specon, flagon);
1252      MaskedArray<Float> specres =   tsysrefscalar*((mon - mref)/mref);
1253      Double resint = onint*refint*smoothref/(onint+refint*smoothref);
1254
1255      //Debug
1256      //cerr<<"Tsys used="<<tsysrefscalar<<endl;
1257      //LogIO os( LogOrigin( "STMath", "dosigref", WHERE ) ) ;
1258      //os<<"Tsys used="<<tsysrefscalar<<LogIO::POST;
1259      // fill the result, replay signal tsys by reference tsys
1260      outintCol.put(i, resint);
1261      outspecCol.put(i, specres.getArray());
1262      outflagCol.put(i, flagsFromMA(specres));
1263      outtsysCol.put(i, tsysref);
1264    }
1265    ++sit;
1266    ++s2it;
1267  }
1268  return out;
1269}
1270
1271CountedPtr< Scantable > STMath::donod(const casa::CountedPtr<Scantable>& s,
1272                                     const std::vector<int>& scans,
1273                                     int smoothref,
1274                                     casa::Float tsysv,
1275                                     casa::Float tau,
1276                                     casa::Float tcal )
1277
1278{
1279  setInsitu(false);
1280  STSelector sel;
1281  std::vector<int> scan1, scan2, beams, types;
1282  std::vector< vector<int> > scanpair;
1283  //std::vector<string> calstate;
1284  std::vector<int> calstate;
1285  String msg;
1286
1287  CountedPtr< Scantable > s1b1on, s1b1off, s1b2on, s1b2off;
1288  CountedPtr< Scantable > s2b1on, s2b1off, s2b2on, s2b2off;
1289
1290  std::vector< CountedPtr< Scantable > > sctables;
1291  sctables.push_back(s1b1on);
1292  sctables.push_back(s1b1off);
1293  sctables.push_back(s1b2on);
1294  sctables.push_back(s1b2off);
1295  sctables.push_back(s2b1on);
1296  sctables.push_back(s2b1off);
1297  sctables.push_back(s2b2on);
1298  sctables.push_back(s2b2off);
1299
1300  //check scanlist
1301  int n=s->checkScanInfo(scans);
1302  if (n==1) {
1303     throw(AipsError("Incorrect scan pairs. "));
1304  }
1305
1306  // Assume scans contain only a pair of consecutive scan numbers.
1307  // It is assumed that first beam, b1,  is on target.
1308  // There is no check if the first beam is on or not.
1309  if ( scans.size()==1 ) {
1310    scan1.push_back(scans[0]);
1311    scan2.push_back(scans[0]+1);
1312  } else if ( scans.size()==2 ) {
1313   scan1.push_back(scans[0]);
1314   scan2.push_back(scans[1]);
1315  } else {
1316    if ( scans.size()%2 == 0 ) {
1317      for (uInt i=0; i<scans.size(); i++) {
1318        if (i%2 == 0) {
1319          scan1.push_back(scans[i]);
1320        }
1321        else {
1322          scan2.push_back(scans[i]);
1323        }
1324      }
1325    } else {
1326      throw(AipsError("Odd numbers of scans, cannot form pairs."));
1327    }
1328  }
1329  scanpair.push_back(scan1);
1330  scanpair.push_back(scan2);
1331  //calstate.push_back("*calon");
1332  //calstate.push_back("*[^calon]");
1333  calstate.push_back(SrcType::NODCAL);
1334  calstate.push_back(SrcType::NOD);
1335  CountedPtr< Scantable > ws = getScantable(s, false);
1336  uInt l=0;
1337  while ( l < sctables.size() ) {
1338    for (uInt i=0; i < 2; i++) {
1339      for (uInt j=0; j < 2; j++) {
1340        for (uInt k=0; k < 2; k++) {
1341          sel.reset();
1342          sel.setScans(scanpair[i]);
1343          //sel.setName(calstate[k]);
1344          types.clear();
1345          types.push_back(calstate[k]);
1346          sel.setTypes(types);
1347          beams.clear();
1348          beams.push_back(j);
1349          sel.setBeams(beams);
1350          ws->setSelection(sel);
1351          sctables[l]= getScantable(ws, false);
1352          l++;
1353        }
1354      }
1355    }
1356  }
1357
1358  // replace here by splitData or getData functionality
1359  CountedPtr< Scantable > sig1;
1360  CountedPtr< Scantable > ref1;
1361  CountedPtr< Scantable > sig2;
1362  CountedPtr< Scantable > ref2;
1363  CountedPtr< Scantable > calb1;
1364  CountedPtr< Scantable > calb2;
1365
1366  msg=String("Processing dototalpower for subset of the data");
1367  ostringstream oss1;
1368  oss1 << msg  << endl;
1369  pushLog(String(oss1));
1370  // Debug for IRC CS data
1371  //float tcal1=7.0;
1372  //float tcal2=4.0;
1373  sig1 = dototalpower(sctables[0], sctables[1], tcal=tcal);
1374  ref1 = dototalpower(sctables[2], sctables[3], tcal=tcal);
1375  ref2 = dototalpower(sctables[4], sctables[5], tcal=tcal);
1376  sig2 = dototalpower(sctables[6], sctables[7], tcal=tcal);
1377
1378  // correction of user-specified tsys for elevation here
1379
1380  // dosigref calibration
1381  msg=String("Processing dosigref for subset of the data");
1382  ostringstream oss2;
1383  oss2 << msg  << endl;
1384  pushLog(String(oss2));
1385  calb1=dosigref(sig1,ref2,smoothref,tsysv,tau);
1386  calb2=dosigref(sig2,ref1,smoothref,tsysv,tau);
1387
1388  // iteration by scanno or cycle no.
1389  Table& tcalb1 = calb1->table();
1390  Table& tcalb2 = calb2->table();
1391  TableIterator sit(tcalb1, "SCANNO");
1392  TableIterator s2it(tcalb2, "SCANNO");
1393  while ( !sit.pastEnd() ) {
1394    Table t1 = sit.table();
1395    Table t2= s2it.table();
1396    ArrayColumn<Float> outspecCol(t1, "SPECTRA");
1397    ArrayColumn<Float> outtsysCol(t1, "TSYS");
1398    ArrayColumn<uChar> outflagCol(t1, "FLAGTRA");
1399    ScalarColumn<Double> outintCol(t1, "INTERVAL");
1400    ArrayColumn<Float> t2specCol(t2, "SPECTRA");
1401    ROArrayColumn<Float> t2tsysCol(t2, "TSYS");
1402    ArrayColumn<uChar> t2flagCol(t2, "FLAGTRA");
1403    ROScalarColumn<Double> t2intCol(t2, "INTERVAL");
1404    for (uInt i=0; i < t1.nrow(); ++i) {
1405      Vector<Float> spec1, spec2;
1406      // to store scalar (mean) tsys
1407      Vector<Float> tsys1, tsys2;
1408      Vector<uChar> flag1, flag2;
1409      Double tint1, tint2;
1410      outspecCol.get(i, spec1);
1411      t2specCol.get(i, spec2);
1412      outflagCol.get(i, flag1);
1413      t2flagCol.get(i, flag2);
1414      outtsysCol.get(i, tsys1);
1415      t2tsysCol.get(i, tsys2);
1416      outintCol.get(i, tint1);
1417      t2intCol.get(i, tint2);
1418      // average
1419      // assume scalar tsys for weights
1420      Float wt1, wt2, tsyssq1, tsyssq2;
1421      tsyssq1 = tsys1[0]*tsys1[0];
1422      tsyssq2 = tsys2[0]*tsys2[0];
1423      wt1 = Float(tint1)/tsyssq1;
1424      wt2 = Float(tint2)/tsyssq2;
1425      Float invsumwt=1/(wt1+wt2);
1426      MaskedArray<Float> mspec1 = maskedArray(spec1, flag1);
1427      MaskedArray<Float> mspec2 = maskedArray(spec2, flag2);
1428      MaskedArray<Float> avspec =  invsumwt * (wt1*mspec1 + wt2*mspec2);
1429      //Array<Float> avtsys =  Float(0.5) * (tsys1 + tsys2);
1430      // cerr<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<endl;
1431      // LogIO os( LogOrigin( "STMath", "donod", WHERE ) ) ;
1432      // os<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<LogIO::POST;
1433      tsys1[0] = sqrt(tsyssq1 + tsyssq2);
1434      Array<Float> avtsys =  tsys1;
1435
1436      outspecCol.put(i, avspec.getArray());
1437      outflagCol.put(i, flagsFromMA(avspec));
1438      outtsysCol.put(i, avtsys);
1439    }
1440    ++sit;
1441    ++s2it;
1442  }
1443  return calb1;
1444}
1445
1446//GBTIDL version of frequency switched data calibration
1447CountedPtr< Scantable > STMath::dofs( const CountedPtr< Scantable >& s,
1448                                      const std::vector<int>& scans,
1449                                      int smoothref,
1450                                      casa::Float tsysv,
1451                                      casa::Float tau,
1452                                      casa::Float tcal )
1453{
1454
1455
1456  (void) scans; //currently unused
1457  STSelector sel;
1458  CountedPtr< Scantable > ws = getScantable(s, false);
1459  CountedPtr< Scantable > sig, sigwcal, ref, refwcal;
1460  CountedPtr< Scantable > calsig, calref, out, out1, out2;
1461  Bool nofold=False;
1462  vector<int> types ;
1463
1464  //split the data
1465  //sel.setName("*_fs");
1466  types.push_back( SrcType::FSON ) ;
1467  sel.setTypes( types ) ;
1468  ws->setSelection(sel);
1469  sig = getScantable(ws,false);
1470  sel.reset();
1471  types.clear() ;
1472  //sel.setName("*_fs_calon");
1473  types.push_back( SrcType::FONCAL ) ;
1474  sel.setTypes( types ) ;
1475  ws->setSelection(sel);
1476  sigwcal = getScantable(ws,false);
1477  sel.reset();
1478  types.clear() ;
1479  //sel.setName("*_fsr");
1480  types.push_back( SrcType::FSOFF ) ;
1481  sel.setTypes( types ) ;
1482  ws->setSelection(sel);
1483  ref = getScantable(ws,false);
1484  sel.reset();
1485  types.clear() ;
1486  //sel.setName("*_fsr_calon");
1487  types.push_back( SrcType::FOFFCAL ) ;
1488  sel.setTypes( types ) ;
1489  ws->setSelection(sel);
1490  refwcal = getScantable(ws,false);
1491  sel.reset() ;
1492  types.clear() ;
1493
1494  calsig = dototalpower(sigwcal, sig, tcal=tcal);
1495  calref = dototalpower(refwcal, ref, tcal=tcal);
1496
1497  out1=dosigref(calsig,calref,smoothref,tsysv,tau);
1498  out2=dosigref(calref,calsig,smoothref,tsysv,tau);
1499
1500  Table& tabout1=out1->table();
1501  Table& tabout2=out2->table();
1502  ROScalarColumn<uInt> freqidCol1(tabout1, "FREQ_ID");
1503  ScalarColumn<uInt> freqidCol2(tabout2, "FREQ_ID");
1504  ROArrayColumn<Float> specCol(tabout2, "SPECTRA");
1505  Vector<Float> spec; specCol.get(0, spec);
1506  uInt nchan = spec.nelements();
1507  uInt freqid1; freqidCol1.get(0,freqid1);
1508  uInt freqid2; freqidCol2.get(0,freqid2);
1509  Double rp1, rp2, rv1, rv2, inc1, inc2;
1510  out1->frequencies().getEntry(rp1, rv1, inc1, freqid1);
1511  out2->frequencies().getEntry(rp2, rv2, inc2, freqid2);
1512  //cerr << out1->frequencies().table().nrow() << " " << out2->frequencies().table().nrow() << endl ;
1513  //LogIO os( LogOrigin( "STMath", "dofs()", WHERE ) ) ;
1514  //os << out1->frequencies().table().nrow() << " " << out2->frequencies().table().nrow() << LogIO::POST ;
1515  if (rp1==rp2) {
1516    Double foffset = rv1 - rv2;
1517    uInt choffset = static_cast<uInt>(foffset/abs(inc2));
1518    if (choffset >= nchan) {
1519      //cerr<<"out-band frequency switching, no folding"<<endl;
1520      LogIO os( LogOrigin( "STMath", "dofs()", WHERE ) ) ;
1521      os<<"out-band frequency switching, no folding"<<LogIO::POST;
1522      nofold = True;
1523    }
1524  }
1525
1526  if (nofold) {
1527    std::vector< CountedPtr< Scantable > > tabs;
1528    tabs.push_back(out1);
1529    tabs.push_back(out2);
1530    out = merge(tabs);
1531  }
1532  else {
1533    //out = out1;
1534    Double choffset = ( rv1 - rv2 ) / inc2 ;
1535    out = dofold( out1, out2, choffset ) ;
1536  }
1537   
1538  return out;
1539}
1540
1541CountedPtr<Scantable> STMath::dofold( const CountedPtr<Scantable> &sig,
1542                                      const CountedPtr<Scantable> &ref,
1543                                      Double choffset,
1544                                      Double choffset2 )
1545{
1546  LogIO os( LogOrigin( "STMath", "dofold", WHERE ) ) ;
1547  os << "choffset=" << choffset << " choffset2=" << choffset2 << LogIO::POST ;
1548
1549  // output scantable
1550  CountedPtr<Scantable> out = getScantable( sig, false ) ;
1551
1552  // separate choffset to integer part and decimal part
1553  Int ioffset = (Int)choffset ;
1554  Double doffset = choffset - ioffset ;
1555  Int ioffset2 = (Int)choffset2 ;
1556  Double doffset2 = choffset2 - ioffset2 ;
1557  os << "ioffset=" << ioffset << " doffset=" << doffset << LogIO::POST ;
1558  os << "ioffset2=" << ioffset2 << " doffset2=" << doffset2 << LogIO::POST ; 
1559
1560  // get column
1561  ROArrayColumn<Float> specCol1( sig->table(), "SPECTRA" ) ;
1562  ROArrayColumn<Float> specCol2( ref->table(), "SPECTRA" ) ;
1563  ROArrayColumn<Float> tsysCol1( sig->table(), "TSYS" ) ;
1564  ROArrayColumn<Float> tsysCol2( ref->table(), "TSYS" ) ;
1565  ROArrayColumn<uChar> flagCol1( sig->table(), "FLAGTRA" ) ;
1566  ROArrayColumn<uChar> flagCol2( ref->table(), "FLAGTRA" ) ;
1567  ROScalarColumn<Double> mjdCol1( sig->table(), "TIME" ) ;
1568  ROScalarColumn<Double> mjdCol2( ref->table(), "TIME" ) ;
1569  ROScalarColumn<Double> intervalCol1( sig->table(), "INTERVAL" ) ;
1570  ROScalarColumn<Double> intervalCol2( ref->table(), "INTERVAL" ) ;
1571
1572  // check
1573  if ( ioffset == 0 ) {
1574    LogIO os( LogOrigin( "STMath", "dofold()", WHERE ) ) ;
1575    os << "channel offset is zero, no folding" << LogIO::POST ;
1576    return out ;
1577  }
1578  int nchan = ref->nchan() ;
1579  if ( abs(ioffset) >= nchan ) {
1580    LogIO os( LogOrigin( "STMath", "dofold()", WHERE ) ) ;
1581    os << "out-band frequency switching, no folding" << LogIO::POST ;
1582    return out ;
1583  }
1584
1585  // attach column for output scantable
1586  ArrayColumn<Float> specColOut( out->table(), "SPECTRA" ) ;
1587  ArrayColumn<uChar> flagColOut( out->table(), "FLAGTRA" ) ;
1588  ArrayColumn<Float> tsysColOut( out->table(), "TSYS" ) ;
1589  ScalarColumn<Double> mjdColOut( out->table(), "TIME" ) ;
1590  ScalarColumn<Double> intervalColOut( out->table(), "INTERVAL" ) ;
1591  ScalarColumn<uInt> fidColOut( out->table(), "FREQ_ID" ) ;
1592
1593  // for each row
1594  // assume that the data order are same between sig and ref
1595  RowAccumulator acc( asap::W_TINTSYS ) ;
1596  for ( int i = 0 ; i < sig->nrow() ; i++ ) {
1597    // get values
1598    Vector<Float> spsig ;
1599    specCol1.get( i, spsig ) ;
1600    Vector<Float> spref ;
1601    specCol2.get( i, spref ) ;
1602    Vector<Float> tsyssig ;
1603    tsysCol1.get( i, tsyssig ) ;
1604    Vector<Float> tsysref ;
1605    tsysCol2.get( i, tsysref ) ;
1606    Vector<uChar> flagsig ;
1607    flagCol1.get( i, flagsig ) ;
1608    Vector<uChar> flagref ;
1609    flagCol2.get( i, flagref ) ;
1610    Double timesig ;
1611    mjdCol1.get( i, timesig ) ;
1612    Double timeref ;
1613    mjdCol2.get( i, timeref ) ;
1614    Double intsig ;
1615    intervalCol1.get( i, intsig ) ;
1616    Double intref ;
1617    intervalCol2.get( i, intref ) ;
1618
1619    // shift reference spectra
1620    int refchan = spref.nelements() ;
1621    Vector<Float> sspref( spref.nelements() ) ;
1622    Vector<Float> stsysref( tsysref.nelements() ) ;
1623    Vector<uChar> sflagref( flagref.nelements() ) ;
1624    if ( ioffset > 0 ) {
1625      // SPECTRA and FLAGTRA
1626      for ( int j = 0 ; j < refchan-ioffset ; j++ ) {
1627        sspref[j] = spref[j+ioffset] ;
1628        sflagref[j] = flagref[j+ioffset] ;
1629      }
1630      for ( int j = refchan-ioffset ; j < refchan ; j++ ) {
1631        sspref[j] = spref[j-refchan+ioffset] ;
1632        sflagref[j] = flagref[j-refchan+ioffset] ;
1633      }
1634      spref = sspref.copy() ;
1635      flagref = sflagref.copy() ;
1636      for ( int j = 0 ; j < refchan - 1 ; j++ ) {
1637        sspref[j] = doffset * spref[j+1] + ( 1.0 - doffset ) * spref[j] ;
1638        sflagref[j] = flagref[j+1] + flagref[j] ;
1639      }
1640      sspref[refchan-1] = doffset * spref[0] + ( 1.0 - doffset ) * spref[refchan-1] ;
1641      sflagref[refchan-1] = flagref[0] + flagref[refchan-1] ;
1642
1643      // TSYS
1644      if ( spref.nelements() == tsysref.nelements() ) {
1645        for ( int j = 0 ; j < refchan-ioffset ; j++ ) {
1646          stsysref[j] = tsysref[j+ioffset] ;
1647        }
1648        for ( int j = refchan-ioffset ; j < refchan ; j++ ) {
1649          stsysref[j] = tsysref[j-refchan+ioffset] ;
1650        }
1651        tsysref = stsysref.copy() ;
1652        for ( int j = 0 ; j < refchan - 1 ; j++ ) {
1653          stsysref[j] = doffset * tsysref[j+1] + ( 1.0 - doffset ) * tsysref[j] ;
1654        }
1655        stsysref[refchan-1] = doffset * tsysref[0] + ( 1.0 - doffset ) * tsysref[refchan-1] ;
1656      }
1657    }
1658    else {
1659      // SPECTRA and FLAGTRA
1660      for ( int j = 0 ; j < abs(ioffset) ; j++ ) {
1661        sspref[j] = spref[refchan+ioffset+j] ;
1662        sflagref[j] = flagref[refchan+ioffset+j] ;
1663      }
1664      for ( int j = abs(ioffset) ; j < refchan ; j++ ) {
1665        sspref[j] = spref[j+ioffset] ;
1666        sflagref[j] = flagref[j+ioffset] ;
1667      }
1668      spref = sspref.copy() ;
1669      flagref = sflagref.copy() ;
1670      sspref[0] = doffset * spref[refchan-1] + ( 1.0 - doffset ) * spref[0] ;
1671      sflagref[0] = flagref[0] + flagref[refchan-1] ;
1672      for ( int j = 1 ; j < refchan ; j++ ) {
1673        sspref[j] = doffset * spref[j-1] + ( 1.0 - doffset ) * spref[j] ;
1674        sflagref[j] = flagref[j-1] + flagref[j] ;
1675      }
1676      // TSYS
1677      if ( spref.nelements() == tsysref.nelements() ) {
1678        for ( int j = 0 ; j < abs(ioffset) ; j++ ) {
1679          stsysref[j] = tsysref[refchan+ioffset+j] ;
1680        }
1681        for ( int j = abs(ioffset) ; j < refchan ; j++ ) {
1682          stsysref[j] = tsysref[j+ioffset] ;
1683        }
1684        tsysref = stsysref.copy() ;
1685        stsysref[0] = doffset * tsysref[refchan-1] + ( 1.0 - doffset ) * tsysref[0] ;
1686        for ( int j = 1 ; j < refchan ; j++ ) {
1687          stsysref[j] = doffset * tsysref[j-1] + ( 1.0 - doffset ) * tsysref[j] ;
1688        }
1689      }
1690    }
1691
1692    // shift signal spectra if necessary (only for APEX?)
1693    if ( choffset2 != 0.0 ) {
1694      int sigchan = spsig.nelements() ;
1695      Vector<Float> sspsig( spsig.nelements() ) ;
1696      Vector<Float> stsyssig( tsyssig.nelements() ) ;
1697      Vector<uChar> sflagsig( flagsig.nelements() ) ;
1698      if ( ioffset2 > 0 ) {
1699        // SPECTRA and FLAGTRA
1700        for ( int j = 0 ; j < sigchan-ioffset2 ; j++ ) {
1701          sspsig[j] = spsig[j+ioffset2] ;
1702          sflagsig[j] = flagsig[j+ioffset2] ;
1703        }
1704        for ( int j = sigchan-ioffset2 ; j < sigchan ; j++ ) {
1705          sspsig[j] = spsig[j-sigchan+ioffset2] ;
1706          sflagsig[j] = flagsig[j-sigchan+ioffset2] ;
1707        }
1708        spsig = sspsig.copy() ;
1709        flagsig = sflagsig.copy() ;
1710        for ( int j = 0 ; j < sigchan - 1 ; j++ ) {
1711          sspsig[j] = doffset2 * spsig[j+1] + ( 1.0 - doffset2 ) * spsig[j] ;
1712          sflagsig[j] = flagsig[j+1] || flagsig[j] ;
1713        }
1714        sspsig[sigchan-1] = doffset2 * spsig[0] + ( 1.0 - doffset2 ) * spsig[sigchan-1] ;
1715        sflagsig[sigchan-1] = flagsig[0] || flagsig[sigchan-1] ;
1716        // TSTS
1717        if ( spsig.nelements() == tsyssig.nelements() ) {
1718          for ( int j = 0 ; j < sigchan-ioffset2 ; j++ ) {
1719            stsyssig[j] = tsyssig[j+ioffset2] ;
1720          }
1721          for ( int j = sigchan-ioffset2 ; j < sigchan ; j++ ) {
1722            stsyssig[j] = tsyssig[j-sigchan+ioffset2] ;
1723          }
1724          tsyssig = stsyssig.copy() ;
1725          for ( int j = 0 ; j < sigchan - 1 ; j++ ) {
1726            stsyssig[j] = doffset2 * tsyssig[j+1] + ( 1.0 - doffset2 ) * tsyssig[j] ;
1727          }
1728          stsyssig[sigchan-1] = doffset2 * tsyssig[0] + ( 1.0 - doffset2 ) * tsyssig[sigchan-1] ;
1729        }
1730      }
1731      else {
1732        // SPECTRA and FLAGTRA
1733        for ( int j = 0 ; j < abs(ioffset2) ; j++ ) {
1734          sspsig[j] = spsig[sigchan+ioffset2+j] ;
1735          sflagsig[j] = flagsig[sigchan+ioffset2+j] ;
1736        }
1737        for ( int j = abs(ioffset2) ; j < sigchan ; j++ ) {
1738          sspsig[j] = spsig[j+ioffset2] ;
1739          sflagsig[j] = flagsig[j+ioffset2] ;
1740        }
1741        spsig = sspsig.copy() ;
1742        flagsig = sflagsig.copy() ;
1743        sspsig[0] = doffset2 * spsig[sigchan-1] + ( 1.0 - doffset2 ) * spsig[0] ;
1744        sflagsig[0] = flagsig[0] + flagsig[sigchan-1] ;
1745        for ( int j = 1 ; j < sigchan ; j++ ) {
1746          sspsig[j] = doffset2 * spsig[j-1] + ( 1.0 - doffset2 ) * spsig[j] ;
1747          sflagsig[j] = flagsig[j-1] + flagsig[j] ;
1748        }
1749        // TSYS
1750        if ( spsig.nelements() == tsyssig.nelements() ) {
1751          for ( int j = 0 ; j < abs(ioffset2) ; j++ ) {
1752            stsyssig[j] = tsyssig[sigchan+ioffset2+j] ;
1753          }
1754          for ( int j = abs(ioffset2) ; j < sigchan ; j++ ) {
1755            stsyssig[j] = tsyssig[j+ioffset2] ;
1756          }
1757          tsyssig = stsyssig.copy() ;
1758          stsyssig[0] = doffset2 * tsyssig[sigchan-1] + ( 1.0 - doffset2 ) * tsyssig[0] ;
1759          for ( int j = 1 ; j < sigchan ; j++ ) {
1760            stsyssig[j] = doffset2 * tsyssig[j-1] + ( 1.0 - doffset2 ) * tsyssig[j] ;
1761          }
1762        }
1763      }
1764    }
1765
1766    // folding
1767    acc.add( spsig, !flagsig, tsyssig, intsig, timesig ) ;
1768    acc.add( sspref, !sflagref, stsysref, intref, timeref ) ;
1769   
1770    // put result
1771    specColOut.put( i, acc.getSpectrum() ) ;
1772    const Vector<Bool> &msk = acc.getMask() ;
1773    Vector<uChar> flg( msk.shape() ) ;
1774    convertArray( flg, !msk ) ;
1775    flagColOut.put( i, flg ) ;
1776    tsysColOut.put( i, acc.getTsys() ) ;
1777    intervalColOut.put( i, acc.getInterval() ) ;
1778    mjdColOut.put( i, acc.getTime() ) ;
1779    // change FREQ_ID to unshifted IF setting (only for APEX?)
1780    if ( choffset2 != 0.0 ) {
1781      uInt freqid = fidColOut( 0 ) ; // assume single-IF data
1782      double refpix, refval, increment ;
1783      out->frequencies().getEntry( refpix, refval, increment, freqid ) ;
1784      refval -= choffset * increment ;
1785      uInt newfreqid = out->frequencies().addEntry( refpix, refval, increment ) ;
1786      Vector<uInt> freqids = fidColOut.getColumn() ;
1787      for ( uInt j = 0 ; j < freqids.nelements() ; j++ ) {
1788        if ( freqids[j] == freqid )
1789          freqids[j] = newfreqid ;
1790      }
1791      fidColOut.putColumn( freqids ) ;
1792    }
1793
1794    acc.reset() ;
1795  }
1796
1797  return out ;
1798}
1799
1800
1801CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
1802{
1803  // make copy or reference
1804  CountedPtr< Scantable > out = getScantable(in, false);
1805  Table& tout = out->table();
1806  Block<String> cols(4);
1807  cols[0] = String("SCANNO");
1808  cols[1] = String("CYCLENO");
1809  cols[2] = String("BEAMNO");
1810  cols[3] = String("POLNO");
1811  TableIterator iter(tout, cols);
1812  while (!iter.pastEnd()) {
1813    Table subt = iter.table();
1814    // this should leave us with two rows for the two IFs....if not ignore
1815    if (subt.nrow() != 2 ) {
1816      continue;
1817    }
1818    ArrayColumn<Float> specCol(subt, "SPECTRA");
1819    ArrayColumn<Float> tsysCol(subt, "TSYS");
1820    ArrayColumn<uChar> flagCol(subt, "FLAGTRA");
1821    Vector<Float> onspec,offspec, ontsys, offtsys;
1822    Vector<uChar> onflag, offflag;
1823    tsysCol.get(0, ontsys);   tsysCol.get(1, offtsys);
1824    specCol.get(0, onspec);   specCol.get(1, offspec);
1825    flagCol.get(0, onflag);   flagCol.get(1, offflag);
1826    MaskedArray<Float> on  = maskedArray(onspec, onflag);
1827    MaskedArray<Float> off = maskedArray(offspec, offflag);
1828    MaskedArray<Float> oncopy = on.copy();
1829
1830    on /= off; on -= 1.0f;
1831    on *= ontsys[0];
1832    off /= oncopy; off -= 1.0f;
1833    off *= offtsys[0];
1834    specCol.put(0, on.getArray());
1835    const Vector<Bool>& m0 = on.getMask();
1836    Vector<uChar> flags0(m0.shape());
1837    convertArray(flags0, !m0);
1838    flagCol.put(0, flags0);
1839
1840    specCol.put(1, off.getArray());
1841    const Vector<Bool>& m1 = off.getMask();
1842    Vector<uChar> flags1(m1.shape());
1843    convertArray(flags1, !m1);
1844    flagCol.put(1, flags1);
1845    ++iter;
1846  }
1847
1848  return out;
1849}
1850
1851std::vector< float > STMath::statistic( const CountedPtr< Scantable > & in,
1852                                        const std::vector< bool > & mask,
1853                                        const std::string& which )
1854{
1855
1856  Vector<Bool> m(mask);
1857  const Table& tab = in->table();
1858  ROArrayColumn<Float> specCol(tab, "SPECTRA");
1859  ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1860  std::vector<float> out;
1861  for (uInt i=0; i < tab.nrow(); ++i ) {
1862    Vector<Float> spec; specCol.get(i, spec);
1863    Vector<uChar> flag; flagCol.get(i, flag);
1864    MaskedArray<Float> ma  = maskedArray(spec, flag);
1865    float outstat = 0.0;
1866    if ( spec.nelements() == m.nelements() ) {
1867      outstat = mathutil::statistics(which, ma(m));
1868    } else {
1869      outstat = mathutil::statistics(which, ma);
1870    }
1871    out.push_back(outstat);
1872  }
1873  return out;
1874}
1875
1876std::vector< float > STMath::statisticRow( const CountedPtr< Scantable > & in,
1877                                        const std::vector< bool > & mask,
1878                                        const std::string& which,
1879                                        int row )
1880{
1881
1882  Vector<Bool> m(mask);
1883  const Table& tab = in->table();
1884  ROArrayColumn<Float> specCol(tab, "SPECTRA");
1885  ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1886  std::vector<float> out;
1887
1888  Vector<Float> spec; specCol.get(row, spec);
1889  Vector<uChar> flag; flagCol.get(row, flag);
1890  MaskedArray<Float> ma  = maskedArray(spec, flag);
1891  float outstat = 0.0;
1892  if ( spec.nelements() == m.nelements() ) {
1893    outstat = mathutil::statistics(which, ma(m));
1894  } else {
1895    outstat = mathutil::statistics(which, ma);
1896  }
1897  out.push_back(outstat);
1898
1899  return out;
1900}
1901
1902std::vector< int > STMath::minMaxChan( const CountedPtr< Scantable > & in,
1903                                        const std::vector< bool > & mask,
1904                                        const std::string& which )
1905{
1906
1907  Vector<Bool> m(mask);
1908  const Table& tab = in->table();
1909  ROArrayColumn<Float> specCol(tab, "SPECTRA");
1910  ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1911  std::vector<int> out;
1912  for (uInt i=0; i < tab.nrow(); ++i ) {
1913    Vector<Float> spec; specCol.get(i, spec);
1914    Vector<uChar> flag; flagCol.get(i, flag);
1915    MaskedArray<Float> ma  = maskedArray(spec, flag);
1916    if (ma.ndim() != 1) {
1917      throw (ArrayError(
1918          "std::vector<int> STMath::minMaxChan("
1919          "ContedPtr<Scantable> &in, std::vector<bool> &mask, "
1920          " std::string &which)"
1921          " - MaskedArray is not 1D"));
1922    }
1923    IPosition outpos(1,0);
1924    if ( spec.nelements() == m.nelements() ) {
1925      outpos = mathutil::minMaxPos(which, ma(m));
1926    } else {
1927      outpos = mathutil::minMaxPos(which, ma);
1928    }
1929    out.push_back(outpos[0]);
1930  }
1931  return out;
1932}
1933
1934CountedPtr< Scantable > STMath::bin( const CountedPtr< Scantable > & in,
1935                                     int width )
1936{
1937  if ( !in->getSelection().empty() ) throw(AipsError("Can't bin subset of the data."));
1938  CountedPtr< Scantable > out = getScantable(in, false);
1939  Table& tout = out->table();
1940  out->frequencies().rescale(width, "BIN");
1941  ArrayColumn<Float> specCol(tout, "SPECTRA");
1942  ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
1943  for (uInt i=0; i < tout.nrow(); ++i ) {
1944    MaskedArray<Float> main  = maskedArray(specCol(i), flagCol(i));
1945    MaskedArray<Float> maout;
1946    LatticeUtilities::bin(maout, main, 0, Int(width));
1947    /// @todo implement channel based tsys binning
1948    specCol.put(i, maout.getArray());
1949    flagCol.put(i, flagsFromMA(maout));
1950    // take only the first binned spectrum's length for the deprecated
1951    // global header item nChan
1952    if (i==0) tout.rwKeywordSet().define(String("nChan"),
1953                                       Int(maout.getArray().nelements()));
1954  }
1955  return out;
1956}
1957
1958CountedPtr< Scantable > STMath::resample( const CountedPtr< Scantable >& in,
1959                                          const std::string& method,
1960                                          float width )
1961//
1962// Should add the possibility of width being specified in km/s. This means
1963// that for each freqID (SpectralCoordinate) we will need to convert to an
1964// average channel width (say at the reference pixel).  Then we would need
1965// to be careful to make sure each spectrum (of different freqID)
1966// is the same length.
1967//
1968{
1969  //InterpolateArray1D<Double,Float>::InterpolationMethod interp;
1970  Int interpMethod(stringToIMethod(method));
1971
1972  CountedPtr< Scantable > out = getScantable(in, false);
1973  Table& tout = out->table();
1974
1975// Resample SpectralCoordinates (one per freqID)
1976  out->frequencies().rescale(width, "RESAMPLE");
1977  TableIterator iter(tout, "IFNO");
1978  TableRow row(tout);
1979  while ( !iter.pastEnd() ) {
1980    Table tab = iter.table();
1981    ArrayColumn<Float> specCol(tab, "SPECTRA");
1982    //ArrayColumn<Float> tsysCol(tout, "TSYS");
1983    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1984    Vector<Float> spec;
1985    Vector<uChar> flag;
1986    specCol.get(0,spec); // the number of channels should be constant per IF
1987    uInt nChanIn = spec.nelements();
1988    Vector<Float> xIn(nChanIn); indgen(xIn);
1989    Int fac =  Int(nChanIn/width);
1990    Vector<Float> xOut(fac+10); // 10 to be safe - resize later
1991    uInt k = 0;
1992    Float x = 0.0;
1993    while (x < Float(nChanIn) ) {
1994      xOut(k) = x;
1995      k++;
1996      x += width;
1997    }
1998    uInt nChanOut = k;
1999    xOut.resize(nChanOut, True);
2000    // process all rows for this IFNO
2001    Vector<Float> specOut;
2002    Vector<Bool> maskOut;
2003    Vector<uChar> flagOut;
2004    for (uInt i=0; i < tab.nrow(); ++i) {
2005      specCol.get(i, spec);
2006      flagCol.get(i, flag);
2007      Vector<Bool> mask(flag.nelements());
2008      convertArray(mask, flag);
2009
2010      IPosition shapeIn(spec.shape());
2011      //sh.nchan = nChanOut;
2012      InterpolateArray1D<Float,Float>::interpolate(specOut, maskOut, xOut,
2013                                                   xIn, spec, mask,
2014                                                   interpMethod, True, True);
2015      /// @todo do the same for channel based Tsys
2016      flagOut.resize(maskOut.nelements());
2017      convertArray(flagOut, maskOut);
2018      specCol.put(i, specOut);
2019      flagCol.put(i, flagOut);
2020    }
2021    ++iter;
2022  }
2023
2024  return out;
2025}
2026
2027STMath::imethod STMath::stringToIMethod(const std::string& in)
2028{
2029  static STMath::imap lookup;
2030
2031  // initialize the lookup table if necessary
2032  if ( lookup.empty() ) {
2033    lookup["nearest"]   = InterpolateArray1D<Double,Float>::nearestNeighbour;
2034    lookup["linear"] = InterpolateArray1D<Double,Float>::linear;
2035    lookup["cubic"]  = InterpolateArray1D<Double,Float>::cubic;
2036    lookup["spline"]  = InterpolateArray1D<Double,Float>::spline;
2037  }
2038
2039  STMath::imap::const_iterator iter = lookup.find(in);
2040
2041  if ( lookup.end() == iter ) {
2042    std::string message = in;
2043    message += " is not a valid interpolation mode";
2044    throw(AipsError(message));
2045  }
2046  return iter->second;
2047}
2048
2049WeightType STMath::stringToWeight(const std::string& in)
2050{
2051  static std::map<std::string, WeightType> lookup;
2052
2053  // initialize the lookup table if necessary
2054  if ( lookup.empty() ) {
2055    lookup["NONE"]   = asap::W_NONE;
2056    lookup["TINT"] = asap::W_TINT;
2057    lookup["TINTSYS"]  = asap::W_TINTSYS;
2058    lookup["TSYS"]  = asap::W_TSYS;
2059    lookup["VAR"]  = asap::W_VAR;
2060  }
2061
2062  std::map<std::string, WeightType>::const_iterator iter = lookup.find(in);
2063
2064  if ( lookup.end() == iter ) {
2065    std::string message = in;
2066    message += " is not a valid weighting mode";
2067    throw(AipsError(message));
2068  }
2069  return iter->second;
2070}
2071
2072CountedPtr< Scantable > STMath::gainElevation( const CountedPtr< Scantable >& in,
2073                                               const vector< float > & coeff,
2074                                               const std::string & filename,
2075                                               const std::string& method)
2076{
2077  // Get elevation data from Scantable and convert to degrees
2078  CountedPtr< Scantable > out = getScantable(in, false);
2079  Table& tab = out->table();
2080  ROScalarColumn<Float> elev(tab, "ELEVATION");
2081  Vector<Float> x = elev.getColumn();
2082  x *= Float(180 / C::pi);                        // Degrees
2083
2084  Vector<Float> coeffs(coeff);
2085  const uInt nc = coeffs.nelements();
2086  if ( filename.length() > 0 && nc > 0 ) {
2087    throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
2088  }
2089
2090  // Correct
2091  if ( nc > 0 || filename.length() == 0 ) {
2092    // Find instrument
2093    Bool throwit = True;
2094    Instrument inst =
2095      STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
2096                                throwit);
2097
2098    // Set polynomial
2099    Polynomial<Float>* ppoly = 0;
2100    Vector<Float> coeff;
2101    String msg;
2102    if ( nc > 0 ) {
2103      ppoly = new Polynomial<Float>(nc-1);
2104      coeff = coeffs;
2105      msg = String("user");
2106    } else {
2107      STAttr sdAttr;
2108      coeff = sdAttr.gainElevationPoly(inst);
2109      ppoly = new Polynomial<Float>(coeff.nelements()-1);
2110      msg = String("built in");
2111    }
2112
2113    if ( coeff.nelements() > 0 ) {
2114      ppoly->setCoefficients(coeff);
2115    } else {
2116      delete ppoly;
2117      throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
2118    }
2119    ostringstream oss;
2120    oss << "Making polynomial correction with " << msg << " coefficients:" << endl;
2121    oss << "   " <<  coeff;
2122    pushLog(String(oss));
2123    const uInt nrow = tab.nrow();
2124    Vector<Float> factor(nrow);
2125    for ( uInt i=0; i < nrow; ++i ) {
2126      factor[i] = 1.0 / (*ppoly)(x[i]);
2127    }
2128    delete ppoly;
2129    scaleByVector(tab, factor, true);
2130
2131  } else {
2132    // Read and correct
2133    pushLog("Making correction from ascii Table");
2134    scaleFromAsciiTable(tab, filename, method, x, true);
2135  }
2136  return out;
2137}
2138
2139void STMath::scaleFromAsciiTable(Table& in, const std::string& filename,
2140                                 const std::string& method,
2141                                 const Vector<Float>& xout, bool dotsys)
2142{
2143
2144// Read gain-elevation ascii file data into a Table.
2145
2146  String formatString;
2147  Table tbl = readAsciiTable(formatString, Table::Memory, filename, "", "", False);
2148  scaleFromTable(in, tbl, method, xout, dotsys);
2149}
2150
2151void STMath::scaleFromTable(Table& in,
2152                            const Table& table,
2153                            const std::string& method,
2154                            const Vector<Float>& xout, bool dotsys)
2155{
2156
2157  ROScalarColumn<Float> geElCol(table, "ELEVATION");
2158  ROScalarColumn<Float> geFacCol(table, "FACTOR");
2159  Vector<Float> xin = geElCol.getColumn();
2160  Vector<Float> yin = geFacCol.getColumn();
2161  Vector<Bool> maskin(xin.nelements(),True);
2162
2163  // Interpolate (and extrapolate) with desired method
2164
2165  InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
2166
2167   Vector<Float> yout;
2168   Vector<Bool> maskout;
2169   InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout,
2170                                                xin, yin, maskin, interp,
2171                                                True, True);
2172
2173   scaleByVector(in, Float(1.0)/yout, dotsys);
2174}
2175
2176void STMath::scaleByVector( Table& in,
2177                            const Vector< Float >& factor,
2178                            bool dotsys )
2179{
2180  uInt nrow = in.nrow();
2181  if ( factor.nelements() != nrow ) {
2182    throw(AipsError("factors.nelements() != table.nelements()"));
2183  }
2184  ArrayColumn<Float> specCol(in, "SPECTRA");
2185  ArrayColumn<uChar> flagCol(in, "FLAGTRA");
2186  ArrayColumn<Float> tsysCol(in, "TSYS");
2187  for (uInt i=0; i < nrow; ++i) {
2188    MaskedArray<Float> ma  = maskedArray(specCol(i), flagCol(i));
2189    ma *= factor[i];
2190    specCol.put(i, ma.getArray());
2191    flagCol.put(i, flagsFromMA(ma));
2192    if ( dotsys ) {
2193      Vector<Float> tsys = tsysCol(i);
2194      tsys *= factor[i];
2195      tsysCol.put(i,tsys);
2196    }
2197  }
2198}
2199
2200CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in,
2201                                             float d, float etaap,
2202                                             float jyperk )
2203{
2204  CountedPtr< Scantable > out = getScantable(in, false);
2205  Table& tab = in->table();
2206  Unit fluxUnit(tab.keywordSet().asString("FluxUnit"));
2207  Unit K(String("K"));
2208  Unit JY(String("Jy"));
2209
2210  bool tokelvin = true;
2211  Double cfac = 1.0;
2212
2213  if ( fluxUnit == JY ) {
2214    pushLog("Converting to K");
2215    Quantum<Double> t(1.0,fluxUnit);
2216    Quantum<Double> t2 = t.get(JY);
2217    cfac = (t2 / t).getValue();               // value to Jy
2218
2219    tokelvin = true;
2220    out->setFluxUnit("K");
2221  } else if ( fluxUnit == K ) {
2222    pushLog("Converting to Jy");
2223    Quantum<Double> t(1.0,fluxUnit);
2224    Quantum<Double> t2 = t.get(K);
2225    cfac = (t2 / t).getValue();              // value to K
2226
2227    tokelvin = false;
2228    out->setFluxUnit("Jy");
2229  } else {
2230    throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
2231  }
2232  // Make sure input values are converted to either Jy or K first...
2233  Float factor = cfac;
2234
2235  // Select method
2236  if (jyperk > 0.0) {
2237    factor *= jyperk;
2238    if ( tokelvin ) factor = 1.0 / jyperk;
2239    ostringstream oss;
2240    oss << "Jy/K = " << jyperk;
2241    pushLog(String(oss));
2242    Vector<Float> factors(tab.nrow(), factor);
2243    scaleByVector(tab,factors, false);
2244  } else if ( etaap > 0.0) {
2245    if (d < 0) {
2246      Instrument inst =
2247        STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
2248                                  True);
2249      STAttr sda;
2250      d = sda.diameter(inst);
2251    }
2252    jyperk = STAttr::findJyPerK(etaap, d);
2253    ostringstream oss;
2254    oss << "Jy/K = " << jyperk;
2255    pushLog(String(oss));
2256    factor *= jyperk;
2257    if ( tokelvin ) {
2258      factor = 1.0 / factor;
2259    }
2260    Vector<Float> factors(tab.nrow(), factor);
2261    scaleByVector(tab, factors, False);
2262  } else {
2263
2264    // OK now we must deal with automatic look up of values.
2265    // We must also deal with the fact that the factors need
2266    // to be computed per IF and may be different and may
2267    // change per integration.
2268
2269    pushLog("Looking up conversion factors");
2270    convertBrightnessUnits(out, tokelvin, cfac);
2271  }
2272
2273  return out;
2274}
2275
2276void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in,
2277                                     bool tokelvin, float cfac )
2278{
2279  Table& table = in->table();
2280  Instrument inst =
2281    STAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True);
2282  TableIterator iter(table, "FREQ_ID");
2283  STFrequencies stfreqs = in->frequencies();
2284  STAttr sdAtt;
2285  while (!iter.pastEnd()) {
2286    Table tab = iter.table();
2287    ArrayColumn<Float> specCol(tab, "SPECTRA");
2288    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2289    ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID");
2290    MEpoch::ROScalarColumn timeCol(tab, "TIME");
2291
2292    uInt freqid; freqidCol.get(0, freqid);
2293    Vector<Float> tmpspec; specCol.get(0, tmpspec);
2294    // STAttr.JyPerK has a Vector interface... change sometime.
2295    Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements()));
2296    for ( uInt i=0; i<tab.nrow(); ++i) {
2297      Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0];
2298      Float factor = cfac * jyperk;
2299      if ( tokelvin ) factor = Float(1.0) / factor;
2300      MaskedArray<Float> ma  = maskedArray(specCol(i), flagCol(i));
2301      ma *= factor;
2302      specCol.put(i, ma.getArray());
2303      flagCol.put(i, flagsFromMA(ma));
2304    }
2305  ++iter;
2306  }
2307}
2308
2309CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
2310                                         const std::vector<float>& tau )
2311{
2312  CountedPtr< Scantable > out = getScantable(in, false);
2313
2314  Table outtab = out->table();
2315
2316  const Int ntau = uInt(tau.size());
2317  std::vector<float>::const_iterator tauit = tau.begin();
2318  AlwaysAssert((ntau == 1 || ntau == in->nif() || ntau == in->nif() * in->npol()),
2319               AipsError);
2320  TableIterator iiter(outtab, "IFNO");
2321  while ( !iiter.pastEnd() ) {
2322    Table itab = iiter.table();
2323    TableIterator piter(itab, "POLNO");
2324    while ( !piter.pastEnd() ) {
2325      Table tab = piter.table();
2326      ROScalarColumn<Float> elev(tab, "ELEVATION");
2327      ArrayColumn<Float> specCol(tab, "SPECTRA");
2328      ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2329      ArrayColumn<Float> tsysCol(tab, "TSYS");
2330      for ( uInt i=0; i<tab.nrow(); ++i) {
2331        Float zdist = Float(C::pi_2) - elev(i);
2332        Float factor = exp(*tauit/cos(zdist));
2333        MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
2334        ma *= factor;
2335        specCol.put(i, ma.getArray());
2336        flagCol.put(i, flagsFromMA(ma));
2337        Vector<Float> tsys;
2338        tsysCol.get(i, tsys);
2339        tsys *= factor;
2340        tsysCol.put(i, tsys);
2341      }
2342      if (ntau == in->nif()*in->npol() ) {
2343        tauit++;
2344      }
2345      piter++;
2346    }
2347    if (ntau >= in->nif() ) {
2348      tauit++;
2349    }
2350    iiter++;
2351  }
2352  return out;
2353}
2354
2355CountedPtr< Scantable > STMath::smoothOther( const CountedPtr< Scantable >& in,
2356                                             const std::string& kernel,
2357                                             float width, int order)
2358{
2359  CountedPtr< Scantable > out = getScantable(in, false);
2360  Table table = out->table();
2361
2362  TableIterator iter(table, "IFNO");
2363  while (!iter.pastEnd()) {
2364    Table tab = iter.table();
2365    ArrayColumn<Float> specCol(tab, "SPECTRA");
2366    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2367    Vector<Float> spec;
2368    Vector<uChar> flag;
2369    for (uInt i = 0; i < tab.nrow(); ++i) {
2370      specCol.get(i, spec);
2371      flagCol.get(i, flag);
2372      Vector<Bool> mask(flag.nelements());
2373      convertArray(mask, flag);
2374      Vector<Float> specout;
2375      Vector<Bool> maskout;
2376      if (kernel == "hanning") {
2377        mathutil::hanning(specout, maskout, spec, !mask);
2378        convertArray(flag, !maskout);
2379      } else if (kernel == "rmedian") {
2380        mathutil::runningMedian(specout, maskout, spec , mask, width);
2381        convertArray(flag, maskout);
2382      } else if (kernel == "poly") {
2383        mathutil::polyfit(specout, maskout, spec, !mask, width, order);
2384        convertArray(flag, !maskout);
2385      }
2386
2387      for (uInt j = 0; j < flag.nelements(); ++j) {
2388        uChar userFlag = 1 << 7;
2389        if (maskout[j]==True) userFlag = 0 << 7;
2390        flag(j) = userFlag;
2391      }
2392
2393      flagCol.put(i, flag);
2394      specCol.put(i, specout);
2395    }
2396  ++iter;
2397  }
2398  return out;
2399}
2400
2401CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
2402                                        const std::string& kernel, float width,
2403                                        int order)
2404{
2405  if (kernel == "rmedian"  || kernel == "hanning" || kernel == "poly") {
2406    return smoothOther(in, kernel, width, order);
2407  }
2408  CountedPtr< Scantable > out = getScantable(in, false);
2409  Table& table = out->table();
2410  VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
2411  // same IFNO should have same no of channels
2412  // this saves overhead
2413  TableIterator iter(table, "IFNO");
2414  while (!iter.pastEnd()) {
2415    Table tab = iter.table();
2416    ArrayColumn<Float> specCol(tab, "SPECTRA");
2417    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2418    Vector<Float> tmpspec; specCol.get(0, tmpspec);
2419    uInt nchan = tmpspec.nelements();
2420    Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
2421    Convolver<Float> conv(kvec, IPosition(1,nchan));
2422    Vector<Float> spec;
2423    Vector<uChar> flag;
2424    for ( uInt i=0; i<tab.nrow(); ++i) {
2425      specCol.get(i, spec);
2426      flagCol.get(i, flag);
2427      Vector<Bool> mask(flag.nelements());
2428      convertArray(mask, flag);
2429      Vector<Float> specout;
2430      mathutil::replaceMaskByZero(specout, mask);
2431      conv.linearConv(specout, spec);
2432      specCol.put(i, specout);
2433    }
2434    ++iter;
2435  }
2436  return out;
2437}
2438
2439CountedPtr< Scantable >
2440  STMath::merge( const std::vector< CountedPtr < Scantable > >& in )
2441{
2442  if ( in.size() < 2 ) {
2443    throw(AipsError("Need at least two scantables to perform a merge."));
2444  }
2445  std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
2446  bool insitu = insitu_;
2447  setInsitu(false);
2448  CountedPtr< Scantable > out = getScantable(*it, false);
2449  setInsitu(insitu);
2450  Table& tout = out->table();
2451  ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
2452  ScalarColumn<uInt> scannocol(tout,"SCANNO"), focusidcol(tout,"FOCUS_ID");
2453  // Renumber SCANNO to be 0-based
2454  Vector<uInt> scannos = scannocol.getColumn();
2455  uInt offset = min(scannos);
2456  scannos -= offset;
2457  scannocol.putColumn(scannos);
2458  uInt newscanno = max(scannos)+1;
2459  ++it;
2460  while ( it != in.end() ){
2461    if ( ! (*it)->conformant(*out) ) {
2462      // non conformant.
2463      //pushLog(String("Warning: Can't merge scantables as header info differs."));
2464      LogIO os( LogOrigin( "STMath", "merge()", WHERE ) ) ;
2465      os << LogIO::SEVERE << "Can't merge scantables as header informations (any one of AntennaName, Equinox, and FluxUnit) differ." << LogIO::EXCEPTION ;
2466    }
2467    out->appendToHistoryTable((*it)->history());
2468    const Table& tab = (*it)->table();
2469
2470    Block<String> cols(3);
2471    cols[0] = String("FREQ_ID");
2472    cols[1] = String("MOLECULE_ID");
2473    cols[2] = String("FOCUS_ID");
2474
2475    TableIterator scanit(tab, "SCANNO");
2476    while (!scanit.pastEnd()) {
2477      ScalarColumn<uInt> thescannocol(scanit.table(),"SCANNO");
2478      Vector<uInt> thescannos(thescannocol.nrow(),newscanno);
2479      thescannocol.putColumn(thescannos);
2480      TableIterator subit(scanit.table(), cols);
2481      while ( !subit.pastEnd() ) {
2482        uInt nrow = tout.nrow();
2483        Table thetab = subit.table();
2484        ROTableRow row(thetab);
2485        Vector<uInt> thecolvals(thetab.nrow());
2486        ScalarColumn<uInt> thefreqidcol(thetab,"FREQ_ID");
2487        ScalarColumn<uInt> themolidcol(thetab, "MOLECULE_ID");
2488        ScalarColumn<uInt> thefocusidcol(thetab,"FOCUS_ID");
2489        // The selected subset of table should have
2490        // the equal FREQ_ID, MOLECULE_ID, and FOCUS_ID values.
2491        const TableRecord& rec = row.get(0);
2492        // Set the proper FREQ_ID
2493        Double rv,rp,inc;
2494        (*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
2495        uInt id;
2496        id = out->frequencies().addEntry(rp, rv, inc);
2497        thecolvals = id;
2498        thefreqidcol.putColumn(thecolvals);
2499        // Set the proper MOLECULE_ID
2500        Vector<String> name,fname;Vector<Double> rf;
2501        (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
2502        id = out->molecules().addEntry(rf, name, fname);
2503        thecolvals = id;
2504        themolidcol.putColumn(thecolvals);
2505        // Set the proper FOCUS_ID
2506        Float fpa,frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
2507        (*it)->focus().getEntry(fpa, fax, ftan, frot, fhand, fmount,fuser,
2508                                fxy, fxyp, rec.asuInt("FOCUS_ID"));
2509        id = out->focus().addEntry(fpa, fax, ftan, frot, fhand, fmount,fuser,
2510                                   fxy, fxyp);
2511        thecolvals = id;
2512        thefocusidcol.putColumn(thecolvals);
2513
2514        tout.addRow(thetab.nrow());
2515        TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
2516
2517        ++subit;
2518      }
2519      ++newscanno;
2520      ++scanit;
2521    }
2522    ++it;
2523  }
2524  return out;
2525}
2526
2527CountedPtr< Scantable >
2528  STMath::invertPhase( const CountedPtr < Scantable >& in )
2529{
2530  return applyToPol(in, &STPol::invertPhase, Float(0.0));
2531}
2532
2533CountedPtr< Scantable >
2534  STMath::rotateXYPhase( const CountedPtr < Scantable >& in, float phase )
2535{
2536   return applyToPol(in, &STPol::rotatePhase, Float(phase));
2537}
2538
2539CountedPtr< Scantable >
2540  STMath::rotateLinPolPhase( const CountedPtr < Scantable >& in, float phase )
2541{
2542  return applyToPol(in, &STPol::rotateLinPolPhase, Float(phase));
2543}
2544
2545CountedPtr< Scantable > STMath::applyToPol( const CountedPtr<Scantable>& in,
2546                                             STPol::polOperation fptr,
2547                                             Float phase )
2548{
2549  CountedPtr< Scantable > out = getScantable(in, false);
2550  Table& tout = out->table();
2551  Block<String> cols(4);
2552  cols[0] = String("SCANNO");
2553  cols[1] = String("BEAMNO");
2554  cols[2] = String("IFNO");
2555  cols[3] = String("CYCLENO");
2556  TableIterator iter(tout, cols);
2557  CountedPtr<STPol> stpol = STPol::getPolClass(out->factories_,
2558                                               out->getPolType() );
2559  while (!iter.pastEnd()) {
2560    Table t = iter.table();
2561    ArrayColumn<Float> speccol(t, "SPECTRA");
2562    ScalarColumn<uInt> focidcol(t, "FOCUS_ID");
2563    Matrix<Float> pols(speccol.getColumn());
2564    try {
2565      stpol->setSpectra(pols);
2566      Float fang,fhand;
2567      fang = in->focusTable_.getTotalAngle(focidcol(0));
2568      fhand = in->focusTable_.getFeedHand(focidcol(0));
2569      stpol->setPhaseCorrections(fang, fhand);
2570      // use a member function pointer in STPol.  This only works on
2571      // the STPol pointer itself, not the Counted Pointer so
2572      // derefernce it.
2573      (&(*(stpol))->*fptr)(phase);
2574      speccol.putColumn(stpol->getSpectra());
2575    } catch (AipsError& e) {
2576      //delete stpol;stpol=0;
2577      throw(e);
2578    }
2579    ++iter;
2580  }
2581  //delete stpol;stpol=0;
2582  return out;
2583}
2584
2585CountedPtr< Scantable >
2586  STMath::swapPolarisations( const CountedPtr< Scantable > & in )
2587{
2588  CountedPtr< Scantable > out = getScantable(in, false);
2589  Table& tout = out->table();
2590  Table t0 = tout(tout.col("POLNO") == 0);
2591  Table t1 = tout(tout.col("POLNO") == 1);
2592  if ( t0.nrow() != t1.nrow() )
2593    throw(AipsError("Inconsistent number of polarisations"));
2594  ArrayColumn<Float> speccol0(t0, "SPECTRA");
2595  ArrayColumn<uChar> flagcol0(t0, "FLAGTRA");
2596  ArrayColumn<Float> speccol1(t1, "SPECTRA");
2597  ArrayColumn<uChar> flagcol1(t1, "FLAGTRA");
2598  Matrix<Float> s0 = speccol0.getColumn();
2599  Matrix<uChar> f0 = flagcol0.getColumn();
2600  speccol0.putColumn(speccol1.getColumn());
2601  flagcol0.putColumn(flagcol1.getColumn());
2602  speccol1.putColumn(s0);
2603  flagcol1.putColumn(f0);
2604  return out;
2605}
2606
2607CountedPtr< Scantable >
2608  STMath::averagePolarisations( const CountedPtr< Scantable > & in,
2609                                const std::vector<bool>& mask,
2610                                const std::string& weight )
2611{
2612  if (in->npol() < 2 )
2613    throw(AipsError("averagePolarisations can only be applied to two or more"
2614                    "polarisations"));
2615  bool insitu = insitu_;
2616  setInsitu(false);
2617  CountedPtr< Scantable > pols = getScantable(in, true);
2618  setInsitu(insitu);
2619  Table& tout = pols->table();
2620  std::string taql = "SELECT FROM $1 WHERE POLNO IN [0,1]";
2621  Table tab = tableCommand(taql, in->table());
2622  if (tab.nrow() == 0 )
2623    throw(AipsError("Could not find  any rows with POLNO==0 and POLNO==1"));
2624  TableCopy::copyRows(tout, tab);
2625  TableVector<uInt> vec(tout, "POLNO");
2626  vec = 0;
2627  pols->table_.rwKeywordSet().define("nPol", Int(1));
2628  pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
2629  //pols->table_.rwKeywordSet().define("POLTYPE", in->getPolType());
2630  std::vector<CountedPtr<Scantable> > vpols;
2631  vpols.push_back(pols);
2632  CountedPtr< Scantable > out = average(vpols, mask, weight, "SCAN");
2633  return out;
2634}
2635
2636CountedPtr< Scantable >
2637  STMath::averageBeams( const CountedPtr< Scantable > & in,
2638                        const std::vector<bool>& mask,
2639                        const std::string& weight )
2640{
2641  bool insitu = insitu_;
2642  setInsitu(false);
2643  CountedPtr< Scantable > beams = getScantable(in, false);
2644  setInsitu(insitu);
2645  Table& tout = beams->table();
2646  // give all rows the same BEAMNO
2647  TableVector<uInt> vec(tout, "BEAMNO");
2648  vec = 0;
2649  beams->table_.rwKeywordSet().define("nBeam", Int(1));
2650  std::vector<CountedPtr<Scantable> > vbeams;
2651  vbeams.push_back(beams);
2652  CountedPtr< Scantable > out = average(vbeams, mask, weight, "SCAN");
2653  return out;
2654}
2655
2656
2657CountedPtr< Scantable >
2658  asap::STMath::frequencyAlign( const CountedPtr< Scantable > & in,
2659                                const std::string & refTime,
2660                                const std::string & method)
2661{
2662  // clone as this is not working insitu
2663  bool insitu = insitu_;
2664  setInsitu(false);
2665  CountedPtr< Scantable > out = getScantable(in, false);
2666  setInsitu(insitu);
2667  Table& tout = out->table();
2668  // Get reference Epoch to time of first row or given String
2669  Unit DAY(String("d"));
2670  MEpoch::Ref epochRef(in->getTimeReference());
2671  MEpoch refEpoch;
2672  if (refTime.length()>0) {
2673    Quantum<Double> qt;
2674    if (MVTime::read(qt,refTime)) {
2675      MVEpoch mv(qt);
2676      refEpoch = MEpoch(mv, epochRef);
2677   } else {
2678      throw(AipsError("Invalid format for Epoch string"));
2679   }
2680  } else {
2681    refEpoch = in->timeCol_(0);
2682  }
2683  MPosition refPos = in->getAntennaPosition();
2684
2685  InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
2686  /*
2687  // Comment from MV.
2688  // the following code has been commented out because different FREQ_IDs have to be aligned together even
2689  // if the frame doesn't change. So far, lack of this check didn't cause any problems.
2690  // test if user frame is different to base frame
2691  if ( in->frequencies().getFrameString(true)
2692       == in->frequencies().getFrameString(false) ) {
2693    throw(AipsError("Can't convert as no output frame has been set"
2694                    " (use set_freqframe) or it is aligned already."));
2695  }
2696  */
2697  MFrequency::Types system = in->frequencies().getFrame();
2698  MVTime mvt(refEpoch.getValue());
2699  String epochout = mvt.string(MVTime::YMD) + String(" (") + refEpoch.getRefString() + String(")");
2700  ostringstream oss;
2701  oss << "Aligned at reference Epoch " << epochout
2702      << " in frame " << MFrequency::showType(system);
2703  pushLog(String(oss));
2704  // set up the iterator
2705  Block<String> cols(4);
2706  // select by constant direction
2707  cols[0] = String("SRCNAME");
2708  cols[1] = String("BEAMNO");
2709  // select by IF ( no of channels varies over this )
2710  cols[2] = String("IFNO");
2711  // select by restfrequency
2712  cols[3] = String("MOLECULE_ID");
2713  TableIterator iter(tout, cols);
2714  while ( !iter.pastEnd() ) {
2715    Table t = iter.table();
2716    MDirection::ROScalarColumn dirCol(t, "DIRECTION");
2717    TableIterator fiter(t, "FREQ_ID");
2718    // determine nchan from the first row. This should work as
2719    // we are iterating over BEAMNO and IFNO    // we should have constant direction
2720
2721    ROArrayColumn<Float> sCol(t, "SPECTRA");
2722    const MDirection direction = dirCol(0);
2723    const uInt nchan = sCol(0).nelements();
2724
2725    // skip operations if there is nothing to align
2726    if (fiter.pastEnd()) {
2727        continue;
2728    }
2729
2730    Table ftab = fiter.table();
2731    // align all frequency ids with respect to the first encountered id
2732    ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
2733    // get the SpectralCoordinate for the freqid, which we are iterating over
2734    SpectralCoordinate sC = in->frequencies().getSpectralCoordinate(freqidCol(0));
2735    FrequencyAligner<Float> fa( sC, nchan, refEpoch,
2736                                direction, refPos, system );
2737    // realign the SpectralCoordinate and put into the output Scantable
2738    Vector<String> units(1);
2739    units = String("Hz");
2740    Bool linear=True;
2741    SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
2742    sc2.setWorldAxisUnits(units);
2743    const uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
2744                                                sc2.referenceValue()[0],
2745                                                sc2.increment()[0]);
2746    while ( !fiter.pastEnd() ) {
2747      ftab = fiter.table();
2748      // spectral coordinate for the current FREQ_ID
2749      ScalarColumn<uInt> freqidCol2(ftab, "FREQ_ID");
2750      sC = in->frequencies().getSpectralCoordinate(freqidCol2(0));
2751      // create the "global" abcissa for alignment with same FREQ_ID
2752      Vector<Double> abc(nchan);
2753      for (uInt i=0; i<nchan; i++) {
2754           Double w;
2755           sC.toWorld(w,Double(i));
2756           abc[i] = w;
2757      }
2758      TableVector<uInt> tvec(ftab, "FREQ_ID");
2759      // assign new frequency id to all rows
2760      tvec = id;
2761      // cache abcissa for same time stamps, so iterate over those
2762      TableIterator timeiter(ftab, "TIME");
2763      while ( !timeiter.pastEnd() ) {
2764        Table tab = timeiter.table();
2765        ArrayColumn<Float> specCol(tab, "SPECTRA");
2766        ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2767        MEpoch::ROScalarColumn timeCol(tab, "TIME");
2768        // use align abcissa cache after the first row
2769        // these rows should be just be POLNO
2770        bool first = true;
2771        for (int i=0; i<int(tab.nrow()); ++i) {
2772          // input values
2773          Vector<uChar> flag = flagCol(i);
2774          Vector<Bool> mask(flag.shape());
2775          Vector<Float> specOut, spec;
2776          spec  = specCol(i);
2777          Vector<Bool> maskOut;Vector<uChar> flagOut;
2778          convertArray(mask, flag);
2779          // alignment
2780          Bool ok = fa.align(specOut, maskOut, abc, spec,
2781                             mask, timeCol(i), !first,
2782                             interp, False);
2783          (void) ok; // unused stop compiler nagging     
2784          // back into scantable
2785          flagOut.resize(maskOut.nelements());
2786          convertArray(flagOut, maskOut);
2787          flagCol.put(i, flagOut);
2788          specCol.put(i, specOut);
2789          // start abcissa caching
2790          first = false;
2791        }
2792        // next timestamp
2793        ++timeiter;
2794      }
2795      // next FREQ_ID
2796      ++fiter;
2797    }
2798    // next aligner
2799    ++iter;
2800  }
2801  // set this afterwards to ensure we are doing insitu correctly.
2802  out->frequencies().setFrame(system, true);
2803  return out;
2804}
2805
2806CountedPtr<Scantable>
2807  asap::STMath::convertPolarisation( const CountedPtr<Scantable>& in,
2808                                     const std::string & newtype )
2809{
2810  if (in->npol() != 2 && in->npol() != 4)
2811    throw(AipsError("Can only convert two or four polarisations."));
2812  if ( in->getPolType() == newtype )
2813    throw(AipsError("No need to convert."));
2814  if ( ! in->selector_.empty() )
2815    throw(AipsError("Can only convert whole scantable. Unset the selection."));
2816  bool insitu = insitu_;
2817  setInsitu(false);
2818  CountedPtr< Scantable > out = getScantable(in, true);
2819  setInsitu(insitu);
2820  Table& tout = out->table();
2821  tout.rwKeywordSet().define("POLTYPE", String(newtype));
2822
2823  Block<String> cols(4);
2824  cols[0] = "SCANNO";
2825  cols[1] = "CYCLENO";
2826  cols[2] = "BEAMNO";
2827  cols[3] = "IFNO";
2828  TableIterator it(in->originalTable_, cols);
2829  String basetype = in->getPolType();
2830  STPol* stpol = STPol::getPolClass(in->factories_, basetype);
2831  try {
2832    while ( !it.pastEnd() ) {
2833      Table tab = it.table();
2834      uInt row = tab.rowNumbers()[0];
2835      stpol->setSpectra(in->getPolMatrix(row));
2836      Float fang,fhand;
2837      fang = in->focusTable_.getTotalAngle(in->mfocusidCol_(row));
2838      fhand = in->focusTable_.getFeedHand(in->mfocusidCol_(row));
2839      stpol->setPhaseCorrections(fang, fhand);
2840      Int npolout = 0;
2841      for (uInt i=0; i<tab.nrow(); ++i) {
2842        Vector<Float> outvec = stpol->getSpectrum(i, newtype);
2843        if ( outvec.nelements() > 0 ) {
2844          tout.addRow();
2845          TableCopy::copyRows(tout, tab, tout.nrow()-1, 0, 1);
2846          ArrayColumn<Float> sCol(tout,"SPECTRA");
2847          ScalarColumn<uInt> pCol(tout,"POLNO");
2848          sCol.put(tout.nrow()-1 ,outvec);
2849          pCol.put(tout.nrow()-1 ,uInt(npolout));
2850          npolout++;
2851       }
2852      }
2853      tout.rwKeywordSet().define("nPol", npolout);
2854      ++it;
2855    }
2856  } catch (AipsError& e) {
2857    delete stpol;
2858    throw(e);
2859  }
2860  delete stpol;
2861  return out;
2862}
2863
2864CountedPtr< Scantable >
2865  asap::STMath::mxExtract( const CountedPtr< Scantable > & in,
2866                           const std::string & scantype )
2867{
2868  bool insitu = insitu_;
2869  setInsitu(false);
2870  CountedPtr< Scantable > out = getScantable(in, true);
2871  setInsitu(insitu);
2872  Table& tout = out->table();
2873  std::string taql = "SELECT FROM $1 WHERE BEAMNO != REFBEAMNO";
2874  if (scantype == "on") {
2875    taql = "SELECT FROM $1 WHERE BEAMNO == REFBEAMNO";
2876  }
2877  Table tab = tableCommand(taql, in->table());
2878  TableCopy::copyRows(tout, tab);
2879  if (scantype == "on") {
2880    // re-index SCANNO to 0
2881    TableVector<uInt> vec(tout, "SCANNO");
2882    vec = 0;
2883  }
2884  return out;
2885}
2886
2887std::vector<float>
2888  asap::STMath::fft( const casa::CountedPtr< Scantable > & in,
2889                     const std::vector<int>& whichrow,
2890                     bool getRealImag )
2891{
2892  std::vector<float> res;
2893  Table tab = in->table();
2894  std::vector<bool> mask;
2895
2896  if (whichrow.size() < 1) {          // for all rows (by default)
2897    int nrow = int(tab.nrow());
2898    for (int i = 0; i < nrow; ++i) {
2899      res = in->execFFT(i, mask, getRealImag);
2900    }
2901  } else {                           // for specified rows
2902    for (uInt i = 0; i < whichrow.size(); ++i) {
2903      res = in->execFFT(i, mask, getRealImag);
2904    }
2905  }
2906
2907  return res;
2908}
2909
2910
2911CountedPtr<Scantable>
2912  asap::STMath::lagFlag( const CountedPtr<Scantable>& in,
2913                         double start, double end,
2914                         const std::string& mode )
2915{
2916  CountedPtr<Scantable> out = getScantable(in, false);
2917  Table& tout = out->table();
2918  TableIterator iter(tout, "FREQ_ID");
2919  FFTServer<Float,Complex> ffts;
2920
2921  while ( !iter.pastEnd() ) {
2922    Table tab = iter.table();
2923    Double rp,rv,inc;
2924    ROTableRow row(tab);
2925    const TableRecord& rec = row.get(0);
2926    uInt freqid = rec.asuInt("FREQ_ID");
2927    out->frequencies().getEntry(rp, rv, inc, freqid);
2928    ArrayColumn<Float> specCol(tab, "SPECTRA");
2929    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2930
2931    for (int i=0; i<int(tab.nrow()); ++i) {
2932      Vector<Float> spec = specCol(i);
2933      Vector<uChar> flag = flagCol(i);
2934      std::vector<bool> mask;
2935      for (uInt j = 0; j < flag.nelements(); ++j) {
2936        mask.push_back(!(flag[j]>0));
2937      }
2938      mathutil::doZeroOrderInterpolation(spec, mask);
2939
2940      Vector<Complex> lags;
2941      ffts.fft0(lags, spec);
2942
2943      Int lag0(start+0.5);
2944      Int lag1(end+0.5);
2945      if (mode == "frequency") {
2946        lag0 = Int(spec.nelements()*abs(inc)/(start)+0.5);
2947        lag1 = Int(spec.nelements()*abs(inc)/(end)+0.5);
2948      }
2949      Int lstart =  max(0, lag0);
2950      Int lend   =  min(Int(lags.nelements()-1), lag1);
2951      if (lstart == lend) {
2952        lags[lstart] = Complex(0.0);
2953      } else {
2954        if (lstart > lend) {
2955          Int tmp = lend;
2956          lend = lstart;
2957          lstart = tmp;
2958        }
2959        for (int j=lstart; j <=lend ;++j) {
2960          lags[j] = Complex(0.0);
2961        }
2962      }
2963
2964      ffts.fft0(spec, lags);
2965
2966      specCol.put(i, spec);
2967    }
2968    ++iter;
2969  }
2970  return out;
2971}
2972
2973// Averaging spectra with different channel/resolution
2974CountedPtr<Scantable>
2975STMath::new_average( const std::vector<CountedPtr<Scantable> >& in,
2976                     const bool& compel,
2977                     const std::vector<bool>& mask,
2978                     const std::string& weight,
2979                     const std::string& avmode )
2980  throw ( casa::AipsError )
2981{
2982  LogIO os( LogOrigin( "STMath", "new_average()", WHERE ) ) ;
2983  if ( avmode == "SCAN" && in.size() != 1 )
2984    throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
2985                    "Use merge first."));
2986 
2987  // 2012/02/17 TN
2988  // Since STGrid is implemented, average doesn't consider direction
2989  // when accumulating
2990  // check if OTF observation
2991//   String obstype = in[0]->getHeader().obstype ;
2992//   Double tol = 0.0 ;
2993//   if ( obstype.find( "OTF" ) != String::npos ) {
2994//     tol = TOL_OTF ;
2995//   }
2996//   else {
2997//     tol = TOL_POINT ;
2998//   }
2999
3000  CountedPtr<Scantable> out ;     // processed result
3001  if ( compel ) {
3002    std::vector< CountedPtr<Scantable> > newin ; // input for average process
3003    uInt insize = in.size() ;    // number of input scantables
3004
3005    // TEST: do normal average in each table before IF grouping
3006    os << "Do preliminary averaging" << LogIO::POST ;
3007    vector< CountedPtr<Scantable> > tmpin( insize ) ;
3008    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3009      vector< CountedPtr<Scantable> > v( 1, in[itable] ) ;
3010      tmpin[itable] = average( v, mask, weight, avmode ) ;
3011    }
3012
3013    // warning
3014    os << "Average spectra with different spectral resolution" << LogIO::POST ;
3015
3016    // temporarily set coordinfo
3017    vector<string> oldinfo( insize ) ;
3018    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3019      vector<string> coordinfo = in[itable]->getCoordInfo() ;
3020      oldinfo[itable] = coordinfo[0] ;
3021      coordinfo[0] = "Hz" ;
3022      tmpin[itable]->setCoordInfo( coordinfo ) ;
3023    }
3024
3025    // columns
3026    ScalarColumn<uInt> freqIDCol ;
3027    ScalarColumn<uInt> ifnoCol ;
3028    ScalarColumn<uInt> scannoCol ;
3029
3030
3031    // check IF frequency coverage
3032    // freqid: list of FREQ_ID, which is used, in each table 
3033    // iffreq: list of minimum and maximum frequency for each FREQ_ID in
3034    //         each table
3035    // freqid[insize][numIF]
3036    // freqid: [[id00, id01, ...],
3037    //          [id10, id11, ...],
3038    //          ...
3039    //          [idn0, idn1, ...]]
3040    // iffreq[insize][numIF*2]
3041    // iffreq: [[min_id00, max_id00, min_id01, max_id01, ...],
3042    //          [min_id10, max_id10, min_id11, max_id11, ...],
3043    //          ...
3044    //          [min_idn0, max_idn0, min_idn1, max_idn1, ...]]
3045    //os << "Check IF settings in each table" << LogIO::POST ;
3046    vector< vector<uInt> > freqid( insize );
3047    vector< vector<double> > iffreq( insize ) ;
3048    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3049      uInt rows = tmpin[itable]->nrow() ;
3050      uInt freqnrows = tmpin[itable]->frequencies().table().nrow() ;
3051      for ( uInt irow = 0 ; irow < rows ; irow++ ) {
3052        if ( freqid[itable].size() == freqnrows ) {
3053          break ;
3054        }
3055        else {
3056          freqIDCol.attach( tmpin[itable]->table(), "FREQ_ID" ) ;
3057          ifnoCol.attach( tmpin[itable]->table(), "IFNO" ) ;
3058          uInt id = freqIDCol( irow ) ;
3059          if ( freqid[itable].size() == 0 || count( freqid[itable].begin(), freqid[itable].end(), id ) == 0 ) {
3060            //os << "itable = " << itable << ": IF " << id << " is included in the list" << LogIO::POST ;
3061            vector<double> abcissa = tmpin[itable]->getAbcissa( irow ) ;
3062            freqid[itable].push_back( id ) ;
3063            iffreq[itable].push_back( abcissa[0] - 0.5 * ( abcissa[1] - abcissa[0] ) ) ;
3064            iffreq[itable].push_back( abcissa[abcissa.size()-1] + 0.5 * ( abcissa[1] - abcissa[0] ) ) ;
3065          }
3066        }
3067      }
3068    }
3069
3070    // debug
3071    //os << "IF settings summary:" << endl ;
3072    //for ( uInt i = 0 ; i < freqid.size() ; i++ ) {
3073    //os << "   Table" << i << endl ;
3074    //for ( uInt j = 0 ; j < freqid[i].size() ; j++ ) {
3075    //os << "      id = " << freqid[i][j] << " (min,max) = (" << iffreq[i][2*j] << "," << iffreq[i][2*j+1] << ")" << endl ;
3076    //}
3077    //}
3078    //os << endl ;
3079    //os.post() ;
3080
3081    // IF grouping based on their frequency coverage
3082    // ifgrp: list of table index and FREQ_ID for all members in each IF group
3083    // ifgfreq: list of minimum and maximum frequency in each IF group
3084    // ifgrp[numgrp][nummember*2]
3085    // ifgrp: [[table00, freqrow00, table01, freqrow01, ...],
3086    //         [table10, freqrow10, table11, freqrow11, ...],
3087    //         ...
3088    //         [tablen0, freqrown0, tablen1, freqrown1, ...]]
3089    // ifgfreq[numgrp*2]
3090    // ifgfreq: [min0_grp0, max0_grp0, min1_grp1, max1_grp1, ...]
3091    //os << "IF grouping based on their frequency coverage" << LogIO::POST ;
3092    vector< vector<uInt> > ifgrp ;
3093    vector<double> ifgfreq ;
3094
3095    // parameter for IF grouping
3096    // groupmode = OR    retrieve all region
3097    //             AND   only retrieve overlaped region
3098    //string groupmode = "AND" ;
3099    string groupmode = "OR" ;
3100    uInt sizecr = 0 ;
3101    if ( groupmode == "AND" )
3102      sizecr = 2 ;
3103    else if ( groupmode == "OR" )
3104      sizecr = 0 ;
3105
3106    vector<double> sortedfreq ;
3107    for ( uInt i = 0 ; i < iffreq.size() ; i++ ) {
3108      for ( uInt j = 0 ; j < iffreq[i].size() ; j++ ) {
3109        if ( count( sortedfreq.begin(), sortedfreq.end(), iffreq[i][j] ) == 0 )
3110          sortedfreq.push_back( iffreq[i][j] ) ;
3111      }
3112    }
3113    sort( sortedfreq.begin(), sortedfreq.end() ) ;
3114    for ( vector<double>::iterator i = sortedfreq.begin() ; i != sortedfreq.end()-1 ; i++ ) {
3115      ifgfreq.push_back( *i ) ;
3116      ifgfreq.push_back( *(i+1) ) ;
3117    }
3118    ifgrp.resize( ifgfreq.size()/2 ) ;
3119    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3120      for ( uInt iif = 0 ; iif < freqid[itable].size() ; iif++ ) {
3121        double range0 = iffreq[itable][2*iif] ;
3122        double range1 = iffreq[itable][2*iif+1] ;
3123        for ( uInt j = 0 ; j < ifgrp.size() ; j++ ) {
3124          double fmin = max( range0, ifgfreq[2*j] ) ;
3125          double fmax = min( range1, ifgfreq[2*j+1] ) ;
3126          if ( fmin < fmax ) {
3127            ifgrp[j].push_back( itable ) ;
3128            ifgrp[j].push_back( freqid[itable][iif] ) ;
3129          }
3130        }
3131      }
3132    }
3133    vector< vector<uInt> >::iterator fiter = ifgrp.begin() ;
3134    vector<double>::iterator giter = ifgfreq.begin() ;
3135    while( fiter != ifgrp.end() ) {
3136      if ( fiter->size() <= sizecr ) {
3137        fiter = ifgrp.erase( fiter ) ;
3138        giter = ifgfreq.erase( giter ) ;
3139        giter = ifgfreq.erase( giter ) ;
3140      }
3141      else {
3142        fiter++ ;
3143        advance( giter, 2 ) ;
3144      }
3145    }
3146
3147    // Grouping continuous IF groups (without frequency gap)
3148    // freqgrp: list of IF group indexes in each frequency group
3149    // freqrange: list of minimum and maximum frequency in each frequency group
3150    // freqgrp[numgrp][nummember]
3151    // freqgrp: [[ifgrp00, ifgrp01, ifgrp02, ...],
3152    //           [ifgrp10, ifgrp11, ifgrp12, ...],
3153    //           ...
3154    //           [ifgrpn0, ifgrpn1, ifgrpn2, ...]]
3155    // freqrange[numgrp*2]
3156    // freqrange: [min_grp0, max_grp0, min_grp1, max_grp1, ...]
3157    vector< vector<uInt> > freqgrp ;
3158    double freqrange = 0.0 ;
3159    uInt grpnum = 0 ;
3160    for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
3161      // Assumed that ifgfreq was sorted
3162      if ( grpnum != 0 && freqrange == ifgfreq[2*i] ) {
3163        freqgrp[grpnum-1].push_back( i ) ;
3164      }
3165      else {
3166        vector<uInt> grp0( 1, i ) ;
3167        freqgrp.push_back( grp0 ) ;
3168        grpnum++ ;
3169      }
3170      freqrange = ifgfreq[2*i+1] ;
3171    }
3172       
3173
3174    // print IF groups
3175    ostringstream oss ;
3176    oss << "IF Group summary: " << endl ;
3177    oss << "   GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
3178    for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
3179      oss << "   GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
3180      for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
3181        oss << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ;
3182      }
3183      oss << endl ;
3184    }
3185    oss << endl ;
3186    os << oss.str() << LogIO::POST ;
3187   
3188    // print frequency group
3189    oss.str("") ;
3190    oss << "Frequency Group summary: " << endl ;
3191    oss << "   GROUP_ID [FREQ_MIN, FREQ_MAX]: IF_GROUP_ID" << endl ;
3192    for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
3193      oss << "   GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*freqgrp[i][0]] << "," << ifgfreq[2*freqgrp[i][freqgrp[i].size()-1]+1] << "]: " ;
3194      for ( uInt j = 0 ; j < freqgrp[i].size() ; j++ ) {
3195        oss << freqgrp[i][j] << " " ;
3196      }
3197      oss << endl ;
3198    }
3199    oss << endl ;
3200    os << oss.str() << LogIO::POST ;
3201
3202    // membership check
3203    // groups: list of IF group indexes whose frequency range overlaps with
3204    //         that of each table and IF
3205    // groups[numtable][numIF][nummembership]
3206    // groups: [[[grp, grp,...], [grp, grp,...],...],
3207    //          [[grp, grp,...], [grp, grp,...],...],
3208    //          ...
3209    //          [[grp, grp,...], [grp, grp,...],...]]
3210    vector< vector< vector<uInt> > > groups( insize ) ;
3211    for ( uInt i = 0 ; i < insize ; i++ ) {
3212      groups[i].resize( freqid[i].size() ) ;
3213    }
3214    for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
3215      for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
3216        uInt tableid = ifgrp[igrp][2*imem] ;
3217        vector<uInt>::iterator iter = find( freqid[tableid].begin(), freqid[tableid].end(), ifgrp[igrp][2*imem+1] ) ;
3218        if ( iter != freqid[tableid].end() ) {
3219          uInt rowid = distance( freqid[tableid].begin(), iter ) ;
3220          groups[tableid][rowid].push_back( igrp ) ;
3221        }
3222      }
3223    }
3224
3225    // print membership
3226    //oss.str("") ;
3227    //for ( uInt i = 0 ; i < insize ; i++ ) {
3228    //oss << "Table " << i << endl ;
3229    //for ( uInt j = 0 ; j < groups[i].size() ; j++ ) {
3230    //oss << "   FREQ_ID " <<  setw( 2 ) << freqid[i][j] << ": " ;
3231    //for ( uInt k = 0 ; k < groups[i][j].size() ; k++ ) {
3232    //oss << setw( 2 ) << groups[i][j][k] << " " ;
3233    //}
3234    //oss << endl ;
3235    //}
3236    //}
3237    //os << oss.str() << LogIO::POST ;
3238
3239    // set back coordinfo
3240    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3241      vector<string> coordinfo = tmpin[itable]->getCoordInfo() ;
3242      coordinfo[0] = oldinfo[itable] ;
3243      tmpin[itable]->setCoordInfo( coordinfo ) ;
3244    }
3245
3246    // Create additional table if needed
3247    bool oldInsitu = insitu_ ;
3248    setInsitu( false ) ;
3249    vector< vector<uInt> > addrow( insize ) ;
3250    vector<uInt> addtable( insize, 0 ) ;
3251    vector<uInt> newtableids( insize ) ;
3252    vector<uInt> newifids( insize, 0 ) ;
3253    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3254      //os << "Table " << itable << ": " ;
3255      for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
3256        addrow[itable].push_back( groups[itable][ifrow].size()-1 ) ;
3257        //os << addrow[itable][ifrow] << " " ;
3258      }
3259      addtable[itable] = *max_element( addrow[itable].begin(), addrow[itable].end() ) ;
3260      //os << "(" << addtable[itable] << ")" << LogIO::POST ;
3261    }
3262    newin.resize( insize ) ;
3263    copy( tmpin.begin(), tmpin.end(), newin.begin() ) ;
3264    for ( uInt i = 0 ; i < insize ; i++ ) {
3265      newtableids[i] = i ;
3266    }
3267    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3268      for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
3269        CountedPtr<Scantable> add = getScantable( newin[itable], false ) ;
3270        vector<int> freqidlist ;
3271        for ( uInt i = 0 ; i < groups[itable].size() ; i++ ) {
3272          if ( groups[itable][i].size() > iadd + 1 ) {
3273            freqidlist.push_back( freqid[itable][i] ) ;
3274          }
3275        }
3276        stringstream taqlstream ;
3277        taqlstream << "SELECT FROM $1 WHERE FREQ_ID IN [" ;
3278        for ( uInt i = 0 ; i < freqidlist.size() ; i++ ) {
3279          taqlstream << freqidlist[i] ;
3280          if ( i < freqidlist.size() - 1 )
3281            taqlstream << "," ;
3282          else
3283            taqlstream << "]" ;
3284        }
3285        string taql = taqlstream.str() ;
3286        //os << "taql = " << taql << LogIO::POST ;
3287        STSelector selector = STSelector() ;
3288        selector.setTaQL( taql ) ;
3289        add->setSelection( selector ) ;
3290        newin.push_back( add ) ;
3291        newtableids.push_back( itable ) ;
3292        newifids.push_back( iadd + 1 ) ;
3293      }
3294    }
3295
3296    // udpate ifgrp
3297    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3298      for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
3299        for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
3300          if ( groups[itable][ifrow].size() > iadd + 1 ) {
3301            uInt igrp = groups[itable][ifrow][iadd+1] ;
3302            for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
3303              if ( ifgrp[igrp][2*imem] == newtableids[iadd+insize] && ifgrp[igrp][2*imem+1] == freqid[newtableids[iadd+insize]][ifrow] ) {
3304                ifgrp[igrp][2*imem] = insize + iadd ;
3305              }
3306            }
3307          }
3308        }
3309      }
3310    }
3311
3312    // print IF groups again for debug
3313    //oss.str( "" ) ;
3314    //oss << "IF Group summary: " << endl ;
3315    //oss << "   GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
3316    //for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
3317    //oss << "   GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
3318    //for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
3319    //oss << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ;
3320    //}
3321    //oss << endl ;
3322    //}
3323    //oss << endl ;
3324    //os << oss.str() << LogIO::POST ;
3325
3326    // reset SCANNO and IFNO/FREQ_ID: IF is reset by the result of sortation
3327    os << "All scan number is set to 0" << LogIO::POST ;
3328    //os << "All IF number is set to IF group index" << LogIO::POST ;
3329    insize = newin.size() ;
3330    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3331      uInt rows = newin[itable]->nrow() ;
3332      Table &tmpt = newin[itable]->table() ;
3333      freqIDCol.attach( tmpt, "FREQ_ID" ) ;
3334      scannoCol.attach( tmpt, "SCANNO" ) ;
3335      ifnoCol.attach( tmpt, "IFNO" ) ;
3336      for ( uInt irow=0 ; irow < rows ; irow++ ) {
3337        scannoCol.put( irow, 0 ) ;
3338        uInt freqID = freqIDCol( irow ) ;
3339        vector<uInt>::iterator iter = find( freqid[newtableids[itable]].begin(), freqid[newtableids[itable]].end(), freqID ) ;
3340        if ( iter != freqid[newtableids[itable]].end() ) {
3341          uInt index = distance( freqid[newtableids[itable]].begin(), iter ) ;
3342          ifnoCol.put( irow, groups[newtableids[itable]][index][newifids[itable]] ) ;
3343        }
3344        else {
3345          throw(AipsError("IF grouping was wrong in additional tables.")) ;
3346        }
3347      }
3348    }
3349    oldinfo.resize( insize ) ;
3350    setInsitu( oldInsitu ) ;
3351
3352    // temporarily set coordinfo
3353    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3354      vector<string> coordinfo = newin[itable]->getCoordInfo() ;
3355      oldinfo[itable] = coordinfo[0] ;
3356      coordinfo[0] = "Hz" ;
3357      newin[itable]->setCoordInfo( coordinfo ) ;
3358    }
3359
3360    // save column values in the vector
3361    vector< vector<uInt> > freqTableIdVec( insize ) ;
3362    vector< vector<uInt> > freqIdVec( insize ) ;
3363    vector< vector<uInt> > ifNoVec( insize ) ;
3364    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3365      ScalarColumn<uInt> freqIDs ;
3366      freqIDs.attach( newin[itable]->frequencies().table(), "ID" ) ;
3367      ifnoCol.attach( newin[itable]->table(), "IFNO" ) ;
3368      freqIDCol.attach( newin[itable]->table(), "FREQ_ID" ) ;
3369      for ( uInt irow = 0 ; irow < newin[itable]->frequencies().table().nrow() ; irow++ ) {
3370        freqTableIdVec[itable].push_back( freqIDs( irow ) ) ;
3371      }
3372      for ( uInt irow = 0 ; irow < newin[itable]->table().nrow() ; irow++ ) {
3373        freqIdVec[itable].push_back( freqIDCol( irow ) ) ;
3374        ifNoVec[itable].push_back( ifnoCol( irow ) ) ;
3375      }
3376    }
3377
3378    // reset spectra and flagtra: pick up common part of frequency coverage
3379    //os << "Pick common frequency range and align resolution" << LogIO::POST ;
3380    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3381      uInt rows = newin[itable]->nrow() ;
3382      int nminchan = -1 ;
3383      int nmaxchan = -1 ;
3384      vector<uInt> freqIdUpdate ;
3385      for ( uInt irow = 0 ; irow < rows ; irow++ ) {
3386        uInt ifno = ifNoVec[itable][irow] ;  // IFNO is reset by group index
3387        double minfreq = ifgfreq[2*ifno] ;
3388        double maxfreq = ifgfreq[2*ifno+1] ;
3389        //os << "frequency range: [" << minfreq << "," << maxfreq << "]" << LogIO::POST ;
3390        vector<double> abcissa = newin[itable]->getAbcissa( irow ) ;
3391        int nchan = abcissa.size() ;
3392        double resol = abcissa[1] - abcissa[0] ;
3393        //os << "abcissa range  : [" << abcissa[0] << "," << abcissa[nchan-1] << "]" << LogIO::POST ;
3394        if ( minfreq <= abcissa[0] )
3395          nminchan = 0 ;
3396        else {
3397          //double cfreq = ( minfreq - abcissa[0] ) / resol ;
3398          double cfreq = ( minfreq - abcissa[0] + 0.5 * resol ) / resol ;
3399          nminchan = int(cfreq) + ( ( cfreq - int(cfreq) <= 0.5 ) ? 0 : 1 ) ;
3400        }
3401        if ( maxfreq >= abcissa[abcissa.size()-1] )
3402          nmaxchan = abcissa.size() - 1 ;
3403        else {
3404          //double cfreq = ( abcissa[abcissa.size()-1] - maxfreq ) / resol ;
3405          double cfreq = ( abcissa[abcissa.size()-1] - maxfreq + 0.5 * resol ) / resol ;
3406          nmaxchan = abcissa.size() - 1 - int(cfreq) - ( ( cfreq - int(cfreq) >= 0.5 ) ? 1 : 0 ) ;
3407        }
3408        //os << "channel range (" << irow << "): [" << nminchan << "," << nmaxchan << "]" << LogIO::POST ;
3409        if ( nmaxchan > nminchan ) {
3410          newin[itable]->reshapeSpectrum( nminchan, nmaxchan, irow ) ;
3411          int newchan = nmaxchan - nminchan + 1 ;
3412          if ( count( freqIdUpdate.begin(), freqIdUpdate.end(), freqIdVec[itable][irow] ) == 0 && newchan < nchan )
3413            freqIdUpdate.push_back( freqIdVec[itable][irow] ) ;
3414        }
3415        else {
3416          throw(AipsError("Failed to pick up common part of frequency range.")) ;
3417        }
3418      }
3419      for ( uInt i = 0 ; i < freqIdUpdate.size() ; i++ ) {
3420        uInt freqId = freqIdUpdate[i] ;
3421        Double refpix ;
3422        Double refval ;
3423        Double increment ;
3424       
3425        // update row
3426        newin[itable]->frequencies().getEntry( refpix, refval, increment, freqId ) ;
3427        refval = refval - ( refpix - nminchan ) * increment ;
3428        refpix = 0 ;
3429        newin[itable]->frequencies().setEntry( refpix, refval, increment, freqId ) ;
3430      }   
3431    }
3432
3433   
3434    // reset spectra and flagtra: align spectral resolution
3435    //os << "Align spectral resolution" << LogIO::POST ;
3436    // gmaxdnu: the coarsest frequency resolution in the frequency group
3437    // gmemid: member index that have a resolution equal to gmaxdnu
3438    // gmaxdnu[numfreqgrp]
3439    // gmaxdnu: [dnu0, dnu1, ...]
3440    // gmemid[numfreqgrp]
3441    // gmemid: [id0, id1, ...]
3442    vector<double> gmaxdnu( freqgrp.size(), 0.0 ) ;
3443    vector<uInt> gmemid( freqgrp.size(), 0 ) ;
3444    for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
3445      double maxdnu = 0.0 ;       // maximum (coarsest) frequency resolution
3446      int minchan = INT_MAX ;     // minimum channel number
3447      Double refpixref = -1 ;     // reference of 'reference pixel'
3448      Double refvalref = -1 ;     // reference of 'reference frequency'
3449      Double refinc = -1 ;        // reference frequency resolution
3450      uInt refreqid ;
3451      uInt reftable = INT_MAX;
3452      // process only if group member > 1
3453      if ( ifgrp[igrp].size() > 2 ) {
3454        // find minchan and maxdnu in each group
3455        for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
3456          uInt tableid = ifgrp[igrp][2*imem] ;
3457          uInt rowid = ifgrp[igrp][2*imem+1] ;
3458          vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
3459          if ( iter != freqIdVec[tableid].end() ) {
3460            uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
3461            vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
3462            int nchan = abcissa.size() ;
3463            double dnu = abcissa[1] - abcissa[0] ;
3464            //os << "GROUP " << igrp << " (" << tableid << "," << rowid << "): nchan = " << nchan << " (minchan = " << minchan << ")" << LogIO::POST ;
3465            if ( nchan < minchan ) {
3466              minchan = nchan ;
3467              maxdnu = dnu ;
3468              newin[tableid]->frequencies().getEntry( refpixref, refvalref, refinc, rowid ) ;
3469              refreqid = rowid ;
3470              reftable = tableid ;
3471            }
3472          }
3473        }
3474        // regrid spectra in each group
3475        os << "GROUP " << igrp << endl ;
3476        os << "   Channel number is adjusted to " << minchan << endl ;
3477        os << "   Corresponding frequency resolution is " << maxdnu << "Hz" << LogIO::POST ;
3478        for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
3479          uInt tableid = ifgrp[igrp][2*imem] ;
3480          uInt rowid = ifgrp[igrp][2*imem+1] ;
3481          freqIDCol.attach( newin[tableid]->table(), "FREQ_ID" ) ;
3482          //os << "tableid = " << tableid << " rowid = " << rowid << ": " << LogIO::POST ;
3483          //os << "   regridChannel applied to " ;
3484          //if ( tableid != reftable )
3485          refreqid = newin[tableid]->frequencies().addEntry( refpixref, refvalref, refinc ) ;
3486          for ( uInt irow = 0 ; irow < newin[tableid]->table().nrow() ; irow++ ) {
3487            uInt tfreqid = freqIdVec[tableid][irow] ;
3488            if ( tfreqid == rowid ) {     
3489              //os << irow << " " ;
3490              newin[tableid]->regridChannel( minchan, maxdnu, irow ) ;
3491              freqIDCol.put( irow, refreqid ) ;
3492              freqIdVec[tableid][irow] = refreqid ;
3493            }
3494          }
3495          //os << LogIO::POST ;
3496        }
3497      }
3498      else {
3499        uInt tableid = ifgrp[igrp][0] ;
3500        uInt rowid = ifgrp[igrp][1] ;
3501        vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
3502        if ( iter != freqIdVec[tableid].end() ) {
3503          uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
3504          vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
3505          minchan = abcissa.size() ;
3506          maxdnu = abcissa[1] - abcissa[0] ;
3507        }
3508      }
3509      for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
3510        if ( count( freqgrp[i].begin(), freqgrp[i].end(), igrp ) > 0 ) {
3511          if ( maxdnu > gmaxdnu[i] ) {
3512            gmaxdnu[i] = maxdnu ;
3513            gmemid[i] = igrp ;
3514          }
3515          break ;
3516        }
3517      }
3518    }
3519
3520    // set back coordinfo
3521    for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3522      vector<string> coordinfo = newin[itable]->getCoordInfo() ;
3523      coordinfo[0] = oldinfo[itable] ;
3524      newin[itable]->setCoordInfo( coordinfo ) ;
3525    }     
3526
3527    // accumulate all rows into the first table
3528    // NOTE: assumed in.size() = 1
3529    vector< CountedPtr<Scantable> > tmp( 1 ) ;
3530    if ( newin.size() == 1 )
3531      tmp[0] = newin[0] ;
3532    else
3533      tmp[0] = merge( newin ) ;
3534
3535    //return tmp[0] ;
3536
3537    // average
3538    CountedPtr<Scantable> tmpout = average( tmp, mask, weight, avmode ) ;
3539
3540    //return tmpout ;
3541
3542    // combine frequency group
3543    os << "Combine spectra based on frequency grouping" << LogIO::POST ;
3544    os << "IFNO is renumbered as frequency group ID (see above)" << LogIO::POST ;
3545    vector<string> coordinfo = tmpout->getCoordInfo() ;
3546    oldinfo[0] = coordinfo[0] ;
3547    coordinfo[0] = "Hz" ;
3548    tmpout->setCoordInfo( coordinfo ) ;
3549    // create proformas of output table
3550    stringstream taqlstream ;
3551    taqlstream << "SELECT FROM $1 WHERE IFNO IN [" ;
3552    for ( uInt i = 0 ; i < gmemid.size() ; i++ ) {
3553      taqlstream << gmemid[i] ;
3554      if ( i < gmemid.size() - 1 )
3555        taqlstream << "," ;
3556      else
3557        taqlstream << "]" ;
3558    }
3559    string taql = taqlstream.str() ;
3560    //os << "taql = " << taql << LogIO::POST ;
3561    STSelector selector = STSelector() ;
3562    selector.setTaQL( taql ) ;
3563    oldInsitu = insitu_ ;
3564    setInsitu( false ) ;
3565    out = getScantable( tmpout, false ) ;
3566    setInsitu( oldInsitu ) ;
3567    out->setSelection( selector ) ;
3568    // regrid rows
3569    ifnoCol.attach( tmpout->table(), "IFNO" ) ;
3570    for ( uInt irow = 0 ; irow < tmpout->table().nrow() ; irow++ ) {
3571      uInt ifno = ifnoCol( irow ) ;
3572      for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3573        if ( count( freqgrp[igrp].begin(), freqgrp[igrp].end(), ifno ) > 0 ) {
3574          vector<double> abcissa = tmpout->getAbcissa( irow ) ;
3575          double bw = ( abcissa[1] - abcissa[0] ) * abcissa.size() ;
3576          int nchan = (int)( bw / gmaxdnu[igrp] ) ;
3577          tmpout->regridChannel( nchan, gmaxdnu[igrp], irow ) ;
3578          break ;
3579        }
3580      }
3581    }
3582    // combine spectra
3583    ArrayColumn<Float> specColOut ;
3584    specColOut.attach( out->table(), "SPECTRA" ) ;
3585    ArrayColumn<uChar> flagColOut ;
3586    flagColOut.attach( out->table(), "FLAGTRA" ) ;
3587    ScalarColumn<uInt> ifnoColOut ;
3588    ifnoColOut.attach( out->table(), "IFNO" ) ;
3589    ScalarColumn<uInt> polnoColOut ;
3590    polnoColOut.attach( out->table(), "POLNO" ) ;
3591    ScalarColumn<uInt> freqidColOut ;
3592    freqidColOut.attach( out->table(), "FREQ_ID" ) ;
3593//     MDirection::ScalarColumn dirColOut ;
3594//     dirColOut.attach( out->table(), "DIRECTION" ) ;
3595    Table &tab = tmpout->table() ;
3596    Block<String> cols(1);
3597    cols[0] = String("POLNO") ;
3598    TableIterator iter( tab, cols ) ;
3599    bool done = false ;
3600    vector< vector<uInt> > sizes( freqgrp.size() ) ;
3601    while( !iter.pastEnd() ) {
3602      vector< vector<Float> > specout( freqgrp.size() ) ;
3603      vector< vector<uChar> > flagout( freqgrp.size() ) ;
3604      ArrayColumn<Float> specCols ;
3605      specCols.attach( iter.table(), "SPECTRA" ) ;
3606      ArrayColumn<uChar> flagCols ;
3607      flagCols.attach( iter.table(), "FLAGTRA" ) ;
3608      ifnoCol.attach( iter.table(), "IFNO" ) ;
3609      ScalarColumn<uInt> polnos ;
3610      polnos.attach( iter.table(), "POLNO" ) ;
3611//       MDirection::ScalarColumn dircol ;
3612//       dircol.attach( iter.table(), "DIRECTION" ) ;
3613      uInt polno = polnos( 0 ) ;
3614      //os << "POLNO iteration: " << polno << LogIO::POST ;
3615//       for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3616//      sizes[igrp].resize( freqgrp[igrp].size() ) ;
3617//      for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
3618//        for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
3619//          uInt ifno = ifnoCol( irow ) ;
3620//          if ( ifno == freqgrp[igrp][imem] ) {
3621//            Vector<Float> spec = specCols( irow ) ;
3622//            Vector<uChar> flag = flagCols( irow ) ;
3623//            vector<Float> svec ;
3624//            spec.tovector( svec ) ;
3625//            vector<uChar> fvec ;
3626//            flag.tovector( fvec ) ;
3627//            //os << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << LogIO::POST ;
3628//            specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ;
3629//            flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ;
3630//            //os << "specout[" << igrp << "].size() = " << specout[igrp].size() << LogIO::POST ;
3631//            sizes[igrp][imem] = spec.nelements() ;
3632//          }
3633//        }
3634//      }
3635//      for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
3636//        uInt ifout = ifnoColOut( irow ) ;
3637//        uInt polout = polnoColOut( irow ) ;
3638//        if ( ifout == gmemid[igrp] && polout == polno ) {
3639//          // set SPECTRA and FRAGTRA
3640//          Vector<Float> newspec( specout[igrp] ) ;
3641//          Vector<uChar> newflag( flagout[igrp] ) ;
3642//          specColOut.put( irow, newspec ) ;
3643//          flagColOut.put( irow, newflag ) ;
3644//          // IFNO renumbering
3645//          ifnoColOut.put( irow, igrp ) ;
3646//        }
3647//      }
3648//       }
3649      // get a list of number of channels for each frequency group member
3650      if ( !done ) {
3651        for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3652          sizes[igrp].resize( freqgrp[igrp].size() ) ;
3653          for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
3654            for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
3655              uInt ifno = ifnoCol( irow ) ;
3656              if ( ifno == freqgrp[igrp][imem] ) {
3657                Vector<Float> spec = specCols( irow ) ;
3658                sizes[igrp][imem] = spec.nelements() ;
3659                break ;
3660              }               
3661            }
3662          }
3663        }
3664        done = true ;
3665      }
3666      // combine spectra
3667      for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
3668        uInt polout = polnoColOut( irow ) ;
3669        if ( polout == polno ) {
3670          uInt ifout = ifnoColOut( irow ) ;
3671//           Vector<Double> direction = dirColOut(irow).getAngle(Unit(String("rad"))).getValue() ;
3672          uInt igrp ;
3673          for ( uInt jgrp = 0 ; jgrp < freqgrp.size() ; jgrp++ ) {
3674            if ( ifout == gmemid[jgrp] ) {
3675              igrp = jgrp ;
3676              break ;
3677            }
3678          }
3679          for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
3680            for ( uInt jrow = 0 ; jrow < iter.table().nrow() ; jrow++ ) {
3681              uInt ifno = ifnoCol( jrow ) ;
3682              // 2012/02/17 TN
3683              // Since STGrid is implemented, average doesn't consider direction
3684              // when accumulating
3685//               Vector<Double> tdir = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
3686//               //if ( ifno == freqgrp[igrp][imem] && allTrue( tdir == direction  ) ) {
3687//               Double dx = tdir[0] - direction[0] ;
3688//               Double dy = tdir[1] - direction[1] ;
3689//               Double dd = sqrt( dx * dx + dy * dy ) ;
3690              //if ( ifno == freqgrp[igrp][imem] && allNearAbs( tdir, direction, tol ) ) {
3691//               if ( ifno == freqgrp[igrp][imem] && dd <= tol ) {
3692              if ( ifno == freqgrp[igrp][imem] ) {
3693                Vector<Float> spec = specCols( jrow ) ;
3694                Vector<uChar> flag = flagCols( jrow ) ;
3695                vector<Float> svec ;
3696                spec.tovector( svec ) ;
3697                vector<uChar> fvec ;
3698                flag.tovector( fvec ) ;
3699                //os << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << LogIO::POST ;
3700                specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ;
3701                flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ;
3702                //os << "specout[" << igrp << "].size() = " << specout[igrp].size() << LogIO::POST ;
3703              }
3704            }
3705          }
3706          // set SPECTRA and FRAGTRA
3707          Vector<Float> newspec( specout[igrp] ) ;
3708          Vector<uChar> newflag( flagout[igrp] ) ;
3709          specColOut.put( irow, newspec ) ;
3710          flagColOut.put( irow, newflag ) ;
3711          // IFNO renumbering
3712          ifnoColOut.put( irow, igrp ) ;
3713        }
3714      }
3715      iter++ ;
3716    }
3717    // update FREQUENCIES subtable
3718    vector<bool> updated( freqgrp.size(), false ) ;
3719    for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3720      uInt index = 0 ;
3721      uInt pixShift = 0 ;
3722      while ( freqgrp[igrp][index] != gmemid[igrp] ) {
3723        pixShift += sizes[igrp][index++] ;
3724      }
3725      for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
3726        if ( ifnoColOut( irow ) == gmemid[igrp] && !updated[igrp] ) {
3727          uInt freqidOut = freqidColOut( irow ) ;
3728          //os << "freqgrp " << igrp << " freqidOut = " << freqidOut << LogIO::POST ;
3729          double refpix ;
3730          double refval ;
3731          double increm ;
3732          out->frequencies().getEntry( refpix, refval, increm, freqidOut ) ;
3733          refpix += pixShift ;
3734          out->frequencies().setEntry( refpix, refval, increm, freqidOut ) ;
3735          updated[igrp] = true ;
3736        }
3737      }
3738    }
3739
3740    //out = tmpout ;
3741
3742    coordinfo = tmpout->getCoordInfo() ;
3743    coordinfo[0] = oldinfo[0] ;
3744    tmpout->setCoordInfo( coordinfo ) ;
3745  }
3746  else {
3747    // simple average
3748    out =  average( in, mask, weight, avmode ) ;
3749  }
3750 
3751  return out;
3752}
3753
3754CountedPtr<Scantable> STMath::cwcal( const CountedPtr<Scantable>& s,
3755                                     const String calmode,
3756                                     const String antname )
3757{
3758  // frequency switch
3759  if ( calmode == "fs" ) {
3760    return cwcalfs( s, antname ) ;
3761  }
3762  else {
3763    vector<bool> masks = s->getMask( 0 ) ;
3764    vector<int> types ;
3765
3766    // sky scan
3767    STSelector sel = STSelector() ;
3768    types.push_back( SrcType::SKY ) ;
3769    sel.setTypes( types ) ;
3770    s->setSelection( sel ) ;
3771    vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ;
3772    CountedPtr<Scantable> asky = average( tmp, masks, "TINT", "SCAN" ) ;
3773    s->unsetSelection() ;
3774    sel.reset() ;
3775    types.clear() ;
3776
3777    // hot scan
3778    types.push_back( SrcType::HOT ) ;
3779    sel.setTypes( types ) ;
3780    s->setSelection( sel ) ;
3781    tmp.clear() ;
3782    tmp.push_back( getScantable( s, false ) ) ;
3783    CountedPtr<Scantable> ahot = average( tmp, masks, "TINT", "SCAN" ) ;
3784    s->unsetSelection() ;
3785    sel.reset() ;
3786    types.clear() ;
3787   
3788    // cold scan
3789    CountedPtr<Scantable> acold ;
3790//     types.push_back( SrcType::COLD ) ;
3791//     sel.setTypes( types ) ;
3792//     s->setSelection( sel ) ;
3793//     tmp.clear() ;
3794//     tmp.push_back( getScantable( s, false ) ) ;
3795//     CountedPtr<Scantable> acold = average( tmp, masks, "TINT", "SCNAN" ) ;
3796//     s->unsetSelection() ;
3797//     sel.reset() ;
3798//     types.clear() ;
3799
3800    // off scan
3801    types.push_back( SrcType::PSOFF ) ;
3802    sel.setTypes( types ) ;
3803    s->setSelection( sel ) ;
3804    tmp.clear() ;
3805    tmp.push_back( getScantable( s, false ) ) ;
3806    CountedPtr<Scantable> aoff = average( tmp, masks, "TINT", "SCAN" ) ;
3807    s->unsetSelection() ;
3808    sel.reset() ;
3809    types.clear() ;
3810   
3811    // on scan
3812    bool insitu = insitu_ ;
3813    insitu_ = false ;
3814    CountedPtr<Scantable> out = getScantable( s, true ) ;
3815    insitu_ = insitu ;
3816    types.push_back( SrcType::PSON ) ;
3817    sel.setTypes( types ) ;
3818    s->setSelection( sel ) ;
3819    TableCopy::copyRows( out->table(), s->table() ) ;
3820    s->unsetSelection() ;
3821    sel.reset() ;
3822    types.clear() ;
3823   
3824    // process each on scan
3825    ArrayColumn<Float> tsysCol ;
3826    tsysCol.attach( out->table(), "TSYS" ) ;
3827    for ( int i = 0 ; i < out->nrow() ; i++ ) {
3828      vector<float> sp = getCalibratedSpectra( out, aoff, asky, ahot, acold, i, antname ) ;
3829      out->setSpectrum( sp, i ) ;
3830      string reftime = out->getTime( i ) ;
3831      vector<int> ii( 1, out->getIF( i ) ) ;
3832      vector<int> ib( 1, out->getBeam( i ) ) ;
3833      vector<int> ip( 1, out->getPol( i ) ) ;
3834      sel.setIFs( ii ) ;
3835      sel.setBeams( ib ) ;
3836      sel.setPolarizations( ip ) ;
3837      asky->setSelection( sel ) ;   
3838      vector<float> sptsys = getTsysFromTime( reftime, asky, "linear" ) ;
3839      const Vector<Float> Vtsys( sptsys ) ;
3840      tsysCol.put( i, Vtsys ) ;
3841      asky->unsetSelection() ;
3842      sel.reset() ;
3843    }
3844
3845    // flux unit
3846    out->setFluxUnit( "K" ) ;
3847
3848    return out ;
3849  }
3850}
3851 
3852CountedPtr<Scantable> STMath::almacal( const CountedPtr<Scantable>& s,
3853                                       const String calmode )
3854{
3855  // frequency switch
3856  if ( calmode == "fs" ) {
3857    return almacalfs( s ) ;
3858  }
3859  else {
3860    vector<bool> masks = s->getMask( 0 ) ;
3861   
3862    // off scan
3863    STSelector sel = STSelector() ;
3864    vector<int> types ;
3865    types.push_back( SrcType::PSOFF ) ;
3866    sel.setTypes( types ) ;
3867    s->setSelection( sel ) ;
3868    // TODO 2010/01/08 TN
3869    // Grouping by time should be needed before averaging.
3870    // Each group must have own unique SCANNO (should be renumbered).
3871    // See PIPELINE/SDCalibration.py
3872    CountedPtr<Scantable> soff = getScantable( s, false ) ;
3873    Table ttab = soff->table() ;
3874    ROScalarColumn<Double> timeCol( ttab, "TIME" ) ;
3875    uInt nrow = timeCol.nrow() ;
3876    Vector<Double> timeSep( nrow - 1 ) ;
3877    for ( uInt i = 0 ; i < nrow - 1 ; i++ ) {
3878      timeSep[i] = timeCol(i+1) - timeCol(i) ;
3879    }
3880    ScalarColumn<Double> intervalCol( ttab, "INTERVAL" ) ;
3881    Vector<Double> interval = intervalCol.getColumn() ;
3882    interval /= 86400.0 ;
3883    ScalarColumn<uInt> scanCol( ttab, "SCANNO" ) ;
3884    vector<uInt> glist ;
3885    for ( uInt i = 0 ; i < nrow - 1 ; i++ ) {
3886      double gap = 2.0 * timeSep[i] / ( interval[i] + interval[i+1] ) ;
3887      //cout << "gap[" << i << "]=" << setw(5) << gap << endl ;
3888      if ( gap > 1.1 ) {
3889        glist.push_back( i ) ;
3890      }
3891    }
3892    Vector<uInt> gaplist( glist ) ;
3893    //cout << "gaplist = " << gaplist << endl ;
3894    uInt newid = 0 ;
3895    for ( uInt i = 0 ; i < nrow ; i++ ) {
3896      scanCol.put( i, newid ) ;
3897      if ( i == gaplist[newid] ) {
3898        newid++ ;
3899      }
3900    }
3901    //cout << "new scancol = " << scanCol.getColumn() << endl ;
3902    vector< CountedPtr<Scantable> > tmp( 1, soff ) ;
3903    CountedPtr<Scantable> aoff = average( tmp, masks, "TINT", "SCAN" ) ;
3904    //cout << "aoff.nrow = " << aoff->nrow() << endl ;
3905    s->unsetSelection() ;
3906    sel.reset() ;
3907    types.clear() ;
3908   
3909    // on scan
3910    bool insitu = insitu_ ;
3911    insitu_ = false ;
3912    CountedPtr<Scantable> out = getScantable( s, true ) ;
3913    insitu_ = insitu ;
3914    types.push_back( SrcType::PSON ) ;
3915    sel.setTypes( types ) ;
3916    s->setSelection( sel ) ;
3917    TableCopy::copyRows( out->table(), s->table() ) ;
3918    s->unsetSelection() ;
3919    sel.reset() ;
3920    types.clear() ;
3921   
3922    // process each on scan
3923    ArrayColumn<Float> tsysCol ;
3924    tsysCol.attach( out->table(), "TSYS" ) ;
3925    for ( int i = 0 ; i < out->nrow() ; i++ ) {
3926      vector<float> sp = getCalibratedSpectra( out, aoff, i ) ;
3927      out->setSpectrum( sp, i ) ;
3928    }
3929
3930    // flux unit
3931    out->setFluxUnit( "K" ) ;
3932
3933    return out ;
3934  }
3935}
3936
3937CountedPtr<Scantable> STMath::cwcalfs( const CountedPtr<Scantable>& s,
3938                                       const String antname )
3939{
3940  vector<int> types ;
3941
3942  // APEX calibration mode
3943  int apexcalmode = 1 ;
3944 
3945  if ( antname.find( "APEX" ) != string::npos ) {
3946    // check if off scan exists or not
3947    STSelector sel = STSelector() ;
3948    //sel.setName( offstr1 ) ;
3949    types.push_back( SrcType::FLOOFF ) ;
3950    sel.setTypes( types ) ;
3951    try {
3952      s->setSelection( sel ) ;
3953    }
3954    catch ( AipsError &e ) {
3955      apexcalmode = 0 ;
3956    }
3957    sel.reset() ;
3958  }
3959  s->unsetSelection() ;
3960  types.clear() ;
3961
3962  vector<bool> masks = s->getMask( 0 ) ;
3963  CountedPtr<Scantable> ssig, sref ;
3964  CountedPtr<Scantable> out ;
3965
3966  if ( antname.find( "APEX" ) != string::npos ) {
3967    // APEX calibration
3968    // sky scan
3969    STSelector sel = STSelector() ;
3970    types.push_back( SrcType::FLOSKY ) ;
3971    sel.setTypes( types ) ;
3972    s->setSelection( sel ) ;
3973    vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ;
3974    CountedPtr<Scantable> askylo = average( tmp, masks, "TINT", "SCAN" ) ;
3975    s->unsetSelection() ;
3976    sel.reset() ;
3977    types.clear() ;
3978    types.push_back( SrcType::FHISKY ) ;
3979    sel.setTypes( types ) ;
3980    s->setSelection( sel ) ;
3981    tmp.clear() ;
3982    tmp.push_back( getScantable( s, false ) ) ;
3983    CountedPtr<Scantable> askyhi = average( tmp, masks, "TINT", "SCAN" ) ;
3984    s->unsetSelection() ;
3985    sel.reset() ;
3986    types.clear() ;
3987   
3988    // hot scan
3989    types.push_back( SrcType::FLOHOT ) ;
3990    sel.setTypes( types ) ;
3991    s->setSelection( sel ) ;
3992    tmp.clear() ;
3993    tmp.push_back( getScantable( s, false ) ) ;
3994    CountedPtr<Scantable> ahotlo = average( tmp, masks, "TINT", "SCAN" ) ;
3995    s->unsetSelection() ;
3996    sel.reset() ;
3997    types.clear() ;
3998    types.push_back( SrcType::FHIHOT ) ;
3999    sel.setTypes( types ) ;
4000    s->setSelection( sel ) ;
4001    tmp.clear() ;
4002    tmp.push_back( getScantable( s, false ) ) ;
4003    CountedPtr<Scantable> ahothi = average( tmp, masks, "TINT", "SCAN" ) ;
4004    s->unsetSelection() ;
4005    sel.reset() ;
4006    types.clear() ;
4007   
4008    // cold scan
4009    CountedPtr<Scantable> acoldlo, acoldhi ;
4010//     types.push_back( SrcType::FLOCOLD ) ;
4011//     sel.setTypes( types ) ;
4012//     s->setSelection( sel ) ;
4013//     tmp.clear() ;
4014//     tmp.push_back( getScantable( s, false ) ) ;
4015//     CountedPtr<Scantable> acoldlo = average( tmp, masks, "TINT", "SCAN" ) ;
4016//     s->unsetSelection() ;
4017//     sel.reset() ;
4018//     types.clear() ;
4019//     types.push_back( SrcType::FHICOLD ) ;
4020//     sel.setTypes( types ) ;
4021//     s->setSelection( sel ) ;
4022//     tmp.clear() ;
4023//     tmp.push_back( getScantable( s, false ) ) ;
4024//     CountedPtr<Scantable> acoldhi = average( tmp, masks, "TINT", "SCAN" ) ;
4025//     s->unsetSelection() ;
4026//     sel.reset() ;
4027//     types.clear() ;
4028
4029    // ref scan
4030    bool insitu = insitu_ ;
4031    insitu_ = false ;
4032    sref = getScantable( s, true ) ;
4033    insitu_ = insitu ;
4034    types.push_back( SrcType::FSLO ) ;
4035    sel.setTypes( types ) ;
4036    s->setSelection( sel ) ;
4037    TableCopy::copyRows( sref->table(), s->table() ) ;
4038    s->unsetSelection() ;
4039    sel.reset() ;
4040    types.clear() ;
4041   
4042    // sig scan
4043    insitu_ = false ;
4044    ssig = getScantable( s, true ) ;
4045    insitu_ = insitu ;
4046    types.push_back( SrcType::FSHI ) ;
4047    sel.setTypes( types ) ;
4048    s->setSelection( sel ) ;
4049    TableCopy::copyRows( ssig->table(), s->table() ) ;
4050    s->unsetSelection() ;
4051    sel.reset() ; 
4052    types.clear() ;
4053         
4054    if ( apexcalmode == 0 ) {
4055      // APEX fs data without off scan
4056      // process each sig and ref scan
4057      ArrayColumn<Float> tsysCollo ;
4058      tsysCollo.attach( ssig->table(), "TSYS" ) ;
4059      ArrayColumn<Float> tsysColhi ;
4060      tsysColhi.attach( sref->table(), "TSYS" ) ;
4061      for ( int i = 0 ; i < ssig->nrow() ; i++ ) {
4062        vector< CountedPtr<Scantable> > sky( 2 ) ;
4063        sky[0] = askylo ;
4064        sky[1] = askyhi ;
4065        vector< CountedPtr<Scantable> > hot( 2 ) ;
4066        hot[0] = ahotlo ;
4067        hot[1] = ahothi ;
4068        vector< CountedPtr<Scantable> > cold( 2 ) ;
4069        //cold[0] = acoldlo ;
4070        //cold[1] = acoldhi ;
4071        vector<float> sp = getFSCalibratedSpectra( ssig, sref, sky, hot, cold, i ) ;
4072        ssig->setSpectrum( sp, i ) ;
4073        string reftime = ssig->getTime( i ) ;
4074        vector<int> ii( 1, ssig->getIF( i ) ) ;
4075        vector<int> ib( 1, ssig->getBeam( i ) ) ;
4076        vector<int> ip( 1, ssig->getPol( i ) ) ;
4077        sel.setIFs( ii ) ;
4078        sel.setBeams( ib ) ;
4079        sel.setPolarizations( ip ) ;
4080        askylo->setSelection( sel ) ;
4081        vector<float> sptsys = getTsysFromTime( reftime, askylo, "linear" ) ;
4082        const Vector<Float> Vtsyslo( sptsys ) ;
4083        tsysCollo.put( i, Vtsyslo ) ;
4084        askylo->unsetSelection() ;
4085        sel.reset() ;
4086        sky[0] = askyhi ;
4087        sky[1] = askylo ;
4088        hot[0] = ahothi ;
4089        hot[1] = ahotlo ;
4090        cold[0] = acoldhi ;
4091        cold[1] = acoldlo ;
4092        sp = getFSCalibratedSpectra( sref, ssig, sky, hot, cold, i ) ;
4093        sref->setSpectrum( sp, i ) ;
4094        reftime = sref->getTime( i ) ;
4095        ii[0] = sref->getIF( i )  ;
4096        ib[0] = sref->getBeam( i ) ;
4097        ip[0] = sref->getPol( i ) ;
4098        sel.setIFs( ii ) ;
4099        sel.setBeams( ib ) ;
4100        sel.setPolarizations( ip ) ;
4101        askyhi->setSelection( sel ) ;   
4102        sptsys = getTsysFromTime( reftime, askyhi, "linear" ) ;
4103        const Vector<Float> Vtsyshi( sptsys ) ;
4104        tsysColhi.put( i, Vtsyshi ) ;
4105        askyhi->unsetSelection() ;
4106        sel.reset() ;
4107      }
4108    }
4109    else if ( apexcalmode == 1 ) {
4110      // APEX fs data with off scan
4111      // off scan
4112      types.push_back( SrcType::FLOOFF ) ;
4113      sel.setTypes( types ) ;
4114      s->setSelection( sel ) ;
4115      tmp.clear() ;
4116      tmp.push_back( getScantable( s, false ) ) ;
4117      CountedPtr<Scantable> aofflo = average( tmp, masks, "TINT", "SCAN" ) ;
4118      s->unsetSelection() ;
4119      sel.reset() ;
4120      types.clear() ;
4121      types.push_back( SrcType::FHIOFF ) ;
4122      sel.setTypes( types ) ;
4123      s->setSelection( sel ) ;
4124      tmp.clear() ;
4125      tmp.push_back( getScantable( s, false ) ) ;
4126      CountedPtr<Scantable> aoffhi = average( tmp, masks, "TINT", "SCAN" ) ;
4127      s->unsetSelection() ;
4128      sel.reset() ;
4129      types.clear() ;
4130     
4131      // process each sig and ref scan
4132      ArrayColumn<Float> tsysCollo ;
4133      tsysCollo.attach( ssig->table(), "TSYS" ) ;
4134      ArrayColumn<Float> tsysColhi ;
4135      tsysColhi.attach( sref->table(), "TSYS" ) ;
4136      for ( int i = 0 ; i < ssig->nrow() ; i++ ) {
4137        vector<float> sp = getCalibratedSpectra( ssig, aofflo, askylo, ahotlo, acoldlo, i, antname ) ;
4138        ssig->setSpectrum( sp, i ) ;
4139        sp = getCalibratedSpectra( sref, aoffhi, askyhi, ahothi, acoldhi, i, antname ) ;
4140        string reftime = ssig->getTime( i ) ;
4141        vector<int> ii( 1, ssig->getIF( i ) ) ;
4142        vector<int> ib( 1, ssig->getBeam( i ) ) ;
4143        vector<int> ip( 1, ssig->getPol( i ) ) ;
4144        sel.setIFs( ii ) ;
4145        sel.setBeams( ib ) ;
4146        sel.setPolarizations( ip ) ;
4147        askylo->setSelection( sel ) ;
4148        vector<float> sptsys = getTsysFromTime( reftime, askylo, "linear" ) ;
4149        const Vector<Float> Vtsyslo( sptsys ) ;
4150        tsysCollo.put( i, Vtsyslo ) ;
4151        askylo->unsetSelection() ;
4152        sel.reset() ;
4153        sref->setSpectrum( sp, i ) ;
4154        reftime = sref->getTime( i ) ;
4155        ii[0] = sref->getIF( i )  ;
4156        ib[0] = sref->getBeam( i ) ;
4157        ip[0] = sref->getPol( i ) ;
4158        sel.setIFs( ii ) ;
4159        sel.setBeams( ib ) ;
4160        sel.setPolarizations( ip ) ;
4161        askyhi->setSelection( sel ) ;   
4162        sptsys = getTsysFromTime( reftime, askyhi, "linear" ) ;
4163        const Vector<Float> Vtsyshi( sptsys ) ;
4164        tsysColhi.put( i, Vtsyshi ) ;
4165        askyhi->unsetSelection() ;
4166        sel.reset() ;
4167      }
4168    }
4169  }
4170  else {
4171    // non-APEX fs data
4172    // sky scan
4173    STSelector sel = STSelector() ;
4174    types.push_back( SrcType::SKY ) ;
4175    sel.setTypes( types ) ;
4176    s->setSelection( sel ) ;
4177    vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ;
4178    CountedPtr<Scantable> asky = average( tmp, masks, "TINT", "SCAN" ) ;
4179    s->unsetSelection() ;
4180    sel.reset() ;
4181    types.clear() ;
4182   
4183    // hot scan
4184    types.push_back( SrcType::HOT ) ;
4185    sel.setTypes( types ) ;
4186    s->setSelection( sel ) ;
4187    tmp.clear() ;
4188    tmp.push_back( getScantable( s, false ) ) ;
4189    CountedPtr<Scantable> ahot = average( tmp, masks, "TINT", "SCAN" ) ;
4190    s->unsetSelection() ;
4191    sel.reset() ;
4192    types.clear() ;
4193
4194    // cold scan
4195    CountedPtr<Scantable> acold ;
4196//     types.push_back( SrcType::COLD ) ;
4197//     sel.setTypes( types ) ;
4198//     s->setSelection( sel ) ;
4199//     tmp.clear() ;
4200//     tmp.push_back( getScantable( s, false ) ) ;
4201//     CountedPtr<Scantable> acold = average( tmp, masks, "TINT", "SCAN" ) ;
4202//     s->unsetSelection() ;
4203//     sel.reset() ;
4204//     types.clear() ;
4205   
4206    // ref scan
4207    bool insitu = insitu_ ;
4208    insitu_ = false ;
4209    sref = getScantable( s, true ) ;
4210    insitu_ = insitu ;
4211    types.push_back( SrcType::FSOFF ) ;
4212    sel.setTypes( types ) ;
4213    s->setSelection( sel ) ;
4214    TableCopy::copyRows( sref->table(), s->table() ) ;
4215    s->unsetSelection() ;
4216    sel.reset() ;
4217    types.clear() ;
4218   
4219    // sig scan
4220    insitu_ = false ;
4221    ssig = getScantable( s, true ) ;
4222    insitu_ = insitu ;
4223    types.push_back( SrcType::FSON ) ;
4224    sel.setTypes( types ) ;
4225    s->setSelection( sel ) ;
4226    TableCopy::copyRows( ssig->table(), s->table() ) ;
4227    s->unsetSelection() ;
4228    sel.reset() ;
4229    types.clear() ;
4230
4231    // process each sig and ref scan
4232    ArrayColumn<Float> tsysColsig ;
4233    tsysColsig.attach( ssig->table(), "TSYS" ) ;
4234    ArrayColumn<Float> tsysColref ;
4235    tsysColref.attach( ssig->table(), "TSYS" ) ;
4236    for ( int i = 0 ; i < ssig->nrow() ; i++ ) {
4237      vector<float> sp = getFSCalibratedSpectra( ssig, sref, asky, ahot, acold, i ) ;
4238      ssig->setSpectrum( sp, i ) ;
4239      string reftime = ssig->getTime( i ) ;
4240      vector<int> ii( 1, ssig->getIF( i ) ) ;
4241      vector<int> ib( 1, ssig->getBeam( i ) ) ;
4242      vector<int> ip( 1, ssig->getPol( i ) ) ;
4243      sel.setIFs( ii ) ;
4244      sel.setBeams( ib ) ;
4245      sel.setPolarizations( ip ) ;
4246      asky->setSelection( sel ) ;
4247      vector<float> sptsys = getTsysFromTime( reftime, asky, "linear" ) ;
4248      const Vector<Float> Vtsys( sptsys ) ;
4249      tsysColsig.put( i, Vtsys ) ;
4250      asky->unsetSelection() ;
4251      sel.reset() ;
4252      sp = getFSCalibratedSpectra( sref, ssig, asky, ahot, acold, i ) ;
4253      sref->setSpectrum( sp, i ) ;
4254      tsysColref.put( i, Vtsys ) ;
4255    }
4256  }
4257
4258  // do folding if necessary
4259  Table sigtab = ssig->table() ;
4260  Table reftab = sref->table() ;
4261  ScalarColumn<uInt> sigifnoCol ;
4262  ScalarColumn<uInt> refifnoCol ;
4263  ScalarColumn<uInt> sigfidCol ;
4264  ScalarColumn<uInt> reffidCol ;
4265  Int nchan = (Int)ssig->nchan() ;
4266  sigifnoCol.attach( sigtab, "IFNO" ) ;
4267  refifnoCol.attach( reftab, "IFNO" ) ;
4268  sigfidCol.attach( sigtab, "FREQ_ID" ) ;
4269  reffidCol.attach( reftab, "FREQ_ID" ) ;
4270  Vector<uInt> sfids( sigfidCol.getColumn() ) ;
4271  Vector<uInt> rfids( reffidCol.getColumn() ) ;
4272  vector<uInt> sfids_unique ;
4273  vector<uInt> rfids_unique ;
4274  vector<uInt> sifno_unique ;
4275  vector<uInt> rifno_unique ;
4276  for ( uInt i = 0 ; i < sfids.nelements() ; i++ ) {
4277    if ( count( sfids_unique.begin(), sfids_unique.end(), sfids[i] ) == 0 ) {
4278      sfids_unique.push_back( sfids[i] ) ;
4279      sifno_unique.push_back( ssig->getIF( i ) ) ;
4280    }
4281    if ( count( rfids_unique.begin(), rfids_unique.end(),  rfids[i] ) == 0 ) {
4282      rfids_unique.push_back( rfids[i] ) ;
4283      rifno_unique.push_back( sref->getIF( i ) ) ;
4284    }
4285  }
4286  double refpix_sig, refval_sig, increment_sig ;
4287  double refpix_ref, refval_ref, increment_ref ;
4288  vector< CountedPtr<Scantable> > tmp( sfids_unique.size() ) ;
4289  for ( uInt i = 0 ; i < sfids_unique.size() ; i++ ) {
4290    ssig->frequencies().getEntry( refpix_sig, refval_sig, increment_sig, sfids_unique[i] ) ;
4291    sref->frequencies().getEntry( refpix_ref, refval_ref, increment_ref, rfids_unique[i] ) ;
4292    if ( refpix_sig == refpix_ref ) {
4293      double foffset = refval_ref - refval_sig ;
4294      int choffset = static_cast<int>(foffset/increment_sig) ;
4295      double doffset = foffset / increment_sig ;
4296      if ( abs(choffset) >= nchan ) {
4297        LogIO os( LogOrigin( "STMath", "cwcalfs", WHERE ) ) ;
4298        os << "FREQ_ID=[" << sfids_unique[i] << "," << rfids_unique[i] << "]: out-band frequency switching, no folding" << LogIO::POST ;
4299        os << "Just return signal data" << LogIO::POST ;
4300        //std::vector< CountedPtr<Scantable> > tabs ;
4301        //tabs.push_back( ssig ) ;
4302        //tabs.push_back( sref ) ;
4303        //out = merge( tabs ) ;
4304        tmp[i] = ssig ;
4305      }
4306      else {
4307        STSelector sel = STSelector() ;
4308        vector<int> v( 1, sifno_unique[i] ) ;
4309        sel.setIFs( v ) ;
4310        ssig->setSelection( sel ) ;
4311        sel.reset() ;
4312        v[0] = rifno_unique[i] ;
4313        sel.setIFs( v ) ;
4314        sref->setSelection( sel ) ;
4315        sel.reset() ;
4316        if ( antname.find( "APEX" ) != string::npos ) {
4317          tmp[i] = dofold( ssig, sref, 0.5*doffset, -0.5*doffset ) ;
4318          //tmp[i] = dofold( ssig, sref, doffset ) ;
4319        }
4320        else {
4321          tmp[i] = dofold( ssig, sref, doffset ) ;
4322        }
4323        ssig->unsetSelection() ;
4324        sref->unsetSelection() ;
4325      }
4326    }
4327  }
4328
4329  if ( tmp.size() > 1 ) {
4330    out = merge( tmp ) ;
4331  }
4332  else {
4333    out = tmp[0] ;
4334  }
4335
4336  // flux unit
4337  out->setFluxUnit( "K" ) ;
4338
4339  return out ;
4340}
4341
4342CountedPtr<Scantable> STMath::almacalfs( const CountedPtr<Scantable>& s )
4343{
4344  (void) s; //currently unused
4345  CountedPtr<Scantable> out ;
4346
4347  return out ;
4348}
4349
4350vector<float> STMath::getSpectrumFromTime( string reftime,
4351                                           CountedPtr<Scantable>& s,
4352                                           string mode )
4353{
4354  LogIO os( LogOrigin( "STMath", "getSpectrumFromTime", WHERE ) ) ;
4355  vector<float> sp ;
4356
4357  if ( s->nrow() == 0 ) {
4358    os << LogIO::SEVERE << "No spectra in the input scantable. Return empty spectrum." << LogIO::POST ;
4359    return sp ;
4360  }
4361  else if ( s->nrow() == 1 ) {
4362    //os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ;
4363    return s->getSpectrum( 0 ) ;
4364  }
4365  else {
4366    vector<int> idx = getRowIdFromTime( reftime, s ) ;
4367    if ( mode == "before" ) {
4368      int id = -1 ;
4369      if ( idx[0] != -1 ) {
4370        id = idx[0] ;
4371      }
4372      else if ( idx[1] != -1 ) {
4373        os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
4374        id = idx[1] ;
4375      }
4376      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4377      sp = s->getSpectrum( id ) ;
4378    }
4379    else if ( mode == "after" ) {
4380      int id = -1 ;
4381      if ( idx[1] != -1 ) {
4382        id = idx[1] ;
4383      }
4384      else if ( idx[0] != -1 ) {
4385        os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
4386        id = idx[1] ;
4387      }
4388      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4389      sp = s->getSpectrum( id ) ;
4390    }
4391    else if ( mode == "nearest" ) {
4392      int id = -1 ;
4393      if ( idx[0] == -1 ) {
4394        id = idx[1] ;
4395      }
4396      else if ( idx[1] == -1 ) {
4397        id = idx[0] ;
4398      }
4399      else if ( idx[0] == idx[1] ) {
4400        id = idx[0] ;
4401      }
4402      else {
4403        //double t0 = getMJD( s->getTime( idx[0] ) ) ;
4404        //double t1 = getMJD( s->getTime( idx[1] ) ) ;
4405        double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
4406        double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
4407        double tref = getMJD( reftime ) ;
4408        if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
4409          id = idx[1] ;
4410        }
4411        else {
4412          id = idx[0] ;
4413        }
4414      }
4415      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4416      sp = s->getSpectrum( id ) ;     
4417    }
4418    else if ( mode == "linear" ) {
4419      if ( idx[0] == -1 ) {
4420        // use after
4421        os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
4422        int id = idx[1] ;
4423        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4424        sp = s->getSpectrum( id ) ;
4425      }
4426      else if ( idx[1] == -1 ) {
4427        // use before
4428        os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
4429        int id = idx[0] ;
4430        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4431        sp = s->getSpectrum( id ) ;
4432      }
4433      else if ( idx[0] == idx[1] ) {
4434        // use before
4435        //os << "No need to interporate." << LogIO::POST ;
4436        int id = idx[0] ;
4437        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4438        sp = s->getSpectrum( id ) ;
4439      }
4440      else {
4441        // do interpolation
4442        //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
4443        //double t0 = getMJD( s->getTime( idx[0] ) ) ;
4444        //double t1 = getMJD( s->getTime( idx[1] ) ) ;
4445        double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
4446        double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
4447        double tref = getMJD( reftime ) ;
4448        vector<float> sp0 = s->getSpectrum( idx[0] ) ;
4449        vector<float> sp1 = s->getSpectrum( idx[1] ) ;
4450        for ( unsigned int i = 0 ; i < sp0.size() ; i++ ) {
4451          float v = ( sp1[i] - sp0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + sp0[i] ;
4452          sp.push_back( v ) ;
4453        }
4454      }
4455    }
4456    else {
4457      os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
4458    }
4459    return sp ;
4460  }
4461}
4462
4463vector<float> STMath::getSpectrumFromTime( double reftime,
4464                                           CountedPtr<Scantable>& s,
4465                                           string mode )
4466{
4467  LogIO os( LogOrigin( "STMath", "getSpectrumFromTime", WHERE ) ) ;
4468  vector<float> sp ;
4469
4470  if ( s->nrow() == 0 ) {
4471    os << LogIO::SEVERE << "No spectra in the input scantable. Return empty spectrum." << LogIO::POST ;
4472    return sp ;
4473  }
4474  else if ( s->nrow() == 1 ) {
4475    //os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ;
4476    return s->getSpectrum( 0 ) ;
4477  }
4478  else {
4479    ROScalarColumn<Double> timeCol( s->table(), "TIME" ) ;
4480    Vector<Double> timeVec = timeCol.getColumn() ;
4481    vector<int> idx = getRowIdFromTime( reftime, timeVec ) ;
4482    if ( mode == "before" ) {
4483      int id = -1 ;
4484      if ( idx[0] != -1 ) {
4485        id = idx[0] ;
4486      }
4487      else if ( idx[1] != -1 ) {
4488        os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
4489        id = idx[1] ;
4490      }
4491      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4492      sp = s->getSpectrum( id ) ;
4493    }
4494    else if ( mode == "after" ) {
4495      int id = -1 ;
4496      if ( idx[1] != -1 ) {
4497        id = idx[1] ;
4498      }
4499      else if ( idx[0] != -1 ) {
4500        os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
4501        id = idx[1] ;
4502      }
4503      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4504      sp = s->getSpectrum( id ) ;
4505    }
4506    else if ( mode == "nearest" ) {
4507      int id = -1 ;
4508      if ( idx[0] == -1 ) {
4509        id = idx[1] ;
4510      }
4511      else if ( idx[1] == -1 ) {
4512        id = idx[0] ;
4513      }
4514      else if ( idx[0] == idx[1] ) {
4515        id = idx[0] ;
4516      }
4517      else {
4518        //double t0 = getMJD( s->getTime( idx[0] ) ) ;
4519        //double t1 = getMJD( s->getTime( idx[1] ) ) ;
4520//         double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
4521//         double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
4522        double t0 = timeVec[idx[0]] ;
4523        double t1 = timeVec[idx[1]] ;
4524//         cout << "t0-t0c=" << t0-t0c << endl ;
4525//         cout << "t1-t1c=" << t1-t1c << endl ;
4526//         double tref = getMJD( reftime ) ;
4527        double tref = reftime ;
4528        if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
4529          id = idx[1] ;
4530        }
4531        else {
4532          id = idx[0] ;
4533        }
4534      }
4535      //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4536      sp = s->getSpectrum( id ) ;     
4537    }
4538    else if ( mode == "linear" ) {
4539      if ( idx[0] == -1 ) {
4540        // use after
4541        os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
4542        int id = idx[1] ;
4543        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4544        sp = s->getSpectrum( id ) ;
4545      }
4546      else if ( idx[1] == -1 ) {
4547        // use before
4548        os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
4549        int id = idx[0] ;
4550        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4551        sp = s->getSpectrum( id ) ;
4552      }
4553      else if ( idx[0] == idx[1] ) {
4554        // use before
4555        //os << "No need to interporate." << LogIO::POST ;
4556        int id = idx[0] ;
4557        //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4558        sp = s->getSpectrum( id ) ;
4559      }
4560      else {
4561        // do interpolation
4562        //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
4563        //double t0 = getMJD( s->getTime( idx[0] ) ) ;
4564        //double t1 = getMJD( s->getTime( idx[1] ) ) ;
4565//         double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
4566//         double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
4567        double t0 = timeVec[idx[0]] ;
4568        double t1 = timeVec[idx[1]] ;
4569//         cout << "t0-t0c=" << t0-t0c << endl ;
4570//         cout << "t1-t1c=" << t1-t1c << endl ;
4571//         double tref = getMJD( reftime ) ;
4572        double tref = reftime ;
4573        vector<float> sp0 = s->getSpectrum( idx[0] ) ;
4574        vector<float> sp1 = s->getSpectrum( idx[1] ) ;
4575        for ( unsigned int i = 0 ; i < sp0.size() ; i++ ) {
4576          float v = ( sp1[i] - sp0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + sp0[i] ;
4577          sp.push_back( v ) ;
4578        }
4579      }
4580    }
4581    else {
4582      os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
4583    }
4584    return sp ;
4585  }
4586}
4587
4588double STMath::getMJD( string strtime )
4589{
4590  if ( strtime.find("/") == string::npos ) {
4591    // MJD time string
4592    return atof( strtime.c_str() ) ;
4593  }
4594  else {
4595    // string in YYYY/MM/DD/HH:MM:SS format
4596    uInt year = atoi( strtime.substr( 0, 4 ).c_str() ) ;
4597    uInt month = atoi( strtime.substr( 5, 2 ).c_str() ) ;
4598    uInt day = atoi( strtime.substr( 8, 2 ).c_str() ) ;
4599    uInt hour = atoi( strtime.substr( 11, 2 ).c_str() ) ;
4600    uInt minute = atoi( strtime.substr( 14, 2 ).c_str() ) ;
4601    uInt sec = atoi( strtime.substr( 17, 2 ).c_str() ) ;
4602    Time t( year, month, day, hour, minute, sec ) ;
4603    return t.modifiedJulianDay() ;
4604  }
4605}
4606
4607vector<int> STMath::getRowIdFromTime( string reftime, CountedPtr<Scantable> &s )
4608{
4609  double reft = getMJD( reftime ) ;
4610  double dtmin = 1.0e100 ;
4611  double dtmax = -1.0e100 ;
4612  vector<double> dt ;
4613  int just_before = -1 ;
4614  int just_after = -1 ;
4615  for ( int i = 0 ; i < s->nrow() ; i++ ) {
4616    dt.push_back( getMJD( s->getTime( i ) ) - reft ) ;
4617  }
4618  for ( unsigned int i = 0 ; i < dt.size() ; i++ ) {
4619    if ( dt[i] > 0.0 ) {
4620      // after reftime
4621      if ( dt[i] < dtmin ) {
4622        just_after = i ;
4623        dtmin = dt[i] ;
4624      }
4625    }
4626    else if ( dt[i] < 0.0 ) {
4627      // before reftime
4628      if ( dt[i] > dtmax ) {
4629        just_before = i ;
4630        dtmax = dt[i] ;
4631      }
4632    }
4633    else {
4634      // just a reftime
4635      just_before = i ;
4636      just_after = i ;
4637      dtmax = 0 ;
4638      dtmin = 0 ;
4639      break ;
4640    }
4641  }
4642
4643  vector<int> v ;
4644  v.push_back( just_before ) ;
4645  v.push_back( just_after ) ;
4646
4647  return v ;
4648}
4649
4650vector<int> STMath::getRowIdFromTime( double reftime, Vector<Double> &t )
4651{
4652//   double reft = reftime ;
4653  double dtmin = 1.0e100 ;
4654  double dtmax = -1.0e100 ;
4655//   vector<double> dt ;
4656  int just_before = -1 ;
4657  int just_after = -1 ;
4658  //cout << setprecision(24) << reft << endl ;
4659//   ROScalarColumn<Double> timeCol( s->table(), "TIME" ) ;
4660//   for ( int i = 0 ; i < s->nrow() ; i++ ) {
4661//     cout << setprecision(24) << timeCol(i) << endl ;
4662//     //dt.push_back( getMJD( s->getTime( i ) ) - reft ) ;
4663//     dt.push_back( timeCol(i) - reft ) ;
4664//   }
4665  Vector<Double> dt = t - reftime ;
4666  for ( unsigned int i = 0 ; i < dt.size() ; i++ ) {
4667    if ( dt[i] > 0.0 ) {
4668      // after reftime
4669      if ( dt[i] < dtmin ) {
4670        just_after = i ;
4671        dtmin = dt[i] ;
4672      }
4673    }
4674    else if ( dt[i] < 0.0 ) {
4675      // before reftime
4676      if ( dt[i] > dtmax ) {
4677        just_before = i ;
4678        dtmax = dt[i] ;
4679      }
4680    }
4681    else {
4682      // just a reftime
4683      just_before = i ;
4684      just_after = i ;
4685      dtmax = 0 ;
4686      dtmin = 0 ;
4687      break ;
4688    }
4689  }
4690
4691  vector<int> v(2) ;
4692  v[0] = just_before ;
4693  v[1] = just_after ;
4694
4695  return v ;
4696}
4697
4698vector<float> STMath::getTcalFromTime( string reftime,
4699                                       CountedPtr<Scantable>& s,
4700                                       string mode )
4701{
4702  LogIO os( LogOrigin( "STMath", "getTcalFromTime", WHERE ) ) ;
4703  vector<float> tcal ;
4704  STTcal tcalTable = s->tcal() ;
4705  String time ;
4706  Vector<Float> tcalval ;
4707  if ( s->nrow() == 0 ) {
4708    os << LogIO::SEVERE << "No row in the input scantable. Return empty tcal." << LogIO::POST ;
4709    return tcal ;
4710  }
4711  else if ( s->nrow() == 1 ) {
4712    uInt tcalid = s->getTcalId( 0 ) ;
4713    //os << "use row " << 0 << " (tcalid = " << tcalid << ")" << LogIO::POST ;
4714    tcalTable.getEntry( time, tcalval, tcalid ) ;
4715    tcalval.tovector( tcal ) ;
4716    return tcal ;
4717  }
4718  else {
4719    vector<int> idx = getRowIdFromTime( reftime, s ) ;
4720    if ( mode == "before" ) {
4721      int id = -1 ;
4722      if ( idx[0] != -1 ) {
4723        id = idx[0] ;
4724      }
4725      else if ( idx[1] != -1 ) {
4726        os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
4727        id = idx[1] ;
4728      }
4729      uInt tcalid = s->getTcalId( id ) ;
4730      //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
4731      tcalTable.getEntry( time, tcalval, tcalid ) ;
4732      tcalval.tovector( tcal ) ;
4733    }
4734    else if ( mode == "after" ) {
4735      int id = -1 ;
4736      if ( idx[1] != -1 ) {
4737        id = idx[1] ;
4738      }
4739      else if ( idx[0] != -1 ) {
4740        os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
4741        id = idx[1] ;
4742      }
4743      uInt tcalid = s->getTcalId( id ) ;
4744      //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
4745      tcalTable.getEntry( time, tcalval, tcalid ) ;
4746      tcalval.tovector( tcal ) ;
4747    }
4748    else if ( mode == "nearest" ) {
4749      int id = -1 ;
4750      if ( idx[0] == -1 ) {
4751        id = idx[1] ;
4752      }
4753      else if ( idx[1] == -1 ) {
4754        id = idx[0] ;
4755      }
4756      else if ( idx[0] == idx[1] ) {
4757        id = idx[0] ;
4758      }
4759      else {
4760        //double t0 = getMJD( s->getTime( idx[0] ) ) ;
4761        //double t1 = getMJD( s->getTime( idx[1] ) ) ;
4762        double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
4763        double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
4764        double tref = getMJD( reftime ) ;
4765        if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
4766          id = idx[1] ;
4767        }
4768        else {
4769          id = idx[0] ;
4770        }
4771      }
4772      uInt tcalid = s->getTcalId( id ) ;
4773      //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
4774      tcalTable.getEntry( time, tcalval, tcalid ) ;
4775      tcalval.tovector( tcal ) ;
4776    }
4777    else if ( mode == "linear" ) {
4778      if ( idx[0] == -1 ) {
4779        // use after
4780        os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
4781        int id = idx[1] ;
4782        uInt tcalid = s->getTcalId( id ) ;
4783        //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
4784        tcalTable.getEntry( time, tcalval, tcalid ) ;
4785        tcalval.tovector( tcal ) ;
4786      }
4787      else if ( idx[1] == -1 ) {
4788        // use before
4789        os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
4790        int id = idx[0] ;
4791        uInt tcalid = s->getTcalId( id ) ;
4792        //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
4793        tcalTable.getEntry( time, tcalval, tcalid ) ;
4794        tcalval.tovector( tcal ) ;
4795      }
4796      else if ( idx[0] == idx[1] ) {
4797        // use before
4798        //os << "No need to interporate." << LogIO::POST ;
4799        int id = idx[0] ;
4800        uInt tcalid = s->getTcalId( id ) ;
4801        //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
4802        tcalTable.getEntry( time, tcalval, tcalid ) ;
4803        tcalval.tovector( tcal ) ;
4804      }
4805      else {
4806        // do interpolation
4807        //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
4808        //double t0 = getMJD( s->getTime( idx[0] ) ) ;
4809        //double t1 = getMJD( s->getTime( idx[1] ) ) ;
4810        double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
4811        double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
4812        double tref = getMJD( reftime ) ;
4813        vector<float> tcal0 ;
4814        vector<float> tcal1 ;
4815        uInt tcalid0 = s->getTcalId( idx[0] ) ;
4816        uInt tcalid1 = s->getTcalId( idx[1] ) ;
4817        tcalTable.getEntry( time, tcalval, tcalid0 ) ;
4818        tcalval.tovector( tcal0 ) ;
4819        tcalTable.getEntry( time, tcalval, tcalid1 ) ;
4820        tcalval.tovector( tcal1 ) ;       
4821        for ( unsigned int i = 0 ; i < tcal0.size() ; i++ ) {
4822          float v = ( tcal1[i] - tcal0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + tcal0[i] ;
4823          tcal.push_back( v ) ;
4824        }
4825      }
4826    }
4827    else {
4828      os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
4829    }
4830    return tcal ;
4831  }
4832}
4833
4834vector<float> STMath::getTsysFromTime( string reftime,
4835                                       CountedPtr<Scantable>& s,
4836                                       string mode )
4837{
4838  LogIO os( LogOrigin( "STMath", "getTsysFromTime", WHERE ) ) ;
4839  ArrayColumn<Float> tsysCol ;
4840  tsysCol.attach( s->table(), "TSYS" ) ;
4841  vector<float> tsys ;
4842  String time ;
4843  Vector<Float> tsysval ;
4844  if ( s->nrow() == 0 ) {
4845    os << LogIO::SEVERE << "No row in the input scantable. Return empty tsys." << LogIO::POST ;
4846    return tsys ;
4847  }
4848  else if ( s->nrow() == 1 ) {
4849    //os << "use row " << 0 << LogIO::POST ;
4850    tsysval = tsysCol( 0 ) ;
4851    tsysval.tovector( tsys ) ;
4852    return tsys ;
4853  }
4854  else {
4855    vector<int> idx = getRowIdFromTime( reftime, s ) ;
4856    if ( mode == "before" ) {
4857      int id = -1 ;
4858      if ( idx[0] != -1 ) {
4859        id = idx[0] ;
4860      }
4861      else if ( idx[1] != -1 ) {
4862        os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
4863        id = idx[1] ;
4864      }
4865      //os << "use row " << id << LogIO::POST ;
4866      tsysval = tsysCol( id ) ;
4867      tsysval.tovector( tsys ) ;
4868    }
4869    else if ( mode == "after" ) {
4870      int id = -1 ;
4871      if ( idx[1] != -1 ) {
4872        id = idx[1] ;
4873      }
4874      else if ( idx[0] != -1 ) {
4875        os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
4876        id = idx[1] ;
4877      }
4878      //os << "use row " << id << LogIO::POST ;
4879      tsysval = tsysCol( id ) ;
4880      tsysval.tovector( tsys ) ;
4881    }
4882    else if ( mode == "nearest" ) {
4883      int id = -1 ;
4884      if ( idx[0] == -1 ) {
4885        id = idx[1] ;
4886      }
4887      else if ( idx[1] == -1 ) {
4888        id = idx[0] ;
4889      }
4890      else if ( idx[0] == idx[1] ) {
4891        id = idx[0] ;
4892      }
4893      else {
4894        //double t0 = getMJD( s->getTime( idx[0] ) ) ;
4895        //double t1 = getMJD( s->getTime( idx[1] ) ) ;
4896        double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
4897        double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
4898        double tref = getMJD( reftime ) ;
4899        if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
4900          id = idx[1] ;
4901        }
4902        else {
4903          id = idx[0] ;
4904        }
4905      }
4906      //os << "use row " << id << LogIO::POST ;
4907      tsysval = tsysCol( id ) ;
4908      tsysval.tovector( tsys ) ;
4909    }
4910    else if ( mode == "linear" ) {
4911      if ( idx[0] == -1 ) {
4912        // use after
4913        os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
4914        int id = idx[1] ;
4915        //os << "use row " << id << LogIO::POST ;
4916        tsysval = tsysCol( id ) ;
4917        tsysval.tovector( tsys ) ;
4918      }
4919      else if ( idx[1] == -1 ) {
4920        // use before
4921        os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
4922        int id = idx[0] ;
4923        //os << "use row " << id << LogIO::POST ;
4924        tsysval = tsysCol( id ) ;
4925        tsysval.tovector( tsys ) ;
4926      }
4927      else if ( idx[0] == idx[1] ) {
4928        // use before
4929        //os << "No need to interporate." << LogIO::POST ;
4930        int id = idx[0] ;
4931        //os << "use row " << id << LogIO::POST ;
4932        tsysval = tsysCol( id ) ;
4933        tsysval.tovector( tsys ) ;
4934      }
4935      else {
4936        // do interpolation
4937        //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
4938        //double t0 = getMJD( s->getTime( idx[0] ) ) ;
4939        //double t1 = getMJD( s->getTime( idx[1] ) ) ;
4940        double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
4941        double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
4942        double tref = getMJD( reftime ) ;
4943        vector<float> tsys0 ;
4944        vector<float> tsys1 ;
4945        tsysval = tsysCol( idx[0] ) ;
4946        tsysval.tovector( tsys0 ) ;
4947        tsysval = tsysCol( idx[1] ) ;
4948        tsysval.tovector( tsys1 ) ;       
4949        for ( unsigned int i = 0 ; i < tsys0.size() ; i++ ) {
4950          float v = ( tsys1[i] - tsys0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + tsys0[i] ;
4951          tsys.push_back( v ) ;
4952        }
4953      }
4954    }
4955    else {
4956      os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
4957    }
4958    return tsys ;
4959  }
4960}
4961
4962vector<float> STMath::getCalibratedSpectra( CountedPtr<Scantable>& on,
4963                                            CountedPtr<Scantable>& off,
4964                                            CountedPtr<Scantable>& sky,
4965                                            CountedPtr<Scantable>& hot,
4966                                            CountedPtr<Scantable>& cold,
4967                                            int index,
4968                                            string antname )
4969{
4970  (void) cold; //currently unused
4971  string reftime = on->getTime( index ) ;
4972  vector<int> ii( 1, on->getIF( index ) ) ;
4973  vector<int> ib( 1, on->getBeam( index ) ) ;
4974  vector<int> ip( 1, on->getPol( index ) ) ;
4975  vector<int> ic( 1, on->getScan( index ) ) ;
4976  STSelector sel = STSelector() ;
4977  sel.setIFs( ii ) ;
4978  sel.setBeams( ib ) ;
4979  sel.setPolarizations( ip ) ;
4980  sky->setSelection( sel ) ;
4981  hot->setSelection( sel ) ;
4982  //cold->setSelection( sel ) ;
4983  off->setSelection( sel ) ;
4984  vector<float> spsky = getSpectrumFromTime( reftime, sky, "linear" ) ;
4985  vector<float> sphot = getSpectrumFromTime( reftime, hot, "linear" ) ;
4986  //vector<float> spcold = getSpectrumFromTime( reftime, cold, "linear" ) ;
4987  vector<float> spoff = getSpectrumFromTime( reftime, off, "linear" ) ;
4988  vector<float> spec = on->getSpectrum( index ) ;
4989  vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ;
4990  vector<float> sp( tcal.size() ) ;
4991  if ( antname.find( "APEX" ) != string::npos ) {
4992    // using gain array
4993    for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
4994      float v = ( ( spec[j] - spoff[j] ) / spoff[j] )
4995        * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
4996      sp[j] = v ;
4997    }
4998  }
4999  else {
5000    // Chopper-Wheel calibration (Ulich & Haas 1976)
5001    for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
5002      float v = ( spec[j] - spoff[j] ) / ( sphot[j] - spsky[j] ) * tcal[j] ;
5003      sp[j] = v ;
5004    }
5005  }
5006  sel.reset() ;
5007  sky->unsetSelection() ;
5008  hot->unsetSelection() ;
5009  //cold->unsetSelection() ;
5010  off->unsetSelection() ;
5011
5012  return sp ;
5013}
5014
5015vector<float> STMath::getCalibratedSpectra( CountedPtr<Scantable>& on,
5016                                            CountedPtr<Scantable>& off,
5017                                            int index )
5018{
5019//   string reftime = on->getTime( index ) ;
5020  ROTableColumn timeCol( on->table(), "TIME" ) ;
5021  double reftime = timeCol.asdouble(index) ;
5022  vector<int> ii( 1, on->getIF( index ) ) ;
5023  vector<int> ib( 1, on->getBeam( index ) ) ;
5024  vector<int> ip( 1, on->getPol( index ) ) ;
5025  vector<int> ic( 1, on->getScan( index ) ) ;
5026  STSelector sel = STSelector() ;
5027  sel.setIFs( ii ) ;
5028  sel.setBeams( ib ) ;
5029  sel.setPolarizations( ip ) ;
5030  off->setSelection( sel ) ;
5031  vector<float> spoff = getSpectrumFromTime( reftime, off, "linear" ) ;
5032  vector<float> spec = on->getSpectrum( index ) ;
5033  //vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ;
5034  //vector<float> tsys = on->getTsysVec( index ) ;
5035  ArrayColumn<Float> tsysCol( on->table(), "TSYS" ) ;
5036  Vector<Float> tsys = tsysCol( index ) ;
5037  vector<float> sp( spec.size() ) ;
5038  // ALMA Calibration
5039  //
5040  // Ta* = Tsys * ( ON - OFF ) / OFF
5041  //
5042  // 2010/01/07 Takeshi Nakazato
5043  unsigned int tsyssize = tsys.nelements() ;
5044  unsigned int spsize = sp.size() ;
5045  for ( unsigned int j = 0 ; j < sp.size() ; j++ ) {
5046    float tscale = 0.0 ;
5047    if ( tsyssize == spsize )
5048      tscale = tsys[j] ;
5049    else
5050      tscale = tsys[0] ;
5051    float v = tscale * ( spec[j] - spoff[j] ) / spoff[j] ;
5052    sp[j] = v ;
5053  }
5054  sel.reset() ;
5055  off->unsetSelection() ;
5056
5057  return sp ;
5058}
5059
5060vector<float> STMath::getFSCalibratedSpectra( CountedPtr<Scantable>& sig,
5061                                              CountedPtr<Scantable>& ref,
5062                                              CountedPtr<Scantable>& sky,
5063                                              CountedPtr<Scantable>& hot,
5064                                              CountedPtr<Scantable>& cold,
5065                                              int index )
5066{
5067  (void) cold; //currently unused
5068  string reftime = sig->getTime( index ) ;
5069  vector<int> ii( 1, sig->getIF( index ) ) ;
5070  vector<int> ib( 1, sig->getBeam( index ) ) ;
5071  vector<int> ip( 1, sig->getPol( index ) ) ;
5072  vector<int> ic( 1, sig->getScan( index ) ) ;
5073  STSelector sel = STSelector() ;
5074  sel.setIFs( ii ) ;
5075  sel.setBeams( ib ) ;
5076  sel.setPolarizations( ip ) ;
5077  sky->setSelection( sel ) ;
5078  hot->setSelection( sel ) ;
5079  //cold->setSelection( sel ) ;
5080  vector<float> spsky = getSpectrumFromTime( reftime, sky, "linear" ) ;
5081  vector<float> sphot = getSpectrumFromTime( reftime, hot, "linear" ) ;
5082  //vector<float> spcold = getSpectrumFromTime( reftime, cold, "linear" ) ;
5083  vector<float> spref = ref->getSpectrum( index ) ;
5084  vector<float> spsig = sig->getSpectrum( index ) ;
5085  vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ;
5086  vector<float> sp( tcal.size() ) ;
5087  for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
5088    float v = tcal[j] * spsky[j] / ( sphot[j] - spsky[j] ) * ( spsig[j] - spref[j] ) / spref[j] ;
5089    sp[j] = v ;
5090  }
5091  sel.reset() ;
5092  sky->unsetSelection() ;
5093  hot->unsetSelection() ;
5094  //cold->unsetSelection() ;
5095
5096  return sp ;
5097}
5098
5099vector<float> STMath::getFSCalibratedSpectra( CountedPtr<Scantable>& sig,
5100                                              CountedPtr<Scantable>& ref,
5101                                              vector< CountedPtr<Scantable> >& sky,
5102                                              vector< CountedPtr<Scantable> >& hot,
5103                                              vector< CountedPtr<Scantable> >& cold,
5104                                              int index )
5105{
5106  (void) cold; //currently unused
5107  string reftime = sig->getTime( index ) ;
5108  vector<int> ii( 1, sig->getIF( index ) ) ;
5109  vector<int> ib( 1, sig->getBeam( index ) ) ;
5110  vector<int> ip( 1, sig->getPol( index ) ) ;
5111  vector<int> ic( 1, sig->getScan( index ) ) ;
5112  STSelector sel = STSelector() ;
5113  sel.setIFs( ii ) ;
5114  sel.setBeams( ib ) ;
5115  sel.setPolarizations( ip ) ;
5116  sky[0]->setSelection( sel ) ;
5117  hot[0]->setSelection( sel ) ;
5118  //cold[0]->setSelection( sel ) ;
5119  vector<float> spskys = getSpectrumFromTime( reftime, sky[0], "linear" ) ;
5120  vector<float> sphots = getSpectrumFromTime( reftime, hot[0], "linear" ) ;
5121  //vector<float> spcolds = getSpectrumFromTime( reftime, cold[0], "linear" ) ;
5122  vector<float> tcals = getTcalFromTime( reftime, sky[0], "linear" ) ;
5123  sel.reset() ;
5124  ii[0] = ref->getIF( index ) ;
5125  sel.setIFs( ii ) ;
5126  sel.setBeams( ib ) ;
5127  sel.setPolarizations( ip ) ;
5128  sky[1]->setSelection( sel ) ;
5129  hot[1]->setSelection( sel ) ;
5130  //cold[1]->setSelection( sel ) ;
5131  vector<float> spskyr = getSpectrumFromTime( reftime, sky[1], "linear" ) ;
5132  vector<float> sphotr = getSpectrumFromTime( reftime, hot[1], "linear" ) ;
5133  //vector<float> spcoldr = getSpectrumFromTime( reftime, cold[1], "linear" ) ;
5134  vector<float> tcalr = getTcalFromTime( reftime, sky[1], "linear" ) ; 
5135  vector<float> spref = ref->getSpectrum( index ) ;
5136  vector<float> spsig = sig->getSpectrum( index ) ;
5137  vector<float> sp( tcals.size() ) ;
5138  for ( unsigned int j = 0 ; j < tcals.size() ; j++ ) {
5139    float v = tcals[j] * spsig[j] / ( sphots[j] - spskys[j] ) - tcalr[j] * spref[j] / ( sphotr[j] - spskyr[j] ) ;
5140    sp[j] = v ;
5141  }
5142  sel.reset() ;
5143  sky[0]->unsetSelection() ;
5144  hot[0]->unsetSelection() ;
5145  //cold[0]->unsetSelection() ;
5146  sky[1]->unsetSelection() ;
5147  hot[1]->unsetSelection() ;
5148  //cold[1]->unsetSelection() ;
5149
5150  return sp ;
5151}
Note: See TracBrowser for help on using the repository browser.