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

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

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