source: trunk/src/STMath.cpp @ 1484

Last change on this file since 1484 was 1484, checked in by Malte Marquarding, 15 years ago

Fix for Ticket #148; forgot to scale TSYS
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 65.6 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 <casa/iomanip.h>
14#include <casa/Exceptions/Error.h>
15#include <casa/Containers/Block.h>
16#include <casa/BasicSL/String.h>
17#include <casa/Arrays/MaskArrLogi.h>
18#include <casa/Arrays/MaskArrMath.h>
19#include <casa/Arrays/ArrayLogical.h>
20#include <casa/Arrays/ArrayMath.h>
21#include <casa/Arrays/Slice.h>
22#include <casa/Arrays/Slicer.h>
23#include <casa/Containers/RecordField.h>
24#include <tables/Tables/TableRow.h>
25#include <tables/Tables/TableVector.h>
26#include <tables/Tables/TabVecMath.h>
27#include <tables/Tables/ExprNode.h>
28#include <tables/Tables/TableRecord.h>
29#include <tables/Tables/TableParse.h>
30#include <tables/Tables/ReadAsciiTable.h>
31#include <tables/Tables/TableIter.h>
32#include <tables/Tables/TableCopy.h>
33#include <scimath/Mathematics/FFTServer.h>
34
35#include <lattices/Lattices/LatticeUtilities.h>
36
37#include <coordinates/Coordinates/SpectralCoordinate.h>
38#include <coordinates/Coordinates/CoordinateSystem.h>
39#include <coordinates/Coordinates/CoordinateUtil.h>
40#include <coordinates/Coordinates/FrequencyAligner.h>
41
42#include <scimath/Mathematics/VectorKernel.h>
43#include <scimath/Mathematics/Convolver.h>
44#include <scimath/Functionals/Polynomial.h>
45
46#include "MathUtils.h"
47#include "RowAccumulator.h"
48#include "STAttr.h"
49#include "STMath.h"
50#include "STSelector.h"
51
52using namespace casa;
53
54using namespace asap;
55
56STMath::STMath(bool insitu) :
57  insitu_(insitu)
58{
59}
60
61
62STMath::~STMath()
63{
64}
65
66CountedPtr<Scantable>
67STMath::average( const std::vector<CountedPtr<Scantable> >& in,
68                 const std::vector<bool>& mask,
69                 const std::string& weight,
70                 const std::string& avmode)
71{
72  if ( avmode == "SCAN" && in.size() != 1 )
73    throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
74                    "Use merge first."));
75  WeightType wtype = stringToWeight(weight);
76
77  // output
78  // clone as this is non insitu
79  bool insitu = insitu_;
80  setInsitu(false);
81  CountedPtr< Scantable > out = getScantable(in[0], true);
82  setInsitu(insitu);
83  std::vector<CountedPtr<Scantable> >::const_iterator stit = in.begin();
84  ++stit;
85  while ( stit != in.end() ) {
86    out->appendToHistoryTable((*stit)->history());
87    ++stit;
88  }
89
90  Table& tout = out->table();
91
92  /// @todo check if all scantables are conformant
93
94  ArrayColumn<Float> specColOut(tout,"SPECTRA");
95  ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
96  ArrayColumn<Float> tsysColOut(tout,"TSYS");
97  ScalarColumn<Double> mjdColOut(tout,"TIME");
98  ScalarColumn<Double> intColOut(tout,"INTERVAL");
99  ScalarColumn<uInt> cycColOut(tout,"CYCLENO");
100  ScalarColumn<uInt> scanColOut(tout,"SCANNO");
101
102  // set up the output table rows. These are based on the structure of the
103  // FIRST scantable in the vector
104  const Table& baset = in[0]->table();
105
106  Block<String> cols(3);
107  cols[0] = String("BEAMNO");
108  cols[1] = String("IFNO");
109  cols[2] = String("POLNO");
110  if ( avmode == "SOURCE" ) {
111    cols.resize(4);
112    cols[3] = String("SRCNAME");
113  }
114  if ( avmode == "SCAN"  && in.size() == 1) {
115    cols.resize(4);
116    cols[3] = String("SCANNO");
117  }
118  uInt outrowCount = 0;
119  TableIterator iter(baset, cols);
120  while (!iter.pastEnd()) {
121    Table subt = iter.table();
122    // copy the first row of this selection into the new table
123    tout.addRow();
124    TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
125    // re-index to 0
126    if ( avmode != "SCAN" && avmode != "SOURCE" ) {
127      scanColOut.put(outrowCount, uInt(0));
128    }
129    ++outrowCount;
130    ++iter;
131  }
132  RowAccumulator acc(wtype);
133  Vector<Bool> cmask(mask);
134  acc.setUserMask(cmask);
135  ROTableRow row(tout);
136  ROArrayColumn<Float> specCol, tsysCol;
137  ROArrayColumn<uChar> flagCol;
138  ROScalarColumn<Double> mjdCol, intCol;
139  ROScalarColumn<Int> scanIDCol;
140
141  Vector<uInt> rowstodelete;
142
143  for (uInt i=0; i < tout.nrow(); ++i) {
144    for ( int j=0; j < int(in.size()); ++j ) {
145      const Table& tin = in[j]->table();
146      const TableRecord& rec = row.get(i);
147      ROScalarColumn<Double> tmp(tin, "TIME");
148      Double td;tmp.get(0,td);
149      Table basesubt = tin(tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
150                       && tin.col("IFNO") == Int(rec.asuInt("IFNO"))
151                       && tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
152      Table subt;
153      if ( avmode == "SOURCE") {
154        subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME") );
155      } else if (avmode == "SCAN") {
156        subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
157      } else {
158        subt = basesubt;
159      }
160      specCol.attach(subt,"SPECTRA");
161      flagCol.attach(subt,"FLAGTRA");
162      tsysCol.attach(subt,"TSYS");
163      intCol.attach(subt,"INTERVAL");
164      mjdCol.attach(subt,"TIME");
165      Vector<Float> spec,tsys;
166      Vector<uChar> flag;
167      Double inter,time;
168      for (uInt k = 0; k < subt.nrow(); ++k ) {
169        flagCol.get(k, flag);
170        Vector<Bool> bflag(flag.shape());
171        convertArray(bflag, flag);
172        /*
173        if ( allEQ(bflag, True) ) {
174        continue;//don't accumulate
175        }
176        */
177        specCol.get(k, spec);
178        tsysCol.get(k, tsys);
179        intCol.get(k, inter);
180        mjdCol.get(k, time);
181        // spectrum has to be added last to enable weighting by the other values
182        acc.add(spec, !bflag, tsys, inter, time);
183      }
184    }
185    const Vector<Bool>& msk = acc.getMask();
186    if ( allEQ(msk, False) ) {
187      uint n = rowstodelete.nelements();
188      rowstodelete.resize(n+1, True);
189      rowstodelete[n] = i;
190      continue;
191    }
192    //write out
193    Vector<uChar> flg(msk.shape());
194    convertArray(flg, !msk);
195    flagColOut.put(i, flg);
196    specColOut.put(i, acc.getSpectrum());
197    tsysColOut.put(i, acc.getTsys());
198    intColOut.put(i, acc.getInterval());
199    mjdColOut.put(i, acc.getTime());
200    // we should only have one cycle now -> reset it to be 0
201    // frequency switched data has different CYCLENO for different IFNO
202    // which requires resetting this value
203    cycColOut.put(i, uInt(0));
204    acc.reset();
205  }
206  if (rowstodelete.nelements() > 0) {
207    tout.removeRow(rowstodelete);
208    if (tout.nrow() == 0) {
209      throw(AipsError("Can't average fully flagged data."));
210    }
211  }
212  return out;
213}
214
215CountedPtr< Scantable >
216  STMath::averageChannel( const CountedPtr < Scantable > & in,
217                          const std::string & mode,
218                          const std::string& avmode )
219{
220  // clone as this is non insitu
221  bool insitu = insitu_;
222  setInsitu(false);
223  CountedPtr< Scantable > out = getScantable(in, true);
224  setInsitu(insitu);
225  Table& tout = out->table();
226  ArrayColumn<Float> specColOut(tout,"SPECTRA");
227  ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
228  ArrayColumn<Float> tsysColOut(tout,"TSYS");
229  ScalarColumn<uInt> scanColOut(tout,"SCANNO");
230  ScalarColumn<Double> intColOut(tout, "INTERVAL");
231  Table tmp = in->table().sort("BEAMNO");
232  Block<String> cols(3);
233  cols[0] = String("BEAMNO");
234  cols[1] = String("IFNO");
235  cols[2] = String("POLNO");
236  if ( avmode == "SCAN") {
237    cols.resize(4);
238    cols[3] = String("SCANNO");
239  }
240  uInt outrowCount = 0;
241  uChar userflag = 1 << 7;
242  TableIterator iter(tmp, cols);
243  while (!iter.pastEnd()) {
244    Table subt = iter.table();
245    ROArrayColumn<Float> specCol, tsysCol;
246    ROArrayColumn<uChar> flagCol;
247    ROScalarColumn<Double> intCol(subt, "INTERVAL");
248    specCol.attach(subt,"SPECTRA");
249    flagCol.attach(subt,"FLAGTRA");
250    tsysCol.attach(subt,"TSYS");
251    tout.addRow();
252    TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
253    if ( avmode != "SCAN") {
254      scanColOut.put(outrowCount, uInt(0));
255    }
256    Vector<Float> tmp;
257    specCol.get(0, tmp);
258    uInt nchan = tmp.nelements();
259    // have to do channel by channel here as MaskedArrMath
260    // doesn't have partialMedians
261    Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
262    Vector<Float> outspec(nchan);
263    Vector<uChar> outflag(nchan,0);
264    Vector<Float> outtsys(1);/// @fixme when tsys is channel based
265    for (uInt i=0; i<nchan; ++i) {
266      Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
267      MaskedArray<Float> ma = maskedArray(specs,flags);
268      outspec[i] = median(ma);
269      if ( allEQ(ma.getMask(), False) )
270        outflag[i] = userflag;// flag data
271    }
272    outtsys[0] = median(tsysCol.getColumn());
273    specColOut.put(outrowCount, outspec);
274    flagColOut.put(outrowCount, outflag);
275    tsysColOut.put(outrowCount, outtsys);
276    Double intsum = sum(intCol.getColumn());
277    intColOut.put(outrowCount, intsum);
278    ++outrowCount;
279    ++iter;
280  }
281  return out;
282}
283
284CountedPtr< Scantable > STMath::getScantable(const CountedPtr< Scantable >& in,
285                                             bool droprows)
286{
287  if (insitu_) return in;
288  else {
289    // clone
290    Scantable* tabp = new Scantable(*in, Bool(droprows));
291    return CountedPtr<Scantable>(tabp);
292  }
293}
294
295CountedPtr< Scantable > STMath::unaryOperate( const CountedPtr< Scantable >& in,
296                                              float val,
297                                              const std::string& mode,
298                                              bool tsys )
299{
300  CountedPtr< Scantable > out = getScantable(in, false);
301  Table& tab = out->table();
302  ArrayColumn<Float> specCol(tab,"SPECTRA");
303  ArrayColumn<Float> tsysCol(tab,"TSYS");
304  for (uInt i=0; i<tab.nrow(); ++i) {
305    Vector<Float> spec;
306    Vector<Float> ts;
307    specCol.get(i, spec);
308    tsysCol.get(i, ts);
309    if (mode == "MUL" || mode == "DIV") {
310      if (mode == "DIV") val = 1.0/val;
311      spec *= val;
312      specCol.put(i, spec);
313      if ( tsys ) {
314        ts *= val;
315        tsysCol.put(i, ts);
316      }
317    } else if ( mode == "ADD"  || mode == "SUB") {
318      if (mode == "SUB") val *= -1.0;
319      spec += val;
320      specCol.put(i, spec);
321      if ( tsys ) {
322        ts += val;
323        tsysCol.put(i, ts);
324      }
325    }
326  }
327  return out;
328}
329
330CountedPtr<Scantable> STMath::binaryOperate(const CountedPtr<Scantable>& left,
331                                            const CountedPtr<Scantable>& right,
332                                            const std::string& mode)
333{
334  bool insitu = insitu_;
335  if ( ! left->conformant(*right) ) {
336    throw(AipsError("'left' and 'right' scantables are not conformant."));
337  }
338  setInsitu(false);
339  CountedPtr< Scantable > out = getScantable(left, false);
340  setInsitu(insitu);
341  Table& tout = out->table();
342  Block<String> coln(5);
343  coln[0] = "SCANNO";  coln[1] = "CYCLENO";  coln[2] = "BEAMNO";
344  coln[3] = "IFNO";  coln[4] = "POLNO";
345  Table tmpl = tout.sort(coln);
346  Table tmpr = right->table().sort(coln);
347  ArrayColumn<Float> lspecCol(tmpl,"SPECTRA");
348  ROArrayColumn<Float> rspecCol(tmpr,"SPECTRA");
349  ArrayColumn<uChar> lflagCol(tmpl,"FLAGTRA");
350  ROArrayColumn<uChar> rflagCol(tmpr,"FLAGTRA");
351
352  for (uInt i=0; i<tout.nrow(); ++i) {
353    Vector<Float> lspecvec, rspecvec;
354    Vector<uChar> lflagvec, rflagvec;
355    lspecvec = lspecCol(i);    rspecvec = rspecCol(i);
356    lflagvec = lflagCol(i);    rflagvec = rflagCol(i);
357    MaskedArray<Float> mleft = maskedArray(lspecvec, lflagvec);
358    MaskedArray<Float> mright = maskedArray(rspecvec, rflagvec);
359    if (mode == "ADD") {
360      mleft += mright;
361    } else if ( mode == "SUB") {
362      mleft -= mright;
363    } else if ( mode == "MUL") {
364      mleft *= mright;
365    } else if ( mode == "DIV") {
366      mleft /= mright;
367    } else {
368      throw(AipsError("Illegal binary operator"));
369    }
370    lspecCol.put(i, mleft.getArray());
371  }
372  return out;
373}
374
375
376
377MaskedArray<Float> STMath::maskedArray( const Vector<Float>& s,
378                                        const Vector<uChar>& f)
379{
380  Vector<Bool> mask;
381  mask.resize(f.shape());
382  convertArray(mask, f);
383  return MaskedArray<Float>(s,!mask);
384}
385
386Vector<uChar> STMath::flagsFromMA(const MaskedArray<Float>& ma)
387{
388  const Vector<Bool>& m = ma.getMask();
389  Vector<uChar> flags(m.shape());
390  convertArray(flags, !m);
391  return flags;
392}
393
394CountedPtr< Scantable > STMath::autoQuotient( const CountedPtr< Scantable >& in,
395                                              const std::string & mode,
396                                              bool preserve )
397{
398  /// @todo make other modes available
399  /// modes should be "nearest", "pair"
400  // make this operation non insitu
401  const Table& tin = in->table();
402  Table ons = tin(tin.col("SRCTYPE") == Int(0));
403  Table offs = tin(tin.col("SRCTYPE") == Int(1));
404  if ( offs.nrow() == 0 )
405    throw(AipsError("No 'off' scans present."));
406  // put all "on" scans into output table
407
408  bool insitu = insitu_;
409  setInsitu(false);
410  CountedPtr< Scantable > out = getScantable(in, true);
411  setInsitu(insitu);
412  Table& tout = out->table();
413
414  TableCopy::copyRows(tout, ons);
415  TableRow row(tout);
416  ROScalarColumn<Double> offtimeCol(offs, "TIME");
417  ArrayColumn<Float> outspecCol(tout, "SPECTRA");
418  ROArrayColumn<Float> outtsysCol(tout, "TSYS");
419  ArrayColumn<uChar> outflagCol(tout, "FLAGTRA");
420  for (uInt i=0; i < tout.nrow(); ++i) {
421    const TableRecord& rec = row.get(i);
422    Double ontime = rec.asDouble("TIME");
423    Table presel = offs(offs.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
424                        && offs.col("IFNO") == Int(rec.asuInt("IFNO"))
425                        && offs.col("POLNO") == Int(rec.asuInt("POLNO")) );
426    ROScalarColumn<Double> offtimeCol(presel, "TIME");
427
428    Double mindeltat = min(abs(offtimeCol.getColumn() - ontime));
429    // Timestamp may vary within a cycle ???!!!
430    // increase this by 0.01 sec in case of rounding errors...
431    // There might be a better way to do this.
432    // fix to this fix. TIME is MJD, so 1.0d not 1.0s
433    mindeltat += 0.01/24./60./60.;
434    Table sel = presel( abs(presel.col("TIME")-ontime) <= mindeltat);
435
436    if ( sel.nrow() < 1 )  {
437      throw(AipsError("No closest in time found... This could be a rounding "
438                      "issue. Try quotient instead."));
439    }
440    TableRow offrow(sel);
441    const TableRecord& offrec = offrow.get(0);//should only be one row
442    RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
443    RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
444    RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
445    /// @fixme this assumes tsys is a scalar not vector
446    Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
447    Vector<Float> specon, tsyson;
448    outtsysCol.get(i, tsyson);
449    outspecCol.get(i, specon);
450    Vector<uChar> flagon;
451    outflagCol.get(i, flagon);
452    MaskedArray<Float> mon = maskedArray(specon, flagon);
453    MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
454    MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
455    if (preserve) {
456      quot -= tsysoffscalar;
457    } else {
458      quot -= tsyson[0];
459    }
460    outspecCol.put(i, quot.getArray());
461    outflagCol.put(i, flagsFromMA(quot));
462  }
463  // renumber scanno
464  TableIterator it(tout, "SCANNO");
465  uInt i = 0;
466  while ( !it.pastEnd() ) {
467    Table t = it.table();
468    TableVector<uInt> vec(t, "SCANNO");
469    vec = i;
470    ++i;
471    ++it;
472  }
473  return out;
474}
475
476
477CountedPtr< Scantable > STMath::quotient( const CountedPtr< Scantable > & on,
478                                          const CountedPtr< Scantable > & off,
479                                          bool preserve )
480{
481  bool insitu = insitu_;
482  if ( ! on->conformant(*off) ) {
483    throw(AipsError("'on' and 'off' scantables are not conformant."));
484  }
485  setInsitu(false);
486  CountedPtr< Scantable > out = getScantable(on, false);
487  setInsitu(insitu);
488  Table& tout = out->table();
489  const Table& toff = off->table();
490  TableIterator sit(tout, "SCANNO");
491  TableIterator s2it(toff, "SCANNO");
492  while ( !sit.pastEnd() ) {
493    Table ton = sit.table();
494    TableRow row(ton);
495    Table t = s2it.table();
496    ArrayColumn<Float> outspecCol(ton, "SPECTRA");
497    ROArrayColumn<Float> outtsysCol(ton, "TSYS");
498    ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
499    for (uInt i=0; i < ton.nrow(); ++i) {
500      const TableRecord& rec = row.get(i);
501      Table offsel = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
502                          && t.col("IFNO") == Int(rec.asuInt("IFNO"))
503                          && t.col("POLNO") == Int(rec.asuInt("POLNO")) );
504      if ( offsel.nrow() == 0 )
505        throw AipsError("STMath::quotient: no matching off");
506      TableRow offrow(offsel);
507      const TableRecord& offrec = offrow.get(0);//should be ncycles - take first
508      RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
509      RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
510      RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
511      Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
512      Vector<Float> specon, tsyson;
513      outtsysCol.get(i, tsyson);
514      outspecCol.get(i, specon);
515      Vector<uChar> flagon;
516      outflagCol.get(i, flagon);
517      MaskedArray<Float> mon = maskedArray(specon, flagon);
518      MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
519      MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
520      if (preserve) {
521        quot -= tsysoffscalar;
522      } else {
523        quot -= tsyson[0];
524      }
525      outspecCol.put(i, quot.getArray());
526      outflagCol.put(i, flagsFromMA(quot));
527    }
528    ++sit;
529    ++s2it;
530    // take the first off for each on scan which doesn't have a
531    // matching off scan
532    // non <= noff:  matching pairs, non > noff matching pairs then first off
533    if ( s2it.pastEnd() ) s2it.reset();
534  }
535  return out;
536}
537
538// dototalpower (migration of GBTIDL procedure dototalpower.pro)
539// calibrate the CAL on-off pair. It calculate Tsys and average CAL on-off subintegrations
540// do it for each cycles in a specific scan.
541CountedPtr< Scantable > STMath::dototalpower( const CountedPtr< Scantable >& calon,
542                                              const CountedPtr< Scantable >& caloff, Float tcal )
543{
544if ( ! calon->conformant(*caloff) ) {
545    throw(AipsError("'CAL on' and 'CAL off' scantables are not conformant."));
546  }
547  setInsitu(false);
548  CountedPtr< Scantable > out = getScantable(caloff, false);
549  Table& tout = out->table();
550  const Table& tcon = calon->table();
551  Vector<Float> tcalout;
552  Vector<Float> tcalout2;  //debug
553
554  if ( tout.nrow() != tcon.nrow() ) {
555    throw(AipsError("Mismatch in number of rows to form cal on - off pair."));
556  }
557  // iteration by scanno or cycle no.
558  TableIterator sit(tout, "SCANNO");
559  TableIterator s2it(tcon, "SCANNO");
560  while ( !sit.pastEnd() ) {
561    Table toff = sit.table();
562    TableRow row(toff);
563    Table t = s2it.table();
564    ScalarColumn<Double> outintCol(toff, "INTERVAL");
565    ArrayColumn<Float> outspecCol(toff, "SPECTRA");
566    ArrayColumn<Float> outtsysCol(toff, "TSYS");
567    ArrayColumn<uChar> outflagCol(toff, "FLAGTRA");
568    ROScalarColumn<uInt> outtcalIdCol(toff, "TCAL_ID");
569    ROScalarColumn<uInt> outpolCol(toff, "POLNO");
570    ROScalarColumn<Double> onintCol(t, "INTERVAL");
571    ROArrayColumn<Float> onspecCol(t, "SPECTRA");
572    ROArrayColumn<Float> ontsysCol(t, "TSYS");
573    ROArrayColumn<uChar> onflagCol(t, "FLAGTRA");
574    //ROScalarColumn<uInt> ontcalIdCol(t, "TCAL_ID");
575
576    for (uInt i=0; i < toff.nrow(); ++i) {
577      //skip these checks -> assumes the data order are the same between the cal on off pairs
578      //
579      Vector<Float> specCalon, specCaloff;
580      // to store scalar (mean) tsys
581      Vector<Float> tsysout(1);
582      uInt tcalId, polno;
583      Double offint, onint;
584      outpolCol.get(i, polno);
585      outspecCol.get(i, specCaloff);
586      onspecCol.get(i, specCalon);
587      Vector<uChar> flagCaloff, flagCalon;
588      outflagCol.get(i, flagCaloff);
589      onflagCol.get(i, flagCalon);
590      outtcalIdCol.get(i, tcalId);
591      outintCol.get(i, offint);
592      onintCol.get(i, onint);
593      // caluculate mean Tsys
594      uInt nchan = specCaloff.nelements();
595      // percentage of edge cut off
596      uInt pc = 10;
597      uInt bchan = nchan/pc;
598      uInt echan = nchan-bchan;
599
600      Slicer chansl(IPosition(1,bchan-1), IPosition(1,echan-1), IPosition(1,1),Slicer::endIsLast);
601      Vector<Float> testsubsp = specCaloff(chansl);
602      MaskedArray<Float> spoff = maskedArray( specCaloff(chansl),flagCaloff(chansl) );
603      MaskedArray<Float> spon = maskedArray( specCalon(chansl),flagCalon(chansl) );
604      MaskedArray<Float> spdiff = spon-spoff;
605      uInt noff = spoff.nelementsValid();
606      //uInt non = spon.nelementsValid();
607      uInt ndiff = spdiff.nelementsValid();
608      Float meantsys;
609
610/**
611      Double subspec, subdiff;
612      uInt usednchan;
613      subspec = 0;
614      subdiff = 0;
615      usednchan = 0;
616      for(uInt k=(bchan-1); k<echan; k++) {
617        subspec += specCaloff[k];
618        subdiff += static_cast<Double>(specCalon[k]-specCaloff[k]);
619        ++usednchan;
620      }
621**/
622      // get tcal if input tcal <= 0
623      String tcalt;
624      Float tcalUsed;
625      tcalUsed = tcal;
626      if ( tcal <= 0.0 ) {
627        caloff->tcal().getEntry(tcalt, tcalout, tcalId);
628        if (polno<=3) {
629          tcalUsed = tcalout[polno];
630        }
631        else {
632          tcalUsed = tcalout[0];
633        }
634      }
635
636      Float meanoff;
637      Float meandiff;
638      if (noff && ndiff) {
639         //Debug
640         //if(noff!=ndiff) cerr<<"noff and ndiff is not equal"<<endl;
641         meanoff = sum(spoff)/noff;
642         meandiff = sum(spdiff)/ndiff;
643         meantsys= (meanoff/meandiff )*tcalUsed + tcalUsed/2;
644      }
645      else {
646         meantsys=1;
647      }
648
649      tsysout[0] = Float(meantsys);
650      MaskedArray<Float> mcaloff = maskedArray(specCaloff, flagCaloff);
651      MaskedArray<Float> mcalon = maskedArray(specCalon, flagCalon);
652      MaskedArray<Float> sig =   Float(0.5) * (mcaloff + mcalon);
653      //uInt ncaloff = mcaloff.nelementsValid();
654      //uInt ncalon = mcalon.nelementsValid();
655
656      outintCol.put(i, offint+onint);
657      outspecCol.put(i, sig.getArray());
658      outflagCol.put(i, flagsFromMA(sig));
659      outtsysCol.put(i, tsysout);
660    }
661    ++sit;
662    ++s2it;
663  }
664  return out;
665}
666
667//dosigref - migrated from GBT IDL's dosigref.pro, do calibration of position switch
668// observatiions.
669// input: sig and ref scantables, and an optional boxcar smoothing width(default width=0,
670//        no smoothing).
671// output: resultant scantable [= (sig-ref/ref)*tsys]
672CountedPtr< Scantable > STMath::dosigref( const CountedPtr < Scantable >& sig,
673                                          const CountedPtr < Scantable >& ref,
674                                          int smoothref,
675                                          casa::Float tsysv,
676                                          casa::Float tau )
677{
678if ( ! ref->conformant(*sig) ) {
679    throw(AipsError("'sig' and 'ref' scantables are not conformant."));
680  }
681  setInsitu(false);
682  CountedPtr< Scantable > out = getScantable(sig, false);
683  CountedPtr< Scantable > smref;
684  if ( smoothref > 1 ) {
685    float fsmoothref = static_cast<float>(smoothref);
686    std::string inkernel = "boxcar";
687    smref = smooth(ref, inkernel, fsmoothref );
688    ostringstream oss;
689    oss<<"Applied smoothing of "<<fsmoothref<<" on the reference."<<endl;
690    pushLog(String(oss));
691  }
692  else {
693    smref = ref;
694  }
695  Table& tout = out->table();
696  const Table& tref = smref->table();
697  if ( tout.nrow() != tref.nrow() ) {
698    throw(AipsError("Mismatch in number of rows to form on-source and reference pair."));
699  }
700  // iteration by scanno? or cycle no.
701  TableIterator sit(tout, "SCANNO");
702  TableIterator s2it(tref, "SCANNO");
703  while ( !sit.pastEnd() ) {
704    Table ton = sit.table();
705    Table t = s2it.table();
706    ScalarColumn<Double> outintCol(ton, "INTERVAL");
707    ArrayColumn<Float> outspecCol(ton, "SPECTRA");
708    ArrayColumn<Float> outtsysCol(ton, "TSYS");
709    ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
710    ArrayColumn<Float> refspecCol(t, "SPECTRA");
711    ROScalarColumn<Double> refintCol(t, "INTERVAL");
712    ROArrayColumn<Float> reftsysCol(t, "TSYS");
713    ArrayColumn<uChar> refflagCol(t, "FLAGTRA");
714    ROScalarColumn<Float> refelevCol(t, "ELEVATION");
715    for (uInt i=0; i < ton.nrow(); ++i) {
716
717      Double onint, refint;
718      Vector<Float> specon, specref;
719      // to store scalar (mean) tsys
720      Vector<Float> tsysref;
721      outintCol.get(i, onint);
722      refintCol.get(i, refint);
723      outspecCol.get(i, specon);
724      refspecCol.get(i, specref);
725      Vector<uChar> flagref, flagon;
726      outflagCol.get(i, flagon);
727      refflagCol.get(i, flagref);
728      reftsysCol.get(i, tsysref);
729
730      Float tsysrefscalar;
731      if ( tsysv > 0.0 ) {
732        ostringstream oss;
733        Float elev;
734        refelevCol.get(i, elev);
735        oss << "user specified Tsys = " << tsysv;
736        // do recalc elevation if EL = 0
737        if ( elev == 0 ) {
738          throw(AipsError("EL=0, elevation data is missing."));
739        } else {
740          if ( tau <= 0.0 ) {
741            throw(AipsError("Valid tau is not supplied."));
742          } else {
743            tsysrefscalar = tsysv * exp(tau/elev);
744          }
745        }
746        oss << ", corrected (for El) tsys= "<<tsysrefscalar;
747        pushLog(String(oss));
748      }
749      else {
750        tsysrefscalar = tsysref[0];
751      }
752      //get quotient spectrum
753      MaskedArray<Float> mref = maskedArray(specref, flagref);
754      MaskedArray<Float> mon = maskedArray(specon, flagon);
755      MaskedArray<Float> specres =   tsysrefscalar*((mon - mref)/mref);
756      Double resint = onint*refint*smoothref/(onint+refint*smoothref);
757
758      //Debug
759      //cerr<<"Tsys used="<<tsysrefscalar<<endl;
760      // fill the result, replay signal tsys by reference tsys
761      outintCol.put(i, resint);
762      outspecCol.put(i, specres.getArray());
763      outflagCol.put(i, flagsFromMA(specres));
764      outtsysCol.put(i, tsysref);
765    }
766    ++sit;
767    ++s2it;
768  }
769  return out;
770}
771
772CountedPtr< Scantable > STMath::donod(const casa::CountedPtr<Scantable>& s,
773                                     const std::vector<int>& scans,
774                                     int smoothref,
775                                     casa::Float tsysv,
776                                     casa::Float tau,
777                                     casa::Float tcal )
778
779{
780  setInsitu(false);
781  STSelector sel;
782  std::vector<int> scan1, scan2, beams;
783  std::vector< vector<int> > scanpair;
784  std::vector<string> calstate;
785  String msg;
786
787  CountedPtr< Scantable > s1b1on, s1b1off, s1b2on, s1b2off;
788  CountedPtr< Scantable > s2b1on, s2b1off, s2b2on, s2b2off;
789
790  std::vector< CountedPtr< Scantable > > sctables;
791  sctables.push_back(s1b1on);
792  sctables.push_back(s1b1off);
793  sctables.push_back(s1b2on);
794  sctables.push_back(s1b2off);
795  sctables.push_back(s2b1on);
796  sctables.push_back(s2b1off);
797  sctables.push_back(s2b2on);
798  sctables.push_back(s2b2off);
799
800  //check scanlist
801  int n=s->checkScanInfo(scans);
802  if (n==1) {
803     throw(AipsError("Incorrect scan pairs. "));
804  }
805
806  // Assume scans contain only a pair of consecutive scan numbers.
807  // It is assumed that first beam, b1,  is on target.
808  // There is no check if the first beam is on or not.
809  if ( scans.size()==1 ) {
810    scan1.push_back(scans[0]);
811    scan2.push_back(scans[0]+1);
812  } else if ( scans.size()==2 ) {
813   scan1.push_back(scans[0]);
814   scan2.push_back(scans[1]);
815  } else {
816    if ( scans.size()%2 == 0 ) {
817      for (uInt i=0; i<scans.size(); i++) {
818        if (i%2 == 0) {
819          scan1.push_back(scans[i]);
820        }
821        else {
822          scan2.push_back(scans[i]);
823        }
824      }
825    } else {
826      throw(AipsError("Odd numbers of scans, cannot form pairs."));
827    }
828  }
829  scanpair.push_back(scan1);
830  scanpair.push_back(scan2);
831  calstate.push_back("*calon");
832  calstate.push_back("*[^calon]");
833  CountedPtr< Scantable > ws = getScantable(s, false);
834  uInt l=0;
835  while ( l < sctables.size() ) {
836    for (uInt i=0; i < 2; i++) {
837      for (uInt j=0; j < 2; j++) {
838        for (uInt k=0; k < 2; k++) {
839          sel.reset();
840          sel.setScans(scanpair[i]);
841          sel.setName(calstate[k]);
842          beams.clear();
843          beams.push_back(j);
844          sel.setBeams(beams);
845          ws->setSelection(sel);
846          sctables[l]= getScantable(ws, false);
847          l++;
848        }
849      }
850    }
851  }
852
853  // replace here by splitData or getData functionality
854  CountedPtr< Scantable > sig1;
855  CountedPtr< Scantable > ref1;
856  CountedPtr< Scantable > sig2;
857  CountedPtr< Scantable > ref2;
858  CountedPtr< Scantable > calb1;
859  CountedPtr< Scantable > calb2;
860
861  msg=String("Processing dototalpower for subset of the data");
862  ostringstream oss1;
863  oss1 << msg  << endl;
864  pushLog(String(oss1));
865  // Debug for IRC CS data
866  //float tcal1=7.0;
867  //float tcal2=4.0;
868  sig1 = dototalpower(sctables[0], sctables[1], tcal=tcal);
869  ref1 = dototalpower(sctables[2], sctables[3], tcal=tcal);
870  ref2 = dototalpower(sctables[4], sctables[5], tcal=tcal);
871  sig2 = dototalpower(sctables[6], sctables[7], tcal=tcal);
872
873  // correction of user-specified tsys for elevation here
874
875  // dosigref calibration
876  msg=String("Processing dosigref for subset of the data");
877  ostringstream oss2;
878  oss2 << msg  << endl;
879  pushLog(String(oss2));
880  calb1=dosigref(sig1,ref2,smoothref,tsysv,tau);
881  calb2=dosigref(sig2,ref1,smoothref,tsysv,tau);
882
883  // iteration by scanno or cycle no.
884  Table& tcalb1 = calb1->table();
885  Table& tcalb2 = calb2->table();
886  TableIterator sit(tcalb1, "SCANNO");
887  TableIterator s2it(tcalb2, "SCANNO");
888  while ( !sit.pastEnd() ) {
889    Table t1 = sit.table();
890    Table t2= s2it.table();
891    ArrayColumn<Float> outspecCol(t1, "SPECTRA");
892    ArrayColumn<Float> outtsysCol(t1, "TSYS");
893    ArrayColumn<uChar> outflagCol(t1, "FLAGTRA");
894    ScalarColumn<Double> outintCol(t1, "INTERVAL");
895    ArrayColumn<Float> t2specCol(t2, "SPECTRA");
896    ROArrayColumn<Float> t2tsysCol(t2, "TSYS");
897    ArrayColumn<uChar> t2flagCol(t2, "FLAGTRA");
898    ROScalarColumn<Double> t2intCol(t2, "INTERVAL");
899    for (uInt i=0; i < t1.nrow(); ++i) {
900      Vector<Float> spec1, spec2;
901      // to store scalar (mean) tsys
902      Vector<Float> tsys1, tsys2;
903      Vector<uChar> flag1, flag2;
904      Double tint1, tint2;
905      outspecCol.get(i, spec1);
906      t2specCol.get(i, spec2);
907      outflagCol.get(i, flag1);
908      t2flagCol.get(i, flag2);
909      outtsysCol.get(i, tsys1);
910      t2tsysCol.get(i, tsys2);
911      outintCol.get(i, tint1);
912      t2intCol.get(i, tint2);
913      // average
914      // assume scalar tsys for weights
915      Float wt1, wt2, tsyssq1, tsyssq2;
916      tsyssq1 = tsys1[0]*tsys1[0];
917      tsyssq2 = tsys2[0]*tsys2[0];
918      wt1 = Float(tint1)/tsyssq1;
919      wt2 = Float(tint2)/tsyssq2;
920      Float invsumwt=1/(wt1+wt2);
921      MaskedArray<Float> mspec1 = maskedArray(spec1, flag1);
922      MaskedArray<Float> mspec2 = maskedArray(spec2, flag2);
923      MaskedArray<Float> avspec =  invsumwt * (wt1*mspec1 + wt2*mspec2);
924      //Array<Float> avtsys =  Float(0.5) * (tsys1 + tsys2);
925      // cerr<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<endl;
926      tsys1[0] = sqrt(tsyssq1 + tsyssq2);
927      Array<Float> avtsys =  tsys1;
928
929      outspecCol.put(i, avspec.getArray());
930      outflagCol.put(i, flagsFromMA(avspec));
931      outtsysCol.put(i, avtsys);
932    }
933    ++sit;
934    ++s2it;
935  }
936  return calb1;
937}
938
939//GBTIDL version of frequency switched data calibration
940CountedPtr< Scantable > STMath::dofs( const CountedPtr< Scantable >& s,
941                                      const std::vector<int>& scans,
942                                      int smoothref,
943                                      casa::Float tsysv,
944                                      casa::Float tau,
945                                      casa::Float tcal )
946{
947
948
949  STSelector sel;
950  CountedPtr< Scantable > ws = getScantable(s, false);
951  CountedPtr< Scantable > sig, sigwcal, ref, refwcal;
952  CountedPtr< Scantable > calsig, calref, out;
953
954  //split the data
955  sel.setName("*_fs");
956  ws->setSelection(sel);
957  sig = getScantable(ws,false);
958  sel.reset();
959  sel.setName("*_fs_calon");
960  ws->setSelection(sel);
961  sigwcal = getScantable(ws,false);
962  sel.reset();
963  sel.setName("*_fsr");
964  ws->setSelection(sel);
965  ref = getScantable(ws,false);
966  sel.reset();
967  sel.setName("*_fsr_calon");
968  ws->setSelection(sel);
969  refwcal = getScantable(ws,false);
970
971  calsig = dototalpower(sigwcal, sig, tcal=tcal);
972  calref = dototalpower(refwcal, ref, tcal=tcal);
973
974  out=dosigref(calsig,calref,smoothref,tsysv,tau);
975
976  return out;
977}
978
979
980CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
981{
982  // make copy or reference
983  CountedPtr< Scantable > out = getScantable(in, false);
984  Table& tout = out->table();
985  Block<String> cols(4);
986  cols[0] = String("SCANNO");
987  cols[1] = String("CYCLENO");
988  cols[2] = String("BEAMNO");
989  cols[3] = String("POLNO");
990  TableIterator iter(tout, cols);
991  while (!iter.pastEnd()) {
992    Table subt = iter.table();
993    // this should leave us with two rows for the two IFs....if not ignore
994    if (subt.nrow() != 2 ) {
995      continue;
996    }
997    ArrayColumn<Float> specCol(subt, "SPECTRA");
998    ArrayColumn<Float> tsysCol(subt, "TSYS");
999    ArrayColumn<uChar> flagCol(subt, "FLAGTRA");
1000    Vector<Float> onspec,offspec, ontsys, offtsys;
1001    Vector<uChar> onflag, offflag;
1002    tsysCol.get(0, ontsys);   tsysCol.get(1, offtsys);
1003    specCol.get(0, onspec);   specCol.get(1, offspec);
1004    flagCol.get(0, onflag);   flagCol.get(1, offflag);
1005    MaskedArray<Float> on  = maskedArray(onspec, onflag);
1006    MaskedArray<Float> off = maskedArray(offspec, offflag);
1007    MaskedArray<Float> oncopy = on.copy();
1008
1009    on /= off; on -= 1.0f;
1010    on *= ontsys[0];
1011    off /= oncopy; off -= 1.0f;
1012    off *= offtsys[0];
1013    specCol.put(0, on.getArray());
1014    const Vector<Bool>& m0 = on.getMask();
1015    Vector<uChar> flags0(m0.shape());
1016    convertArray(flags0, !m0);
1017    flagCol.put(0, flags0);
1018
1019    specCol.put(1, off.getArray());
1020    const Vector<Bool>& m1 = off.getMask();
1021    Vector<uChar> flags1(m1.shape());
1022    convertArray(flags1, !m1);
1023    flagCol.put(1, flags1);
1024    ++iter;
1025  }
1026
1027  return out;
1028}
1029
1030std::vector< float > STMath::statistic( const CountedPtr< Scantable > & in,
1031                                        const std::vector< bool > & mask,
1032                                        const std::string& which )
1033{
1034
1035  Vector<Bool> m(mask);
1036  const Table& tab = in->table();
1037  ROArrayColumn<Float> specCol(tab, "SPECTRA");
1038  ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1039  std::vector<float> out;
1040  for (uInt i=0; i < tab.nrow(); ++i ) {
1041    Vector<Float> spec; specCol.get(i, spec);
1042    Vector<uChar> flag; flagCol.get(i, flag);
1043    MaskedArray<Float> ma  = maskedArray(spec, flag);
1044    float outstat = 0.0;
1045    if ( spec.nelements() == m.nelements() ) {
1046      outstat = mathutil::statistics(which, ma(m));
1047    } else {
1048      outstat = mathutil::statistics(which, ma);
1049    }
1050    out.push_back(outstat);
1051  }
1052  return out;
1053}
1054
1055CountedPtr< Scantable > STMath::bin( const CountedPtr< Scantable > & in,
1056                                     int width )
1057{
1058  if ( !in->getSelection().empty() ) throw(AipsError("Can't bin subset of the data."));
1059  CountedPtr< Scantable > out = getScantable(in, false);
1060  Table& tout = out->table();
1061  out->frequencies().rescale(width, "BIN");
1062  ArrayColumn<Float> specCol(tout, "SPECTRA");
1063  ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
1064  for (uInt i=0; i < tout.nrow(); ++i ) {
1065    MaskedArray<Float> main  = maskedArray(specCol(i), flagCol(i));
1066    MaskedArray<Float> maout;
1067    LatticeUtilities::bin(maout, main, 0, Int(width));
1068    /// @todo implement channel based tsys binning
1069    specCol.put(i, maout.getArray());
1070    flagCol.put(i, flagsFromMA(maout));
1071    // take only the first binned spectrum's length for the deprecated
1072    // global header item nChan
1073    if (i==0) tout.rwKeywordSet().define(String("nChan"),
1074                                       Int(maout.getArray().nelements()));
1075  }
1076  return out;
1077}
1078
1079CountedPtr< Scantable > STMath::resample( const CountedPtr< Scantable >& in,
1080                                          const std::string& method,
1081                                          float width )
1082//
1083// Should add the possibility of width being specified in km/s. This means
1084// that for each freqID (SpectralCoordinate) we will need to convert to an
1085// average channel width (say at the reference pixel).  Then we would need
1086// to be careful to make sure each spectrum (of different freqID)
1087// is the same length.
1088//
1089{
1090  //InterpolateArray1D<Double,Float>::InterpolationMethod interp;
1091  Int interpMethod(stringToIMethod(method));
1092
1093  CountedPtr< Scantable > out = getScantable(in, false);
1094  Table& tout = out->table();
1095
1096// Resample SpectralCoordinates (one per freqID)
1097  out->frequencies().rescale(width, "RESAMPLE");
1098  TableIterator iter(tout, "IFNO");
1099  TableRow row(tout);
1100  while ( !iter.pastEnd() ) {
1101    Table tab = iter.table();
1102    ArrayColumn<Float> specCol(tab, "SPECTRA");
1103    //ArrayColumn<Float> tsysCol(tout, "TSYS");
1104    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1105    Vector<Float> spec;
1106    Vector<uChar> flag;
1107    specCol.get(0,spec); // the number of channels should be constant per IF
1108    uInt nChanIn = spec.nelements();
1109    Vector<Float> xIn(nChanIn); indgen(xIn);
1110    Int fac =  Int(nChanIn/width);
1111    Vector<Float> xOut(fac+10); // 10 to be safe - resize later
1112    uInt k = 0;
1113    Float x = 0.0;
1114    while (x < Float(nChanIn) ) {
1115      xOut(k) = x;
1116      k++;
1117      x += width;
1118    }
1119    uInt nChanOut = k;
1120    xOut.resize(nChanOut, True);
1121    // process all rows for this IFNO
1122    Vector<Float> specOut;
1123    Vector<Bool> maskOut;
1124    Vector<uChar> flagOut;
1125    for (uInt i=0; i < tab.nrow(); ++i) {
1126      specCol.get(i, spec);
1127      flagCol.get(i, flag);
1128      Vector<Bool> mask(flag.nelements());
1129      convertArray(mask, flag);
1130
1131      IPosition shapeIn(spec.shape());
1132      //sh.nchan = nChanOut;
1133      InterpolateArray1D<Float,Float>::interpolate(specOut, maskOut, xOut,
1134                                                   xIn, spec, mask,
1135                                                   interpMethod, True, True);
1136      /// @todo do the same for channel based Tsys
1137      flagOut.resize(maskOut.nelements());
1138      convertArray(flagOut, maskOut);
1139      specCol.put(i, specOut);
1140      flagCol.put(i, flagOut);
1141    }
1142    ++iter;
1143  }
1144
1145  return out;
1146}
1147
1148STMath::imethod STMath::stringToIMethod(const std::string& in)
1149{
1150  static STMath::imap lookup;
1151
1152  // initialize the lookup table if necessary
1153  if ( lookup.empty() ) {
1154    lookup["nearest"]   = InterpolateArray1D<Double,Float>::nearestNeighbour;
1155    lookup["linear"] = InterpolateArray1D<Double,Float>::linear;
1156    lookup["cubic"]  = InterpolateArray1D<Double,Float>::cubic;
1157    lookup["spline"]  = InterpolateArray1D<Double,Float>::spline;
1158  }
1159
1160  STMath::imap::const_iterator iter = lookup.find(in);
1161
1162  if ( lookup.end() == iter ) {
1163    std::string message = in;
1164    message += " is not a valid interpolation mode";
1165    throw(AipsError(message));
1166  }
1167  return iter->second;
1168}
1169
1170WeightType STMath::stringToWeight(const std::string& in)
1171{
1172  static std::map<std::string, WeightType> lookup;
1173
1174  // initialize the lookup table if necessary
1175  if ( lookup.empty() ) {
1176    lookup["NONE"]   = asap::NONE;
1177    lookup["TINT"] = asap::TINT;
1178    lookup["TINTSYS"]  = asap::TINTSYS;
1179    lookup["TSYS"]  = asap::TSYS;
1180    lookup["VAR"]  = asap::VAR;
1181  }
1182
1183  std::map<std::string, WeightType>::const_iterator iter = lookup.find(in);
1184
1185  if ( lookup.end() == iter ) {
1186    std::string message = in;
1187    message += " is not a valid weighting mode";
1188    throw(AipsError(message));
1189  }
1190  return iter->second;
1191}
1192
1193CountedPtr< Scantable > STMath::gainElevation( const CountedPtr< Scantable >& in,
1194                                               const vector< float > & coeff,
1195                                               const std::string & filename,
1196                                               const std::string& method)
1197{
1198  // Get elevation data from Scantable and convert to degrees
1199  CountedPtr< Scantable > out = getScantable(in, false);
1200  Table& tab = out->table();
1201  ROScalarColumn<Float> elev(tab, "ELEVATION");
1202  Vector<Float> x = elev.getColumn();
1203  x *= Float(180 / C::pi);                        // Degrees
1204
1205  Vector<Float> coeffs(coeff);
1206  const uInt nc = coeffs.nelements();
1207  if ( filename.length() > 0 && nc > 0 ) {
1208    throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
1209  }
1210
1211  // Correct
1212  if ( nc > 0 || filename.length() == 0 ) {
1213    // Find instrument
1214    Bool throwit = True;
1215    Instrument inst =
1216      STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
1217                                throwit);
1218
1219    // Set polynomial
1220    Polynomial<Float>* ppoly = 0;
1221    Vector<Float> coeff;
1222    String msg;
1223    if ( nc > 0 ) {
1224      ppoly = new Polynomial<Float>(nc);
1225      coeff = coeffs;
1226      msg = String("user");
1227    } else {
1228      STAttr sdAttr;
1229      coeff = sdAttr.gainElevationPoly(inst);
1230      ppoly = new Polynomial<Float>(3);
1231      msg = String("built in");
1232    }
1233
1234    if ( coeff.nelements() > 0 ) {
1235      ppoly->setCoefficients(coeff);
1236    } else {
1237      delete ppoly;
1238      throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
1239    }
1240    ostringstream oss;
1241    oss << "Making polynomial correction with " << msg << " coefficients:" << endl;
1242    oss << "   " <<  coeff;
1243    pushLog(String(oss));
1244    const uInt nrow = tab.nrow();
1245    Vector<Float> factor(nrow);
1246    for ( uInt i=0; i < nrow; ++i ) {
1247      factor[i] = 1.0 / (*ppoly)(x[i]);
1248    }
1249    delete ppoly;
1250    scaleByVector(tab, factor, true);
1251
1252  } else {
1253    // Read and correct
1254    pushLog("Making correction from ascii Table");
1255    scaleFromAsciiTable(tab, filename, method, x, true);
1256  }
1257  return out;
1258}
1259
1260void STMath::scaleFromAsciiTable(Table& in, const std::string& filename,
1261                                 const std::string& method,
1262                                 const Vector<Float>& xout, bool dotsys)
1263{
1264
1265// Read gain-elevation ascii file data into a Table.
1266
1267  String formatString;
1268  Table tbl = readAsciiTable(formatString, Table::Memory, filename, "", "", False);
1269  scaleFromTable(in, tbl, method, xout, dotsys);
1270}
1271
1272void STMath::scaleFromTable(Table& in,
1273                            const Table& table,
1274                            const std::string& method,
1275                            const Vector<Float>& xout, bool dotsys)
1276{
1277
1278  ROScalarColumn<Float> geElCol(table, "ELEVATION");
1279  ROScalarColumn<Float> geFacCol(table, "FACTOR");
1280  Vector<Float> xin = geElCol.getColumn();
1281  Vector<Float> yin = geFacCol.getColumn();
1282  Vector<Bool> maskin(xin.nelements(),True);
1283
1284  // Interpolate (and extrapolate) with desired method
1285
1286  InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
1287
1288   Vector<Float> yout;
1289   Vector<Bool> maskout;
1290   InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout,
1291                                                xin, yin, maskin, interp,
1292                                                True, True);
1293
1294   scaleByVector(in, Float(1.0)/yout, dotsys);
1295}
1296
1297void STMath::scaleByVector( Table& in,
1298                            const Vector< Float >& factor,
1299                            bool dotsys )
1300{
1301  uInt nrow = in.nrow();
1302  if ( factor.nelements() != nrow ) {
1303    throw(AipsError("factors.nelements() != table.nelements()"));
1304  }
1305  ArrayColumn<Float> specCol(in, "SPECTRA");
1306  ArrayColumn<uChar> flagCol(in, "FLAGTRA");
1307  ArrayColumn<Float> tsysCol(in, "TSYS");
1308  for (uInt i=0; i < nrow; ++i) {
1309    MaskedArray<Float> ma  = maskedArray(specCol(i), flagCol(i));
1310    ma *= factor[i];
1311    specCol.put(i, ma.getArray());
1312    flagCol.put(i, flagsFromMA(ma));
1313    if ( dotsys ) {
1314      Vector<Float> tsys = tsysCol(i);
1315      tsys *= factor[i];
1316      tsysCol.put(i,tsys);
1317    }
1318  }
1319}
1320
1321CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in,
1322                                             float d, float etaap,
1323                                             float jyperk )
1324{
1325  CountedPtr< Scantable > out = getScantable(in, false);
1326  Table& tab = in->table();
1327  Unit fluxUnit(tab.keywordSet().asString("FluxUnit"));
1328  Unit K(String("K"));
1329  Unit JY(String("Jy"));
1330
1331  bool tokelvin = true;
1332  Double cfac = 1.0;
1333
1334  if ( fluxUnit == JY ) {
1335    pushLog("Converting to K");
1336    Quantum<Double> t(1.0,fluxUnit);
1337    Quantum<Double> t2 = t.get(JY);
1338    cfac = (t2 / t).getValue();               // value to Jy
1339
1340    tokelvin = true;
1341    out->setFluxUnit("K");
1342  } else if ( fluxUnit == K ) {
1343    pushLog("Converting to Jy");
1344    Quantum<Double> t(1.0,fluxUnit);
1345    Quantum<Double> t2 = t.get(K);
1346    cfac = (t2 / t).getValue();              // value to K
1347
1348    tokelvin = false;
1349    out->setFluxUnit("Jy");
1350  } else {
1351    throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
1352  }
1353  // Make sure input values are converted to either Jy or K first...
1354  Float factor = cfac;
1355
1356  // Select method
1357  if (jyperk > 0.0) {
1358    factor *= jyperk;
1359    if ( tokelvin ) factor = 1.0 / jyperk;
1360    ostringstream oss;
1361    oss << "Jy/K = " << jyperk;
1362    pushLog(String(oss));
1363    Vector<Float> factors(tab.nrow(), factor);
1364    scaleByVector(tab,factors, false);
1365  } else if ( etaap > 0.0) {
1366    if (d < 0) {
1367      Instrument inst =
1368        STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
1369                                  True);
1370      STAttr sda;
1371      d = sda.diameter(inst);
1372    }
1373    jyperk = STAttr::findJyPerK(etaap, d);
1374    ostringstream oss;
1375    oss << "Jy/K = " << jyperk;
1376    pushLog(String(oss));
1377    factor *= jyperk;
1378    if ( tokelvin ) {
1379      factor = 1.0 / factor;
1380    }
1381    Vector<Float> factors(tab.nrow(), factor);
1382    scaleByVector(tab, factors, False);
1383  } else {
1384
1385    // OK now we must deal with automatic look up of values.
1386    // We must also deal with the fact that the factors need
1387    // to be computed per IF and may be different and may
1388    // change per integration.
1389
1390    pushLog("Looking up conversion factors");
1391    convertBrightnessUnits(out, tokelvin, cfac);
1392  }
1393
1394  return out;
1395}
1396
1397void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in,
1398                                     bool tokelvin, float cfac )
1399{
1400  Table& table = in->table();
1401  Instrument inst =
1402    STAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True);
1403  TableIterator iter(table, "FREQ_ID");
1404  STFrequencies stfreqs = in->frequencies();
1405  STAttr sdAtt;
1406  while (!iter.pastEnd()) {
1407    Table tab = iter.table();
1408    ArrayColumn<Float> specCol(tab, "SPECTRA");
1409    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1410    ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID");
1411    MEpoch::ROScalarColumn timeCol(tab, "TIME");
1412
1413    uInt freqid; freqidCol.get(0, freqid);
1414    Vector<Float> tmpspec; specCol.get(0, tmpspec);
1415    // STAttr.JyPerK has a Vector interface... change sometime.
1416    Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements()));
1417    for ( uInt i=0; i<tab.nrow(); ++i) {
1418      Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0];
1419      Float factor = cfac * jyperk;
1420      if ( tokelvin ) factor = Float(1.0) / factor;
1421      MaskedArray<Float> ma  = maskedArray(specCol(i), flagCol(i));
1422      ma *= factor;
1423      specCol.put(i, ma.getArray());
1424      flagCol.put(i, flagsFromMA(ma));
1425    }
1426  ++iter;
1427  }
1428}
1429
1430CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
1431                                         float tau )
1432{
1433  CountedPtr< Scantable > out = getScantable(in, false);
1434
1435  Table tab = out->table();
1436  ROScalarColumn<Float> elev(tab, "ELEVATION");
1437  ArrayColumn<Float> specCol(tab, "SPECTRA");
1438  ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1439  ArrayColumn<Float> tsysCol(tab, "TSYS");
1440  for ( uInt i=0; i<tab.nrow(); ++i) {
1441    Float zdist = Float(C::pi_2) - elev(i);
1442    Float factor = exp(tau/cos(zdist));
1443    MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
1444    ma *= factor;
1445    specCol.put(i, ma.getArray());
1446    flagCol.put(i, flagsFromMA(ma));
1447    Vector<Float> tsys;
1448    tsysCol.get(i, tsys);
1449    tsys *= factor;
1450    tsysCol.put(i, tsys);
1451  }
1452  return out;
1453}
1454
1455CountedPtr< Scantable > STMath::smoothOther( const CountedPtr< Scantable >& in,
1456                                             const std::string& kernel,
1457                                             float width )
1458{
1459  CountedPtr< Scantable > out = getScantable(in, false);
1460  Table& table = out->table();
1461  ArrayColumn<Float> specCol(table, "SPECTRA");
1462  ArrayColumn<uChar> flagCol(table, "FLAGTRA");
1463  Vector<Float> spec;
1464  Vector<uChar> flag;
1465  for ( uInt i=0; i<table.nrow(); ++i) {
1466    specCol.get(i, spec);
1467    flagCol.get(i, flag);
1468    Vector<Bool> mask(flag.nelements());
1469    convertArray(mask, flag);
1470    Vector<Float> specout;
1471    Vector<Bool> maskout;
1472    if ( kernel == "hanning" ) {
1473      mathutil::hanning(specout, maskout, spec , !mask);
1474      convertArray(flag, !maskout);
1475    } else if (  kernel == "rmedian" ) {
1476      mathutil::runningMedian(specout, maskout, spec , mask, width);
1477      convertArray(flag, maskout);
1478    }
1479    flagCol.put(i, flag);
1480    specCol.put(i, specout);
1481  }
1482  return out;
1483}
1484
1485CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
1486                                        const std::string& kernel, float width )
1487{
1488  if (kernel == "rmedian"  || kernel == "hanning") {
1489    return smoothOther(in, kernel, width);
1490  }
1491  CountedPtr< Scantable > out = getScantable(in, false);
1492  Table& table = out->table();
1493  VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
1494  // same IFNO should have same no of channels
1495  // this saves overhead
1496  TableIterator iter(table, "IFNO");
1497  while (!iter.pastEnd()) {
1498    Table tab = iter.table();
1499    ArrayColumn<Float> specCol(tab, "SPECTRA");
1500    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1501    Vector<Float> tmpspec; specCol.get(0, tmpspec);
1502    uInt nchan = tmpspec.nelements();
1503    Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
1504    Convolver<Float> conv(kvec, IPosition(1,nchan));
1505    Vector<Float> spec;
1506    Vector<uChar> flag;
1507    for ( uInt i=0; i<tab.nrow(); ++i) {
1508      specCol.get(i, spec);
1509      flagCol.get(i, flag);
1510      Vector<Bool> mask(flag.nelements());
1511      convertArray(mask, flag);
1512      Vector<Float> specout;
1513      mathutil::replaceMaskByZero(specout, mask);
1514      conv.linearConv(specout, spec);
1515      specCol.put(i, specout);
1516    }
1517    ++iter;
1518  }
1519  return out;
1520}
1521
1522CountedPtr< Scantable >
1523  STMath::merge( const std::vector< CountedPtr < Scantable > >& in )
1524{
1525  if ( in.size() < 2 ) {
1526    throw(AipsError("Need at least two scantables to perform a merge."));
1527  }
1528  std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
1529  bool insitu = insitu_;
1530  setInsitu(false);
1531  CountedPtr< Scantable > out = getScantable(*it, false);
1532  setInsitu(insitu);
1533  Table& tout = out->table();
1534  ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
1535  ScalarColumn<uInt> scannocol(tout,"SCANNO"), focusidcol(tout,"FOCUS_ID");
1536  // Renumber SCANNO to be 0-based
1537  Vector<uInt> scannos = scannocol.getColumn();
1538  uInt offset = min(scannos);
1539  scannos -= offset;
1540  scannocol.putColumn(scannos);
1541  uInt newscanno = max(scannos)+1;
1542  ++it;
1543  while ( it != in.end() ){
1544    if ( ! (*it)->conformant(*out) ) {
1545      // non conformant.
1546      pushLog(String("Warning: Can't merge scantables as header info differs."));
1547    }
1548    out->appendToHistoryTable((*it)->history());
1549    const Table& tab = (*it)->table();
1550    TableIterator scanit(tab, "SCANNO");
1551    while (!scanit.pastEnd()) {
1552      TableIterator freqit(scanit.table(), "FREQ_ID");
1553      while ( !freqit.pastEnd() ) {
1554        Table thetab = freqit.table();
1555        uInt nrow = tout.nrow();
1556        //tout.addRow(thetab.nrow());
1557        TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
1558        ROTableRow row(thetab);
1559        for ( uInt i=0; i<thetab.nrow(); ++i) {
1560          uInt k = nrow+i;
1561          scannocol.put(k, newscanno);
1562          const TableRecord& rec = row.get(i);
1563          Double rv,rp,inc;
1564          (*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
1565          uInt id;
1566          id = out->frequencies().addEntry(rp, rv, inc);
1567          freqidcol.put(k,id);
1568          String name,fname;Double rf;
1569          (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
1570          id = out->molecules().addEntry(rf, name, fname);
1571          molidcol.put(k, id);
1572          Float frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
1573          (*it)->focus().getEntry(fax, ftan, frot, fhand,
1574                                  fmount,fuser, fxy, fxyp,
1575                                  rec.asuInt("FOCUS_ID"));
1576          id = out->focus().addEntry(fax, ftan, frot, fhand,
1577                                     fmount,fuser, fxy, fxyp);
1578          focusidcol.put(k, id);
1579        }
1580        ++freqit;
1581      }
1582      ++newscanno;
1583      ++scanit;
1584    }
1585    ++it;
1586  }
1587  return out;
1588}
1589
1590CountedPtr< Scantable >
1591  STMath::invertPhase( const CountedPtr < Scantable >& in )
1592{
1593  return applyToPol(in, &STPol::invertPhase, Float(0.0));
1594}
1595
1596CountedPtr< Scantable >
1597  STMath::rotateXYPhase( const CountedPtr < Scantable >& in, float phase )
1598{
1599   return applyToPol(in, &STPol::rotatePhase, Float(phase));
1600}
1601
1602CountedPtr< Scantable >
1603  STMath::rotateLinPolPhase( const CountedPtr < Scantable >& in, float phase )
1604{
1605  return applyToPol(in, &STPol::rotateLinPolPhase, Float(phase));
1606}
1607
1608CountedPtr< Scantable > STMath::applyToPol( const CountedPtr<Scantable>& in,
1609                                             STPol::polOperation fptr,
1610                                             Float phase )
1611{
1612  CountedPtr< Scantable > out = getScantable(in, false);
1613  Table& tout = out->table();
1614  Block<String> cols(4);
1615  cols[0] = String("SCANNO");
1616  cols[1] = String("BEAMNO");
1617  cols[2] = String("IFNO");
1618  cols[3] = String("CYCLENO");
1619  TableIterator iter(tout, cols);
1620  CountedPtr<STPol> stpol = STPol::getPolClass(out->factories_,
1621                                               out->getPolType() );
1622  while (!iter.pastEnd()) {
1623    Table t = iter.table();
1624    ArrayColumn<Float> speccol(t, "SPECTRA");
1625    ScalarColumn<uInt> focidcol(t, "FOCUS_ID");
1626    ScalarColumn<Float> parancol(t, "PARANGLE");
1627    Matrix<Float> pols(speccol.getColumn());
1628    try {
1629      stpol->setSpectra(pols);
1630      Float fang,fhand,parang;
1631      fang = in->focusTable_.getTotalFeedAngle(focidcol(0));
1632      fhand = in->focusTable_.getFeedHand(focidcol(0));
1633      parang = parancol(0);
1634      /// @todo re-enable this
1635      // disable total feed angle to support paralactifying Caswell style
1636      stpol->setPhaseCorrections(parang, -parang, fhand);
1637      // use a member function pointer in STPol.  This only works on
1638      // the STPol pointer itself, not the Counted Pointer so
1639      // derefernce it.
1640      (&(*(stpol))->*fptr)(phase);
1641      speccol.putColumn(stpol->getSpectra());
1642    } catch (AipsError& e) {
1643      //delete stpol;stpol=0;
1644      throw(e);
1645    }
1646    ++iter;
1647  }
1648  //delete stpol;stpol=0;
1649  return out;
1650}
1651
1652CountedPtr< Scantable >
1653  STMath::swapPolarisations( const CountedPtr< Scantable > & in )
1654{
1655  CountedPtr< Scantable > out = getScantable(in, false);
1656  Table& tout = out->table();
1657  Table t0 = tout(tout.col("POLNO") == 0);
1658  Table t1 = tout(tout.col("POLNO") == 1);
1659  if ( t0.nrow() != t1.nrow() )
1660    throw(AipsError("Inconsistent number of polarisations"));
1661  ArrayColumn<Float> speccol0(t0, "SPECTRA");
1662  ArrayColumn<uChar> flagcol0(t0, "FLAGTRA");
1663  ArrayColumn<Float> speccol1(t1, "SPECTRA");
1664  ArrayColumn<uChar> flagcol1(t1, "FLAGTRA");
1665  Matrix<Float> s0 = speccol0.getColumn();
1666  Matrix<uChar> f0 = flagcol0.getColumn();
1667  speccol0.putColumn(speccol1.getColumn());
1668  flagcol0.putColumn(flagcol1.getColumn());
1669  speccol1.putColumn(s0);
1670  flagcol1.putColumn(f0);
1671  return out;
1672}
1673
1674CountedPtr< Scantable >
1675  STMath::averagePolarisations( const CountedPtr< Scantable > & in,
1676                                const std::vector<bool>& mask,
1677                                const std::string& weight )
1678{
1679  if (in->npol() < 2 )
1680    throw(AipsError("averagePolarisations can only be applied to two or more"
1681                    "polarisations"));
1682  bool insitu = insitu_;
1683  setInsitu(false);
1684  CountedPtr< Scantable > pols = getScantable(in, true);
1685  setInsitu(insitu);
1686  Table& tout = pols->table();
1687  std::string taql = "SELECT FROM $1 WHERE POLNO IN [0,1]";
1688  Table tab = tableCommand(taql, in->table());
1689  if (tab.nrow() == 0 )
1690    throw(AipsError("Could not find  any rows with POLNO==0 and POLNO==1"));
1691  TableCopy::copyRows(tout, tab);
1692  TableVector<uInt> vec(tout, "POLNO");
1693  vec = 0;
1694  pols->table_.rwKeywordSet().define("nPol", Int(1));
1695  //pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
1696  pols->table_.rwKeywordSet().define("POLTYPE", in->getPolType());
1697  std::vector<CountedPtr<Scantable> > vpols;
1698  vpols.push_back(pols);
1699  CountedPtr< Scantable > out = average(vpols, mask, weight, "SCAN");
1700  return out;
1701}
1702
1703CountedPtr< Scantable >
1704  STMath::averageBeams( const CountedPtr< Scantable > & in,
1705                        const std::vector<bool>& mask,
1706                        const std::string& weight )
1707{
1708  bool insitu = insitu_;
1709  setInsitu(false);
1710  CountedPtr< Scantable > beams = getScantable(in, false);
1711  setInsitu(insitu);
1712  Table& tout = beams->table();
1713  // give all rows the same BEAMNO
1714  TableVector<uInt> vec(tout, "BEAMNO");
1715  vec = 0;
1716  beams->table_.rwKeywordSet().define("nBeam", Int(1));
1717  std::vector<CountedPtr<Scantable> > vbeams;
1718  vbeams.push_back(beams);
1719  CountedPtr< Scantable > out = average(vbeams, mask, weight, "SCAN");
1720  return out;
1721}
1722
1723
1724CountedPtr< Scantable >
1725  asap::STMath::frequencyAlign( const CountedPtr< Scantable > & in,
1726                                const std::string & refTime,
1727                                const std::string & method)
1728{
1729  // clone as this is not working insitu
1730  bool insitu = insitu_;
1731  setInsitu(false);
1732  CountedPtr< Scantable > out = getScantable(in, false);
1733  setInsitu(insitu);
1734  Table& tout = out->table();
1735  // Get reference Epoch to time of first row or given String
1736  Unit DAY(String("d"));
1737  MEpoch::Ref epochRef(in->getTimeReference());
1738  MEpoch refEpoch;
1739  if (refTime.length()>0) {
1740    Quantum<Double> qt;
1741    if (MVTime::read(qt,refTime)) {
1742      MVEpoch mv(qt);
1743      refEpoch = MEpoch(mv, epochRef);
1744   } else {
1745      throw(AipsError("Invalid format for Epoch string"));
1746   }
1747  } else {
1748    refEpoch = in->timeCol_(0);
1749  }
1750  MPosition refPos = in->getAntennaPosition();
1751
1752  InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
1753  /*
1754  // Comment from MV.
1755  // the following code has been commented out because different FREQ_IDs have to be aligned together even
1756  // if the frame doesn't change. So far, lack of this check didn't cause any problems.
1757  // test if user frame is different to base frame
1758  if ( in->frequencies().getFrameString(true)
1759       == in->frequencies().getFrameString(false) ) {
1760    throw(AipsError("Can't convert as no output frame has been set"
1761                    " (use set_freqframe) or it is aligned already."));
1762  }
1763  */
1764  MFrequency::Types system = in->frequencies().getFrame();
1765  MVTime mvt(refEpoch.getValue());
1766  String epochout = mvt.string(MVTime::YMD) + String(" (") + refEpoch.getRefString() + String(")");
1767  ostringstream oss;
1768  oss << "Aligned at reference Epoch " << epochout
1769      << " in frame " << MFrequency::showType(system);
1770  pushLog(String(oss));
1771  // set up the iterator
1772  Block<String> cols(4);
1773  // select by constant direction
1774  cols[0] = String("SRCNAME");
1775  cols[1] = String("BEAMNO");
1776  // select by IF ( no of channels varies over this )
1777  cols[2] = String("IFNO");
1778  // select by restfrequency
1779  cols[3] = String("MOLECULE_ID");
1780  TableIterator iter(tout, cols);
1781  while ( !iter.pastEnd() ) {
1782    Table t = iter.table();
1783    MDirection::ROScalarColumn dirCol(t, "DIRECTION");
1784    TableIterator fiter(t, "FREQ_ID");
1785    // determine nchan from the first row. This should work as
1786    // we are iterating over BEAMNO and IFNO    // we should have constant direction
1787
1788    ROArrayColumn<Float> sCol(t, "SPECTRA");
1789    const MDirection direction = dirCol(0);
1790    const uInt nchan = sCol(0).nelements();
1791
1792    // skip operations if there is nothing to align
1793    if (fiter.pastEnd()) {
1794        continue;
1795    }
1796
1797    Table ftab = fiter.table();
1798    // align all frequency ids with respect to the first encountered id
1799    ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
1800    // get the SpectralCoordinate for the freqid, which we are iterating over
1801    SpectralCoordinate sC = in->frequencies().getSpectralCoordinate(freqidCol(0));
1802    FrequencyAligner<Float> fa( sC, nchan, refEpoch,
1803                                direction, refPos, system );
1804    // realign the SpectralCoordinate and put into the output Scantable
1805    Vector<String> units(1);
1806    units = String("Hz");
1807    Bool linear=True;
1808    SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
1809    sc2.setWorldAxisUnits(units);
1810    const uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
1811                                                sc2.referenceValue()[0],
1812                                                sc2.increment()[0]);
1813    while ( !fiter.pastEnd() ) {
1814      ftab = fiter.table();
1815      // spectral coordinate for the current FREQ_ID
1816      ScalarColumn<uInt> freqidCol2(ftab, "FREQ_ID");
1817      sC = in->frequencies().getSpectralCoordinate(freqidCol2(0));
1818      // create the "global" abcissa for alignment with same FREQ_ID
1819      Vector<Double> abc(nchan);
1820      for (uInt i=0; i<nchan; i++) {
1821           Double w;
1822           sC.toWorld(w,Double(i));
1823           abc[i] = w;
1824      }
1825      TableVector<uInt> tvec(ftab, "FREQ_ID");
1826      // assign new frequency id to all rows
1827      tvec = id;
1828      // cache abcissa for same time stamps, so iterate over those
1829      TableIterator timeiter(ftab, "TIME");
1830      while ( !timeiter.pastEnd() ) {
1831        Table tab = timeiter.table();
1832        ArrayColumn<Float> specCol(tab, "SPECTRA");
1833        ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1834        MEpoch::ROScalarColumn timeCol(tab, "TIME");
1835        // use align abcissa cache after the first row
1836        // these rows should be just be POLNO
1837        bool first = true;
1838        for (int i=0; i<int(tab.nrow()); ++i) {
1839          // input values
1840          Vector<uChar> flag = flagCol(i);
1841          Vector<Bool> mask(flag.shape());
1842          Vector<Float> specOut, spec;
1843          spec  = specCol(i);
1844          Vector<Bool> maskOut;Vector<uChar> flagOut;
1845          convertArray(mask, flag);
1846          // alignment
1847          Bool ok = fa.align(specOut, maskOut, abc, spec,
1848                             mask, timeCol(i), !first,
1849                             interp, False);
1850          // back into scantable
1851          flagOut.resize(maskOut.nelements());
1852          convertArray(flagOut, maskOut);
1853          flagCol.put(i, flagOut);
1854          specCol.put(i, specOut);
1855          // start abcissa caching
1856          first = false;
1857        }
1858        // next timestamp
1859        ++timeiter;
1860      }
1861      // next FREQ_ID
1862      ++fiter;
1863    }
1864    // next aligner
1865    ++iter;
1866  }
1867  // set this afterwards to ensure we are doing insitu correctly.
1868  out->frequencies().setFrame(system, true);
1869  return out;
1870}
1871
1872CountedPtr<Scantable>
1873  asap::STMath::convertPolarisation( const CountedPtr<Scantable>& in,
1874                                     const std::string & newtype )
1875{
1876  if (in->npol() != 2 && in->npol() != 4)
1877    throw(AipsError("Can only convert two or four polarisations."));
1878  if ( in->getPolType() == newtype )
1879    throw(AipsError("No need to convert."));
1880  if ( ! in->selector_.empty() )
1881    throw(AipsError("Can only convert whole scantable. Unset the selection."));
1882  bool insitu = insitu_;
1883  setInsitu(false);
1884  CountedPtr< Scantable > out = getScantable(in, true);
1885  setInsitu(insitu);
1886  Table& tout = out->table();
1887  tout.rwKeywordSet().define("POLTYPE", String(newtype));
1888
1889  Block<String> cols(4);
1890  cols[0] = "SCANNO";
1891  cols[1] = "CYCLENO";
1892  cols[2] = "BEAMNO";
1893  cols[3] = "IFNO";
1894  TableIterator it(in->originalTable_, cols);
1895  String basetype = in->getPolType();
1896  STPol* stpol = STPol::getPolClass(in->factories_, basetype);
1897  try {
1898    while ( !it.pastEnd() ) {
1899      Table tab = it.table();
1900      uInt row = tab.rowNumbers()[0];
1901      stpol->setSpectra(in->getPolMatrix(row));
1902      Float fang,fhand,parang;
1903      fang = in->focusTable_.getTotalFeedAngle(in->mfocusidCol_(row));
1904      fhand = in->focusTable_.getFeedHand(in->mfocusidCol_(row));
1905      parang = in->paraCol_(row);
1906      /// @todo re-enable this
1907      // disable total feed angle to support paralactifying Caswell style
1908      stpol->setPhaseCorrections(parang, -parang, fhand);
1909      Int npolout = 0;
1910      for (uInt i=0; i<tab.nrow(); ++i) {
1911        Vector<Float> outvec = stpol->getSpectrum(i, newtype);
1912        if ( outvec.nelements() > 0 ) {
1913          tout.addRow();
1914          TableCopy::copyRows(tout, tab, tout.nrow()-1, 0, 1);
1915          ArrayColumn<Float> sCol(tout,"SPECTRA");
1916          ScalarColumn<uInt> pCol(tout,"POLNO");
1917          sCol.put(tout.nrow()-1 ,outvec);
1918          pCol.put(tout.nrow()-1 ,uInt(npolout));
1919          npolout++;
1920       }
1921      }
1922      tout.rwKeywordSet().define("nPol", npolout);
1923      ++it;
1924    }
1925  } catch (AipsError& e) {
1926    delete stpol;
1927    throw(e);
1928  }
1929  delete stpol;
1930  return out;
1931}
1932
1933CountedPtr< Scantable >
1934  asap::STMath::mxExtract( const CountedPtr< Scantable > & in,
1935                           const std::string & scantype )
1936{
1937  bool insitu = insitu_;
1938  setInsitu(false);
1939  CountedPtr< Scantable > out = getScantable(in, true);
1940  setInsitu(insitu);
1941  Table& tout = out->table();
1942  std::string taql = "SELECT FROM $1 WHERE BEAMNO != REFBEAMNO";
1943  if (scantype == "on") {
1944    taql = "SELECT FROM $1 WHERE BEAMNO == REFBEAMNO";
1945  }
1946  Table tab = tableCommand(taql, in->table());
1947  TableCopy::copyRows(tout, tab);
1948  if (scantype == "on") {
1949    // re-index SCANNO to 0
1950    TableVector<uInt> vec(tout, "SCANNO");
1951    vec = 0;
1952  }
1953  return out;
1954}
1955
1956CountedPtr< Scantable >
1957  asap::STMath::lagFlag( const CountedPtr< Scantable > & in,
1958                          double frequency, double width )
1959{
1960  CountedPtr< Scantable > out = getScantable(in, false);
1961  Table& tout = out->table();
1962  TableIterator iter(tout, "FREQ_ID");
1963  FFTServer<Float,Complex> ffts;
1964  while ( !iter.pastEnd() ) {
1965    Table tab = iter.table();
1966    Double rp,rv,inc;
1967    ROTableRow row(tab);
1968    const TableRecord& rec = row.get(0);
1969    uInt freqid = rec.asuInt("FREQ_ID");
1970    out->frequencies().getEntry(rp, rv, inc, freqid);
1971    ArrayColumn<Float> specCol(tab, "SPECTRA");
1972    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1973    for (int i=0; i<int(tab.nrow()); ++i) {
1974      Vector<Float> spec = specCol(i);
1975      Vector<uChar> flag = flagCol(i);
1976      Int lag0 = Int(spec.nelements()*abs(inc)/(frequency+width)+0.5);
1977      Int lag1 = Int(spec.nelements()*abs(inc)/(frequency-width)+0.5);
1978      for (int k=0; k < flag.nelements(); ++k ) {
1979        if (flag[k] > 0) {
1980          spec[k] = 0.0;
1981        }
1982      }
1983      Vector<Complex> lags;
1984      ffts.fft0(lags, spec);
1985      Int start =  max(0, lag0);
1986      Int end =  min(Int(lags.nelements()-1), lag1);
1987      if (start == end) {
1988        lags[start] = Complex(0.0);
1989      } else {
1990        for (int j=start; j <=end ;++j) {
1991          lags[j] = Complex(0.0);
1992        }
1993      }
1994      ffts.fft0(spec, lags);
1995      specCol.put(i, spec);
1996    }
1997    ++iter;
1998  }
1999  return out;
2000}
Note: See TracBrowser for help on using the repository browser.