source: trunk/src/STMath.cpp @ 1321

Last change on this file since 1321 was 1321, checked in by mar637, 17 years ago

Hopefully final fix to ticket #78;had to do a sub selection of BEAMNO,IFNO,POLNO fisrt, and also TIME column is MJD, which is in datys not seconds.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 49.1 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
51using namespace casa;
52
53using namespace asap;
54
55STMath::STMath(bool insitu) :
56  insitu_(insitu)
57{
58}
59
60
61STMath::~STMath()
62{
63}
64
65CountedPtr<Scantable>
66STMath::average( const std::vector<CountedPtr<Scantable> >& in,
67                 const std::vector<bool>& mask,
68                 const std::string& weight,
69                 const std::string& avmode)
70{
71  if ( avmode == "SCAN" && in.size() != 1 )
72    throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
73                    "Use merge first."));
74  WeightType wtype = stringToWeight(weight);
75
76  // output
77  // clone as this is non insitu
78  bool insitu = insitu_;
79  setInsitu(false);
80  CountedPtr< Scantable > out = getScantable(in[0], true);
81  setInsitu(insitu);
82  std::vector<CountedPtr<Scantable> >::const_iterator stit = in.begin();
83  ++stit;
84  while ( stit != in.end() ) {
85    out->appendToHistoryTable((*stit)->history());
86    ++stit;
87  }
88
89  Table& tout = out->table();
90
91  /// @todo check if all scantables are conformant
92
93  ArrayColumn<Float> specColOut(tout,"SPECTRA");
94  ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
95  ArrayColumn<Float> tsysColOut(tout,"TSYS");
96  ScalarColumn<Double> mjdColOut(tout,"TIME");
97  ScalarColumn<Double> intColOut(tout,"INTERVAL");
98  ScalarColumn<uInt> cycColOut(tout,"CYCLENO");
99  ScalarColumn<uInt> scanColOut(tout,"SCANNO");
100
101  // set up the output table rows. These are based on the structure of the
102  // FIRST scantable in the vector
103  const Table& baset = in[0]->table();
104
105  Block<String> cols(3);
106  cols[0] = String("BEAMNO");
107  cols[1] = String("IFNO");
108  cols[2] = String("POLNO");
109  if ( avmode == "SOURCE" ) {
110    cols.resize(4);
111    cols[3] = String("SRCNAME");
112  }
113  if ( avmode == "SCAN"  && in.size() == 1) {
114    cols.resize(4);
115    cols[3] = String("SCANNO");
116  }
117  uInt outrowCount = 0;
118  TableIterator iter(baset, cols);
119  while (!iter.pastEnd()) {
120    Table subt = iter.table();
121    // copy the first row of this selection into the new table
122    tout.addRow();
123    TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
124    // re-index to 0
125    if ( avmode != "SCAN" && avmode != "SOURCE" ) {
126      scanColOut.put(outrowCount, uInt(0));
127    }
128    ++outrowCount;
129    ++iter;
130  }
131  RowAccumulator acc(wtype);
132  Vector<Bool> cmask(mask);
133  acc.setUserMask(cmask);
134  ROTableRow row(tout);
135  ROArrayColumn<Float> specCol, tsysCol;
136  ROArrayColumn<uChar> flagCol;
137  ROScalarColumn<Double> mjdCol, intCol;
138  ROScalarColumn<Int> scanIDCol;
139
140  for (uInt i=0; i < tout.nrow(); ++i) {
141    for ( int j=0; j < int(in.size()); ++j ) {
142      const Table& tin = in[j]->table();
143      const TableRecord& rec = row.get(i);
144      ROScalarColumn<Double> tmp(tin, "TIME");
145      Double td;tmp.get(0,td);
146      Table basesubt = tin(tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
147                       && tin.col("IFNO") == Int(rec.asuInt("IFNO"))
148                       && tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
149      Table subt;
150      if ( avmode == "SOURCE") {
151        subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME") );
152      } else if (avmode == "SCAN") {
153        subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
154      } else {
155        subt = basesubt;
156      }
157      specCol.attach(subt,"SPECTRA");
158      flagCol.attach(subt,"FLAGTRA");
159      tsysCol.attach(subt,"TSYS");
160      intCol.attach(subt,"INTERVAL");
161      mjdCol.attach(subt,"TIME");
162      Vector<Float> spec,tsys;
163      Vector<uChar> flag;
164      Double inter,time;
165      for (uInt k = 0; k < subt.nrow(); ++k ) {
166        flagCol.get(k, flag);
167        Vector<Bool> bflag(flag.shape());
168        convertArray(bflag, flag);
169        /*
170        if ( allEQ(bflag, True) ) {
171        continue;//don't accumulate
172        }
173        */
174        specCol.get(k, spec);
175        tsysCol.get(k, tsys);
176        intCol.get(k, inter);
177        mjdCol.get(k, time);
178        // spectrum has to be added last to enable weighting by the other values
179        acc.add(spec, !bflag, tsys, inter, time);
180      }
181    }
182    //write out
183    specColOut.put(i, acc.getSpectrum());
184    const Vector<Bool>& msk = acc.getMask();
185    Vector<uChar> flg(msk.shape());
186    convertArray(flg, !msk);
187    flagColOut.put(i, flg);
188    tsysColOut.put(i, acc.getTsys());
189    intColOut.put(i, acc.getInterval());
190    mjdColOut.put(i, acc.getTime());
191    // we should only have one cycle now -> reset it to be 0
192    // frequency switched data has different CYCLENO for different IFNO
193    // which requires resetting this value
194    cycColOut.put(i, uInt(0));
195    acc.reset();
196  }
197  return out;
198}
199
200CountedPtr< Scantable >
201  STMath::averageChannel( const CountedPtr < Scantable > & in,
202                          const std::string & mode,
203                          const std::string& avmode )
204{
205  // clone as this is non insitu
206  bool insitu = insitu_;
207  setInsitu(false);
208  CountedPtr< Scantable > out = getScantable(in, true);
209  setInsitu(insitu);
210  Table& tout = out->table();
211  ArrayColumn<Float> specColOut(tout,"SPECTRA");
212  ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
213  ArrayColumn<Float> tsysColOut(tout,"TSYS");
214  ScalarColumn<uInt> scanColOut(tout,"SCANNO");
215  ScalarColumn<Double> intColOut(tout, "INTERVAL");
216  Table tmp = in->table().sort("BEAMNO");
217  Block<String> cols(3);
218  cols[0] = String("BEAMNO");
219  cols[1] = String("IFNO");
220  cols[2] = String("POLNO");
221  if ( avmode == "SCAN") {
222    cols.resize(4);
223    cols[3] = String("SCANNO");
224  }
225  uInt outrowCount = 0;
226  uChar userflag = 1 << 7;
227  TableIterator iter(tmp, cols);
228  while (!iter.pastEnd()) {
229    Table subt = iter.table();
230    ROArrayColumn<Float> specCol, tsysCol;
231    ROArrayColumn<uChar> flagCol;
232    ROScalarColumn<Double> intCol(subt, "INTERVAL");
233    specCol.attach(subt,"SPECTRA");
234    flagCol.attach(subt,"FLAGTRA");
235    tsysCol.attach(subt,"TSYS");
236    tout.addRow();
237    TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
238    if ( avmode != "SCAN") {
239      scanColOut.put(outrowCount, uInt(0));
240    }
241    Vector<Float> tmp;
242    specCol.get(0, tmp);
243    uInt nchan = tmp.nelements();
244    // have to do channel by channel here as MaskedArrMath
245    // doesn't have partialMedians
246    Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
247    Vector<Float> outspec(nchan);
248    Vector<uChar> outflag(nchan,0);
249    Vector<Float> outtsys(1);/// @fixme when tsys is channel based
250    for (uInt i=0; i<nchan; ++i) {
251      Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
252      MaskedArray<Float> ma = maskedArray(specs,flags);
253      outspec[i] = median(ma);
254      if ( allEQ(ma.getMask(), False) )
255        outflag[i] = userflag;// flag data
256    }
257    outtsys[0] = median(tsysCol.getColumn());
258    specColOut.put(outrowCount, outspec);
259    flagColOut.put(outrowCount, outflag);
260    tsysColOut.put(outrowCount, outtsys);
261    Double intsum = sum(intCol.getColumn());
262    intColOut.put(outrowCount, intsum);
263    ++outrowCount;
264    ++iter;
265  }
266  return out;
267}
268
269CountedPtr< Scantable > STMath::getScantable(const CountedPtr< Scantable >& in,
270                                             bool droprows)
271{
272  if (insitu_) return in;
273  else {
274    // clone
275    Scantable* tabp = new Scantable(*in, Bool(droprows));
276    return CountedPtr<Scantable>(tabp);
277  }
278}
279
280CountedPtr< Scantable > STMath::unaryOperate( const CountedPtr< Scantable >& in,
281                                              float val,
282                                              const std::string& mode,
283                                              bool tsys )
284{
285  CountedPtr< Scantable > out = getScantable(in, false);
286  Table& tab = out->table();
287  ArrayColumn<Float> specCol(tab,"SPECTRA");
288  ArrayColumn<Float> tsysCol(tab,"TSYS");
289  for (uInt i=0; i<tab.nrow(); ++i) {
290    Vector<Float> spec;
291    Vector<Float> ts;
292    specCol.get(i, spec);
293    tsysCol.get(i, ts);
294    if (mode == "MUL" || mode == "DIV") {
295      if (mode == "DIV") val = 1.0/val;
296      spec *= val;
297      specCol.put(i, spec);
298      if ( tsys ) {
299        ts *= val;
300        tsysCol.put(i, ts);
301      }
302    } else if ( mode == "ADD"  || mode == "SUB") {
303      if (mode == "SUB") val *= -1.0;
304      spec += val;
305      specCol.put(i, spec);
306      if ( tsys ) {
307        ts += val;
308        tsysCol.put(i, ts);
309      }
310    }
311  }
312  return out;
313}
314
315CountedPtr<Scantable> STMath::binaryOperate(const CountedPtr<Scantable>& left,
316                                            const CountedPtr<Scantable>& right,
317                                            const std::string& mode)
318{
319  bool insitu = insitu_;
320  if ( ! left->conformant(*right) ) {
321    throw(AipsError("'left' and 'right' scantables are not conformant."));
322  }
323  setInsitu(false);
324  CountedPtr< Scantable > out = getScantable(left, false);
325  setInsitu(insitu);
326  Table& tout = out->table();
327  Block<String> coln(5);
328  coln[0] = "SCANNO";  coln[1] = "CYCLENO";  coln[2] = "BEAMNO";
329  coln[3] = "IFNO";  coln[4] = "POLNO";
330  Table tmpl = tout.sort(coln);
331  Table tmpr = right->table().sort(coln);
332  ArrayColumn<Float> lspecCol(tmpl,"SPECTRA");
333  ROArrayColumn<Float> rspecCol(tmpr,"SPECTRA");
334  ArrayColumn<uChar> lflagCol(tmpl,"FLAGTRA");
335  ROArrayColumn<uChar> rflagCol(tmpr,"FLAGTRA");
336
337  for (uInt i=0; i<tout.nrow(); ++i) {
338    Vector<Float> lspecvec, rspecvec;
339    Vector<uChar> lflagvec, rflagvec;
340    lspecvec = lspecCol(i);    rspecvec = rspecCol(i);
341    lflagvec = lflagCol(i);    rflagvec = rflagCol(i);
342    MaskedArray<Float> mleft = maskedArray(lspecvec, lflagvec);
343    MaskedArray<Float> mright = maskedArray(rspecvec, rflagvec);
344    if (mode == "ADD") {
345      mleft += mright;
346    } else if ( mode == "SUB") {
347      mleft -= mright;
348    } else if ( mode == "MUL") {
349      mleft *= mright;
350    } else if ( mode == "DIV") {
351      mleft /= mright;
352    } else {
353      throw(AipsError("Illegal binary operator"));
354    }
355    lspecCol.put(i, mleft.getArray());
356  }
357  return out;
358}
359
360
361
362MaskedArray<Float> STMath::maskedArray( const Vector<Float>& s,
363                                        const Vector<uChar>& f)
364{
365  Vector<Bool> mask;
366  mask.resize(f.shape());
367  convertArray(mask, f);
368  return MaskedArray<Float>(s,!mask);
369}
370
371Vector<uChar> STMath::flagsFromMA(const MaskedArray<Float>& ma)
372{
373  const Vector<Bool>& m = ma.getMask();
374  Vector<uChar> flags(m.shape());
375  convertArray(flags, !m);
376  return flags;
377}
378
379CountedPtr< Scantable > STMath::autoQuotient( const CountedPtr< Scantable >& in,
380                                              const std::string & mode,
381                                              bool preserve )
382{
383  /// @todo make other modes available
384  /// modes should be "nearest", "pair"
385  // make this operation non insitu
386  const Table& tin = in->table();
387  Table ons = tin(tin.col("SRCTYPE") == Int(0));
388  Table offs = tin(tin.col("SRCTYPE") == Int(1));
389  if ( offs.nrow() == 0 )
390    throw(AipsError("No 'off' scans present."));
391  // put all "on" scans into output table
392
393  bool insitu = insitu_;
394  setInsitu(false);
395  CountedPtr< Scantable > out = getScantable(in, true);
396  setInsitu(insitu);
397  Table& tout = out->table();
398
399  TableCopy::copyRows(tout, ons);
400  TableRow row(tout);
401  ROScalarColumn<Double> offtimeCol(offs, "TIME");
402  ArrayColumn<Float> outspecCol(tout, "SPECTRA");
403  ROArrayColumn<Float> outtsysCol(tout, "TSYS");
404  ArrayColumn<uChar> outflagCol(tout, "FLAGTRA");
405  for (uInt i=0; i < tout.nrow(); ++i) {
406    const TableRecord& rec = row.get(i);
407    Double ontime = rec.asDouble("TIME");
408    Table presel = offs(offs.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
409                        && offs.col("IFNO") == Int(rec.asuInt("IFNO"))
410                        && offs.col("POLNO") == Int(rec.asuInt("POLNO")) );
411    ROScalarColumn<Double> offtimeCol(presel, "TIME");
412
413    Double mindeltat = min(abs(offtimeCol.getColumn() - ontime));
414    // Timestamp may vary within a cycle ???!!!
415    // increase this by 0.01 sec in case of rounding errors...
416    // There might be a better way to do this.
417    // fix to this fix. TIME is MJD, so 1.0d not 1.0s
418    mindeltat += 0.01/24./60./60.;
419    Table sel = presel( abs(presel.col("TIME")-ontime) <= mindeltat);
420
421    if ( sel.nrow() < 1 )  {
422      throw(AipsError("No closest in time found... This could be a rounding "
423                      "issue. Try quotient instead."));
424    }
425    TableRow offrow(sel);
426    const TableRecord& offrec = offrow.get(0);//should only be one row
427    RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
428    RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
429    RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
430    /// @fixme this assumes tsys is a scalar not vector
431    Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
432    Vector<Float> specon, tsyson;
433    outtsysCol.get(i, tsyson);
434    outspecCol.get(i, specon);
435    Vector<uChar> flagon;
436    outflagCol.get(i, flagon);
437    MaskedArray<Float> mon = maskedArray(specon, flagon);
438    MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
439    MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
440    if (preserve) {
441      quot -= tsysoffscalar;
442    } else {
443      quot -= tsyson[0];
444    }
445    outspecCol.put(i, quot.getArray());
446    outflagCol.put(i, flagsFromMA(quot));
447  }
448  // renumber scanno
449  TableIterator it(tout, "SCANNO");
450  uInt i = 0;
451  while ( !it.pastEnd() ) {
452    Table t = it.table();
453    TableVector<uInt> vec(t, "SCANNO");
454    vec = i;
455    ++i;
456    ++it;
457  }
458  return out;
459}
460
461
462CountedPtr< Scantable > STMath::quotient( const CountedPtr< Scantable > & on,
463                                          const CountedPtr< Scantable > & off,
464                                          bool preserve )
465{
466  bool insitu = insitu_;
467  if ( ! on->conformant(*off) ) {
468    throw(AipsError("'on' and 'off' scantables are not conformant."));
469  }
470  setInsitu(false);
471  CountedPtr< Scantable > out = getScantable(on, false);
472  setInsitu(insitu);
473  Table& tout = out->table();
474  const Table& toff = off->table();
475  TableIterator sit(tout, "SCANNO");
476  TableIterator s2it(toff, "SCANNO");
477  while ( !sit.pastEnd() ) {
478    Table ton = sit.table();
479    TableRow row(ton);
480    Table t = s2it.table();
481    ArrayColumn<Float> outspecCol(ton, "SPECTRA");
482    ROArrayColumn<Float> outtsysCol(ton, "TSYS");
483    ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
484    for (uInt i=0; i < ton.nrow(); ++i) {
485      const TableRecord& rec = row.get(i);
486      Table offsel = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
487                          && t.col("IFNO") == Int(rec.asuInt("IFNO"))
488                          && t.col("POLNO") == Int(rec.asuInt("POLNO")) );
489      if ( offsel.nrow() == 0 )
490        throw AipsError("STMath::quotient: no matching off");
491      TableRow offrow(offsel);
492      const TableRecord& offrec = offrow.get(0);//should be ncycles - take first
493      RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
494      RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
495      RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
496      Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
497      Vector<Float> specon, tsyson;
498      outtsysCol.get(i, tsyson);
499      outspecCol.get(i, specon);
500      Vector<uChar> flagon;
501      outflagCol.get(i, flagon);
502      MaskedArray<Float> mon = maskedArray(specon, flagon);
503      MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
504      MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
505      if (preserve) {
506        quot -= tsysoffscalar;
507      } else {
508        quot -= tsyson[0];
509      }
510      outspecCol.put(i, quot.getArray());
511      outflagCol.put(i, flagsFromMA(quot));
512    }
513    ++sit;
514    ++s2it;
515    // take the first off for each on scan which doesn't have a
516    // matching off scan
517    // non <= noff:  matching pairs, non > noff matching pairs then first off
518    if ( s2it.pastEnd() ) s2it.reset();
519  }
520  return out;
521}
522
523
524CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
525{
526  // make copy or reference
527  CountedPtr< Scantable > out = getScantable(in, false);
528  Table& tout = out->table();
529  Block<String> cols(4);
530  cols[0] = String("SCANNO");
531  cols[1] = String("CYCLENO");
532  cols[2] = String("BEAMNO");
533  cols[3] = String("POLNO");
534  TableIterator iter(tout, cols);
535  while (!iter.pastEnd()) {
536    Table subt = iter.table();
537    // this should leave us with two rows for the two IFs....if not ignore
538    if (subt.nrow() != 2 ) {
539      continue;
540    }
541    ArrayColumn<Float> specCol(subt, "SPECTRA");
542    ArrayColumn<Float> tsysCol(subt, "TSYS");
543    ArrayColumn<uChar> flagCol(subt, "FLAGTRA");
544    Vector<Float> onspec,offspec, ontsys, offtsys;
545    Vector<uChar> onflag, offflag;
546    tsysCol.get(0, ontsys);   tsysCol.get(1, offtsys);
547    specCol.get(0, onspec);   specCol.get(1, offspec);
548    flagCol.get(0, onflag);   flagCol.get(1, offflag);
549    MaskedArray<Float> on  = maskedArray(onspec, onflag);
550    MaskedArray<Float> off = maskedArray(offspec, offflag);
551    MaskedArray<Float> oncopy = on.copy();
552
553    on /= off; on -= 1.0f;
554    on *= ontsys[0];
555    off /= oncopy; off -= 1.0f;
556    off *= offtsys[0];
557    specCol.put(0, on.getArray());
558    const Vector<Bool>& m0 = on.getMask();
559    Vector<uChar> flags0(m0.shape());
560    convertArray(flags0, !m0);
561    flagCol.put(0, flags0);
562
563    specCol.put(1, off.getArray());
564    const Vector<Bool>& m1 = off.getMask();
565    Vector<uChar> flags1(m1.shape());
566    convertArray(flags1, !m1);
567    flagCol.put(1, flags1);
568    ++iter;
569  }
570
571  return out;
572}
573
574std::vector< float > STMath::statistic( const CountedPtr< Scantable > & in,
575                                        const std::vector< bool > & mask,
576                                        const std::string& which )
577{
578
579  Vector<Bool> m(mask);
580  const Table& tab = in->table();
581  ROArrayColumn<Float> specCol(tab, "SPECTRA");
582  ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
583  std::vector<float> out;
584  for (uInt i=0; i < tab.nrow(); ++i ) {
585    Vector<Float> spec; specCol.get(i, spec);
586    Vector<uChar> flag; flagCol.get(i, flag);
587    MaskedArray<Float> ma  = maskedArray(spec, flag);
588    float outstat = 0.0;
589    if ( spec.nelements() == m.nelements() ) {
590      outstat = mathutil::statistics(which, ma(m));
591    } else {
592      outstat = mathutil::statistics(which, ma);
593    }
594    out.push_back(outstat);
595  }
596  return out;
597}
598
599CountedPtr< Scantable > STMath::bin( const CountedPtr< Scantable > & in,
600                                     int width )
601{
602  if ( !in->getSelection().empty() ) throw(AipsError("Can't bin subset of the data."));
603  CountedPtr< Scantable > out = getScantable(in, false);
604  Table& tout = out->table();
605  out->frequencies().rescale(width, "BIN");
606  ArrayColumn<Float> specCol(tout, "SPECTRA");
607  ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
608  for (uInt i=0; i < tout.nrow(); ++i ) {
609    MaskedArray<Float> main  = maskedArray(specCol(i), flagCol(i));
610    MaskedArray<Float> maout;
611    LatticeUtilities::bin(maout, main, 0, Int(width));
612    /// @todo implement channel based tsys binning
613    specCol.put(i, maout.getArray());
614    flagCol.put(i, flagsFromMA(maout));
615    // take only the first binned spectrum's length for the deprecated
616    // global header item nChan
617    if (i==0) tout.rwKeywordSet().define(String("nChan"),
618                                       Int(maout.getArray().nelements()));
619  }
620  return out;
621}
622
623CountedPtr< Scantable > STMath::resample( const CountedPtr< Scantable >& in,
624                                          const std::string& method,
625                                          float width )
626//
627// Should add the possibility of width being specified in km/s. This means
628// that for each freqID (SpectralCoordinate) we will need to convert to an
629// average channel width (say at the reference pixel).  Then we would need
630// to be careful to make sure each spectrum (of different freqID)
631// is the same length.
632//
633{
634  //InterpolateArray1D<Double,Float>::InterpolationMethod interp;
635  Int interpMethod(stringToIMethod(method));
636
637  CountedPtr< Scantable > out = getScantable(in, false);
638  Table& tout = out->table();
639
640// Resample SpectralCoordinates (one per freqID)
641  out->frequencies().rescale(width, "RESAMPLE");
642  TableIterator iter(tout, "IFNO");
643  TableRow row(tout);
644  while ( !iter.pastEnd() ) {
645    Table tab = iter.table();
646    ArrayColumn<Float> specCol(tab, "SPECTRA");
647    //ArrayColumn<Float> tsysCol(tout, "TSYS");
648    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
649    Vector<Float> spec;
650    Vector<uChar> flag;
651    specCol.get(0,spec); // the number of channels should be constant per IF
652    uInt nChanIn = spec.nelements();
653    Vector<Float> xIn(nChanIn); indgen(xIn);
654    Int fac =  Int(nChanIn/width);
655    Vector<Float> xOut(fac+10); // 10 to be safe - resize later
656    uInt k = 0;
657    Float x = 0.0;
658    while (x < Float(nChanIn) ) {
659      xOut(k) = x;
660      k++;
661      x += width;
662    }
663    uInt nChanOut = k;
664    xOut.resize(nChanOut, True);
665    // process all rows for this IFNO
666    Vector<Float> specOut;
667    Vector<Bool> maskOut;
668    Vector<uChar> flagOut;
669    for (uInt i=0; i < tab.nrow(); ++i) {
670      specCol.get(i, spec);
671      flagCol.get(i, flag);
672      Vector<Bool> mask(flag.nelements());
673      convertArray(mask, flag);
674
675      IPosition shapeIn(spec.shape());
676      //sh.nchan = nChanOut;
677      InterpolateArray1D<Float,Float>::interpolate(specOut, maskOut, xOut,
678                                                   xIn, spec, mask,
679                                                   interpMethod, True, True);
680      /// @todo do the same for channel based Tsys
681      flagOut.resize(maskOut.nelements());
682      convertArray(flagOut, maskOut);
683      specCol.put(i, specOut);
684      flagCol.put(i, flagOut);
685    }
686    ++iter;
687  }
688
689  return out;
690}
691
692STMath::imethod STMath::stringToIMethod(const std::string& in)
693{
694  static STMath::imap lookup;
695
696  // initialize the lookup table if necessary
697  if ( lookup.empty() ) {
698    lookup["nearest"]   = InterpolateArray1D<Double,Float>::nearestNeighbour;
699    lookup["linear"] = InterpolateArray1D<Double,Float>::linear;
700    lookup["cubic"]  = InterpolateArray1D<Double,Float>::cubic;
701    lookup["spline"]  = InterpolateArray1D<Double,Float>::spline;
702  }
703
704  STMath::imap::const_iterator iter = lookup.find(in);
705
706  if ( lookup.end() == iter ) {
707    std::string message = in;
708    message += " is not a valid interpolation mode";
709    throw(AipsError(message));
710  }
711  return iter->second;
712}
713
714WeightType STMath::stringToWeight(const std::string& in)
715{
716  static std::map<std::string, WeightType> lookup;
717
718  // initialize the lookup table if necessary
719  if ( lookup.empty() ) {
720    lookup["NONE"]   = asap::NONE;
721    lookup["TINT"] = asap::TINT;
722    lookup["TINTSYS"]  = asap::TINTSYS;
723    lookup["TSYS"]  = asap::TSYS;
724    lookup["VAR"]  = asap::VAR;
725  }
726
727  std::map<std::string, WeightType>::const_iterator iter = lookup.find(in);
728
729  if ( lookup.end() == iter ) {
730    std::string message = in;
731    message += " is not a valid weighting mode";
732    throw(AipsError(message));
733  }
734  return iter->second;
735}
736
737CountedPtr< Scantable > STMath::gainElevation( const CountedPtr< Scantable >& in,
738                                               const vector< float > & coeff,
739                                               const std::string & filename,
740                                               const std::string& method)
741{
742  // Get elevation data from Scantable and convert to degrees
743  CountedPtr< Scantable > out = getScantable(in, false);
744  Table& tab = out->table();
745  ROScalarColumn<Float> elev(tab, "ELEVATION");
746  Vector<Float> x = elev.getColumn();
747  x *= Float(180 / C::pi);                        // Degrees
748
749  Vector<Float> coeffs(coeff);
750  const uInt nc = coeffs.nelements();
751  if ( filename.length() > 0 && nc > 0 ) {
752    throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
753  }
754
755  // Correct
756  if ( nc > 0 || filename.length() == 0 ) {
757    // Find instrument
758    Bool throwit = True;
759    Instrument inst =
760      STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
761                                throwit);
762
763    // Set polynomial
764    Polynomial<Float>* ppoly = 0;
765    Vector<Float> coeff;
766    String msg;
767    if ( nc > 0 ) {
768      ppoly = new Polynomial<Float>(nc);
769      coeff = coeffs;
770      msg = String("user");
771    } else {
772      STAttr sdAttr;
773      coeff = sdAttr.gainElevationPoly(inst);
774      ppoly = new Polynomial<Float>(3);
775      msg = String("built in");
776    }
777
778    if ( coeff.nelements() > 0 ) {
779      ppoly->setCoefficients(coeff);
780    } else {
781      delete ppoly;
782      throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
783    }
784    ostringstream oss;
785    oss << "Making polynomial correction with " << msg << " coefficients:" << endl;
786    oss << "   " <<  coeff;
787    pushLog(String(oss));
788    const uInt nrow = tab.nrow();
789    Vector<Float> factor(nrow);
790    for ( uInt i=0; i < nrow; ++i ) {
791      factor[i] = 1.0 / (*ppoly)(x[i]);
792    }
793    delete ppoly;
794    scaleByVector(tab, factor, true);
795
796  } else {
797    // Read and correct
798    pushLog("Making correction from ascii Table");
799    scaleFromAsciiTable(tab, filename, method, x, true);
800  }
801  return out;
802}
803
804void STMath::scaleFromAsciiTable(Table& in, const std::string& filename,
805                                 const std::string& method,
806                                 const Vector<Float>& xout, bool dotsys)
807{
808
809// Read gain-elevation ascii file data into a Table.
810
811  String formatString;
812  Table tbl = readAsciiTable(formatString, Table::Memory, filename, "", "", False);
813  scaleFromTable(in, tbl, method, xout, dotsys);
814}
815
816void STMath::scaleFromTable(Table& in,
817                            const Table& table,
818                            const std::string& method,
819                            const Vector<Float>& xout, bool dotsys)
820{
821
822  ROScalarColumn<Float> geElCol(table, "ELEVATION");
823  ROScalarColumn<Float> geFacCol(table, "FACTOR");
824  Vector<Float> xin = geElCol.getColumn();
825  Vector<Float> yin = geFacCol.getColumn();
826  Vector<Bool> maskin(xin.nelements(),True);
827
828  // Interpolate (and extrapolate) with desired method
829
830  InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
831
832   Vector<Float> yout;
833   Vector<Bool> maskout;
834   InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout,
835                                                xin, yin, maskin, interp,
836                                                True, True);
837
838   scaleByVector(in, Float(1.0)/yout, dotsys);
839}
840
841void STMath::scaleByVector( Table& in,
842                            const Vector< Float >& factor,
843                            bool dotsys )
844{
845  uInt nrow = in.nrow();
846  if ( factor.nelements() != nrow ) {
847    throw(AipsError("factors.nelements() != table.nelements()"));
848  }
849  ArrayColumn<Float> specCol(in, "SPECTRA");
850  ArrayColumn<uChar> flagCol(in, "FLAGTRA");
851  ArrayColumn<Float> tsysCol(in, "TSYS");
852  for (uInt i=0; i < nrow; ++i) {
853    MaskedArray<Float> ma  = maskedArray(specCol(i), flagCol(i));
854    ma *= factor[i];
855    specCol.put(i, ma.getArray());
856    flagCol.put(i, flagsFromMA(ma));
857    if ( dotsys ) {
858      Vector<Float> tsys = tsysCol(i);
859      tsys *= factor[i];
860      tsysCol.put(i,tsys);
861    }
862  }
863}
864
865CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in,
866                                             float d, float etaap,
867                                             float jyperk )
868{
869  CountedPtr< Scantable > out = getScantable(in, false);
870  Table& tab = in->table();
871  Unit fluxUnit(tab.keywordSet().asString("FluxUnit"));
872  Unit K(String("K"));
873  Unit JY(String("Jy"));
874
875  bool tokelvin = true;
876  Double cfac = 1.0;
877
878  if ( fluxUnit == JY ) {
879    pushLog("Converting to K");
880    Quantum<Double> t(1.0,fluxUnit);
881    Quantum<Double> t2 = t.get(JY);
882    cfac = (t2 / t).getValue();               // value to Jy
883
884    tokelvin = true;
885    out->setFluxUnit("K");
886  } else if ( fluxUnit == K ) {
887    pushLog("Converting to Jy");
888    Quantum<Double> t(1.0,fluxUnit);
889    Quantum<Double> t2 = t.get(K);
890    cfac = (t2 / t).getValue();              // value to K
891
892    tokelvin = false;
893    out->setFluxUnit("Jy");
894  } else {
895    throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
896  }
897  // Make sure input values are converted to either Jy or K first...
898  Float factor = cfac;
899
900  // Select method
901  if (jyperk > 0.0) {
902    factor *= jyperk;
903    if ( tokelvin ) factor = 1.0 / jyperk;
904    ostringstream oss;
905    oss << "Jy/K = " << jyperk;
906    pushLog(String(oss));
907    Vector<Float> factors(tab.nrow(), factor);
908    scaleByVector(tab,factors, false);
909  } else if ( etaap > 0.0) {
910    if (d < 0) {
911      Instrument inst =
912        STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
913                                  True);
914      STAttr sda;
915      d = sda.diameter(inst);
916    }
917    jyperk = STAttr::findJyPerK(etaap, d);
918    ostringstream oss;
919    oss << "Jy/K = " << jyperk;
920    pushLog(String(oss));
921    factor *= jyperk;
922    if ( tokelvin ) {
923      factor = 1.0 / factor;
924    }
925    Vector<Float> factors(tab.nrow(), factor);
926    scaleByVector(tab, factors, False);
927  } else {
928
929    // OK now we must deal with automatic look up of values.
930    // We must also deal with the fact that the factors need
931    // to be computed per IF and may be different and may
932    // change per integration.
933
934    pushLog("Looking up conversion factors");
935    convertBrightnessUnits(out, tokelvin, cfac);
936  }
937
938  return out;
939}
940
941void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in,
942                                     bool tokelvin, float cfac )
943{
944  Table& table = in->table();
945  Instrument inst =
946    STAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True);
947  TableIterator iter(table, "FREQ_ID");
948  STFrequencies stfreqs = in->frequencies();
949  STAttr sdAtt;
950  while (!iter.pastEnd()) {
951    Table tab = iter.table();
952    ArrayColumn<Float> specCol(tab, "SPECTRA");
953    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
954    ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID");
955    MEpoch::ROScalarColumn timeCol(tab, "TIME");
956
957    uInt freqid; freqidCol.get(0, freqid);
958    Vector<Float> tmpspec; specCol.get(0, tmpspec);
959    // STAttr.JyPerK has a Vector interface... change sometime.
960    Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements()));
961    for ( uInt i=0; i<tab.nrow(); ++i) {
962      Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0];
963      Float factor = cfac * jyperk;
964      if ( tokelvin ) factor = Float(1.0) / factor;
965      MaskedArray<Float> ma  = maskedArray(specCol(i), flagCol(i));
966      ma *= factor;
967      specCol.put(i, ma.getArray());
968      flagCol.put(i, flagsFromMA(ma));
969    }
970  ++iter;
971  }
972}
973
974CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
975                                         float tau )
976{
977  CountedPtr< Scantable > out = getScantable(in, false);
978
979  Table tab = out->table();
980  ROScalarColumn<Float> elev(tab, "ELEVATION");
981  ArrayColumn<Float> specCol(tab, "SPECTRA");
982  ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
983  for ( uInt i=0; i<tab.nrow(); ++i) {
984    Float zdist = Float(C::pi_2) - elev(i);
985    Float factor = exp(tau/cos(zdist));
986    MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
987    ma *= factor;
988    specCol.put(i, ma.getArray());
989    flagCol.put(i, flagsFromMA(ma));
990  }
991  return out;
992}
993
994CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
995                                        const std::string& kernel, float width )
996{
997  CountedPtr< Scantable > out = getScantable(in, false);
998  Table& table = out->table();
999  VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
1000  // same IFNO should have same no of channels
1001  // this saves overhead
1002  TableIterator iter(table, "IFNO");
1003  while (!iter.pastEnd()) {
1004    Table tab = iter.table();
1005    ArrayColumn<Float> specCol(tab, "SPECTRA");
1006    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1007    Vector<Float> tmpspec; specCol.get(0, tmpspec);
1008    uInt nchan = tmpspec.nelements();
1009    Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
1010    Convolver<Float> conv(kvec, IPosition(1,nchan));
1011    Vector<Float> spec;
1012    Vector<uChar> flag;
1013    for ( uInt i=0; i<tab.nrow(); ++i) {
1014      specCol.get(i, spec);
1015      flagCol.get(i, flag);
1016      Vector<Bool> mask(flag.nelements());
1017      convertArray(mask, flag);
1018      Vector<Float> specout;
1019      if ( type == VectorKernel::HANNING ) {
1020        Vector<Bool> maskout;
1021        mathutil::hanning(specout, maskout, spec , !mask);
1022        convertArray(flag, !maskout);
1023        flagCol.put(i, flag);
1024        specCol.put(i, specout);
1025     } else {
1026        mathutil::replaceMaskByZero(specout, mask);
1027        conv.linearConv(specout, spec);
1028        specCol.put(i, specout);
1029      }
1030    }
1031    ++iter;
1032  }
1033  return out;
1034}
1035
1036CountedPtr< Scantable >
1037  STMath::merge( const std::vector< CountedPtr < Scantable > >& in )
1038{
1039  if ( in.size() < 2 ) {
1040    throw(AipsError("Need at least two scantables to perform a merge."));
1041  }
1042  std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
1043  bool insitu = insitu_;
1044  setInsitu(false);
1045  CountedPtr< Scantable > out = getScantable(*it, false);
1046  setInsitu(insitu);
1047  Table& tout = out->table();
1048  ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
1049  ScalarColumn<uInt> scannocol(tout,"SCANNO"), focusidcol(tout,"FOCUS_ID");
1050  // Renumber SCANNO to be 0-based
1051  Vector<uInt> scannos = scannocol.getColumn();
1052  uInt offset = min(scannos);
1053  scannos -= offset;
1054  scannocol.putColumn(scannos);
1055  uInt newscanno = max(scannos)+1;
1056  ++it;
1057  while ( it != in.end() ){
1058    if ( ! (*it)->conformant(*out) ) {
1059      // log message: "ignoring scantable i, as it isn't
1060      // conformant with the other(s)"
1061      cerr << "oh oh" << endl;
1062      ++it;
1063      continue;
1064    }
1065    out->appendToHistoryTable((*it)->history());
1066    const Table& tab = (*it)->table();
1067    TableIterator scanit(tab, "SCANNO");
1068    while (!scanit.pastEnd()) {
1069      TableIterator freqit(scanit.table(), "FREQ_ID");
1070      while ( !freqit.pastEnd() ) {
1071        Table thetab = freqit.table();
1072        uInt nrow = tout.nrow();
1073        //tout.addRow(thetab.nrow());
1074        TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
1075        ROTableRow row(thetab);
1076        for ( uInt i=0; i<thetab.nrow(); ++i) {
1077          uInt k = nrow+i;
1078          scannocol.put(k, newscanno);
1079          const TableRecord& rec = row.get(i);
1080          Double rv,rp,inc;
1081          (*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
1082          uInt id;
1083          id = out->frequencies().addEntry(rp, rv, inc);
1084          freqidcol.put(k,id);
1085          String name,fname;Double rf;
1086          (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
1087          id = out->molecules().addEntry(rf, name, fname);
1088          molidcol.put(k, id);
1089          Float frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
1090          (*it)->focus().getEntry(fax, ftan, frot, fhand,
1091                                  fmount,fuser, fxy, fxyp,
1092                                  rec.asuInt("FOCUS_ID"));
1093          id = out->focus().addEntry(fax, ftan, frot, fhand,
1094                                     fmount,fuser, fxy, fxyp);
1095          focusidcol.put(k, id);
1096        }
1097        ++freqit;
1098      }
1099      ++newscanno;
1100      ++scanit;
1101    }
1102    ++it;
1103  }
1104  return out;
1105}
1106
1107CountedPtr< Scantable >
1108  STMath::invertPhase( const CountedPtr < Scantable >& in )
1109{
1110  return applyToPol(in, &STPol::invertPhase, Float(0.0));
1111}
1112
1113CountedPtr< Scantable >
1114  STMath::rotateXYPhase( const CountedPtr < Scantable >& in, float phase )
1115{
1116   return applyToPol(in, &STPol::rotatePhase, Float(phase));
1117}
1118
1119CountedPtr< Scantable >
1120  STMath::rotateLinPolPhase( const CountedPtr < Scantable >& in, float phase )
1121{
1122  return applyToPol(in, &STPol::rotateLinPolPhase, Float(phase));
1123}
1124
1125CountedPtr< Scantable > STMath::applyToPol( const CountedPtr<Scantable>& in,
1126                                             STPol::polOperation fptr,
1127                                             Float phase )
1128{
1129  CountedPtr< Scantable > out = getScantable(in, false);
1130  Table& tout = out->table();
1131  Block<String> cols(4);
1132  cols[0] = String("SCANNO");
1133  cols[1] = String("BEAMNO");
1134  cols[2] = String("IFNO");
1135  cols[3] = String("CYCLENO");
1136  TableIterator iter(tout, cols);
1137  STPol* stpol = STPol::getPolClass(out->factories_, out->getPolType() );
1138  while (!iter.pastEnd()) {
1139    Table t = iter.table();
1140    ArrayColumn<Float> speccol(t, "SPECTRA");
1141    ScalarColumn<uInt> focidcol(t, "FOCUS_ID");
1142    ScalarColumn<Float> parancol(t, "PARANGLE");
1143    Matrix<Float> pols = speccol.getColumn();
1144    try {
1145      stpol->setSpectra(pols);
1146      Float fang,fhand,parang;
1147      fang = in->focusTable_.getTotalFeedAngle(focidcol(0));
1148      fhand = in->focusTable_.getFeedHand(focidcol(0));
1149      parang = parancol(0);
1150      /// @todo re-enable this
1151      // disable total feed angle to support paralactifying Caswell style
1152      stpol->setPhaseCorrections(parang, -parang, fhand);
1153      (stpol->*fptr)(phase);
1154      speccol.putColumn(stpol->getSpectra());
1155      Matrix<Float> tmp = stpol->getSpectra();
1156    } catch (AipsError& e) {
1157      delete stpol;stpol=0;
1158      throw(e);
1159    }
1160    ++iter;
1161  }
1162  delete stpol;stpol=0;
1163  return out;
1164}
1165
1166CountedPtr< Scantable >
1167  STMath::swapPolarisations( const CountedPtr< Scantable > & in )
1168{
1169  CountedPtr< Scantable > out = getScantable(in, false);
1170  Table& tout = out->table();
1171  Table t0 = tout(tout.col("POLNO") == 0);
1172  Table t1 = tout(tout.col("POLNO") == 1);
1173  if ( t0.nrow() != t1.nrow() )
1174    throw(AipsError("Inconsistent number of polarisations"));
1175  ArrayColumn<Float> speccol0(t0, "SPECTRA");
1176  ArrayColumn<uChar> flagcol0(t0, "FLAGTRA");
1177  ArrayColumn<Float> speccol1(t1, "SPECTRA");
1178  ArrayColumn<uChar> flagcol1(t1, "FLAGTRA");
1179  Matrix<Float> s0 = speccol0.getColumn();
1180  Matrix<uChar> f0 = flagcol0.getColumn();
1181  speccol0.putColumn(speccol1.getColumn());
1182  flagcol0.putColumn(flagcol1.getColumn());
1183  speccol1.putColumn(s0);
1184  flagcol1.putColumn(f0);
1185  return out;
1186}
1187
1188CountedPtr< Scantable >
1189  STMath::averagePolarisations( const CountedPtr< Scantable > & in,
1190                                const std::vector<bool>& mask,
1191                                const std::string& weight )
1192{
1193  if (in->npol() < 2 )
1194    throw(AipsError("averagePolarisations can only be applied to two or more"
1195                    "polarisations"));
1196  bool insitu = insitu_;
1197  setInsitu(false);
1198  CountedPtr< Scantable > pols = getScantable(in, true);
1199  setInsitu(insitu);
1200  Table& tout = pols->table();
1201  std::string taql = "SELECT FROM $1 WHERE POLNO IN [0,1]";
1202  Table tab = tableCommand(taql, in->table());
1203  if (tab.nrow() == 0 )
1204    throw(AipsError("Could not find  any rows with POLNO==0 and POLNO==1"));
1205  TableCopy::copyRows(tout, tab);
1206  TableVector<uInt> vec(tout, "POLNO");
1207  vec = 0;
1208  pols->table_.rwKeywordSet().define("nPol", Int(1));
1209  pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
1210  std::vector<CountedPtr<Scantable> > vpols;
1211  vpols.push_back(pols);
1212  CountedPtr< Scantable > out = average(vpols, mask, weight, "SCAN");
1213  return out;
1214}
1215
1216CountedPtr< Scantable >
1217  STMath::averageBeams( const CountedPtr< Scantable > & in,
1218                        const std::vector<bool>& mask,
1219                        const std::string& weight )
1220{
1221  bool insitu = insitu_;
1222  setInsitu(false);
1223  CountedPtr< Scantable > beams = getScantable(in, false);
1224  setInsitu(insitu);
1225  Table& tout = beams->table();
1226  // give all rows the same BEAMNO
1227  TableVector<uInt> vec(tout, "BEAMNO");
1228  vec = 0;
1229  beams->table_.rwKeywordSet().define("nBeam", Int(1));
1230  std::vector<CountedPtr<Scantable> > vbeams;
1231  vbeams.push_back(beams);
1232  CountedPtr< Scantable > out = average(vbeams, mask, weight, "SCAN");
1233  return out;
1234}
1235
1236
1237CountedPtr< Scantable >
1238  asap::STMath::frequencyAlign( const CountedPtr< Scantable > & in,
1239                                const std::string & refTime,
1240                                const std::string & method)
1241{
1242  // clone as this is not working insitu
1243  bool insitu = insitu_;
1244  setInsitu(false);
1245  CountedPtr< Scantable > out = getScantable(in, false);
1246  setInsitu(insitu);
1247  Table& tout = out->table();
1248  // Get reference Epoch to time of first row or given String
1249  Unit DAY(String("d"));
1250  MEpoch::Ref epochRef(in->getTimeReference());
1251  MEpoch refEpoch;
1252  if (refTime.length()>0) {
1253    Quantum<Double> qt;
1254    if (MVTime::read(qt,refTime)) {
1255      MVEpoch mv(qt);
1256      refEpoch = MEpoch(mv, epochRef);
1257   } else {
1258      throw(AipsError("Invalid format for Epoch string"));
1259   }
1260  } else {
1261    refEpoch = in->timeCol_(0);
1262  }
1263  MPosition refPos = in->getAntennaPosition();
1264
1265  InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
1266  // test if user frame is different to base frame
1267  if ( in->frequencies().getFrameString(true)
1268       == in->frequencies().getFrameString(false) ) {
1269    throw(AipsError("Can't convert as no output frame has been set"
1270                    " (use set_freqframe) or it is aligned already."));
1271  }
1272  MFrequency::Types system = in->frequencies().getFrame();
1273  MVTime mvt(refEpoch.getValue());
1274  String epochout = mvt.string(MVTime::YMD) + String(" (") + refEpoch.getRefString() + String(")");
1275  ostringstream oss;
1276  oss << "Aligned at reference Epoch " << epochout
1277      << " in frame " << MFrequency::showType(system);
1278  pushLog(String(oss));
1279  // set up the iterator
1280  Block<String> cols(4);
1281  // select by constant direction
1282  cols[0] = String("SRCNAME");
1283  cols[1] = String("BEAMNO");
1284  // select by IF ( no of channels varies over this )
1285  cols[2] = String("IFNO");
1286  // select by restfrequency
1287  cols[3] = String("MOLECULE_ID");
1288  TableIterator iter(tout, cols);
1289  while ( !iter.pastEnd() ) {
1290    Table t = iter.table();
1291    MDirection::ROScalarColumn dirCol(t, "DIRECTION");
1292    TableIterator fiter(t, "FREQ_ID");
1293    // determine nchan from the first row. This should work as
1294    // we are iterating over BEAMNO and IFNO    // we should have constant direction
1295
1296    ROArrayColumn<Float> sCol(t, "SPECTRA");
1297    MDirection direction = dirCol(0);
1298    uInt nchan = sCol(0).nelements();
1299    while ( !fiter.pastEnd() ) {
1300      Table ftab = fiter.table();
1301      ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
1302      // get the SpectralCoordinate for the freqid, which we are iterating over
1303      SpectralCoordinate sC = in->frequencies().getSpectralCoordinate(freqidCol(0));
1304      FrequencyAligner<Float> fa( sC, nchan, refEpoch,
1305                                  direction, refPos, system );
1306      // realign the SpectralCoordinate and put into the output Scantable
1307      Vector<String> units(1);
1308      units = String("Hz");
1309      Bool linear=True;
1310      SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
1311      sc2.setWorldAxisUnits(units);
1312      uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
1313                                            sc2.referenceValue()[0],
1314                                            sc2.increment()[0]);
1315      TableVector<uInt> tvec(ftab, "FREQ_ID");
1316      tvec = id;
1317      // create the "global" abcissa for alignment with same FREQ_ID
1318      Vector<Double> abc(nchan);
1319      Double w;
1320      for (uInt i=0; i<nchan; i++) {
1321        sC.toWorld(w,Double(i));
1322        abc[i] = w;
1323      }
1324      // cache abcissa for same time stamps, so iterate over those
1325      TableIterator timeiter(ftab, "TIME");
1326      while ( !timeiter.pastEnd() ) {
1327        Table tab = timeiter.table();
1328        ArrayColumn<Float> specCol(tab, "SPECTRA");
1329        ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1330        MEpoch::ROScalarColumn timeCol(tab, "TIME");
1331        // use align abcissa cache after the first row
1332        bool first = true;
1333        // these rows should be just be POLNO
1334        for (int i=0; i<int(tab.nrow()); ++i) {
1335          // input values
1336          Vector<uChar> flag = flagCol(i);
1337          Vector<Bool> mask(flag.shape());
1338          Vector<Float> specOut, spec;
1339          spec  = specCol(i);
1340          Vector<Bool> maskOut;Vector<uChar> flagOut;
1341          convertArray(mask, flag);
1342          // alignment
1343          Bool ok = fa.align(specOut, maskOut, abc, spec,
1344                             mask, timeCol(i), !first,
1345                             interp, False);
1346          // back into scantable
1347          flagOut.resize(maskOut.nelements());
1348          convertArray(flagOut, maskOut);
1349          flagCol.put(i, flagOut);
1350          specCol.put(i, specOut);
1351          // start abcissa caching
1352          first = false;
1353        }
1354        // next timestamp
1355        ++timeiter;
1356      }
1357      // next FREQ_ID
1358      ++fiter;
1359    }
1360    // next aligner
1361    ++iter;
1362  }
1363  // set this afterwards to ensure we are doing insitu correctly.
1364  out->frequencies().setFrame(system, true);
1365  return out;
1366}
1367
1368CountedPtr<Scantable>
1369  asap::STMath::convertPolarisation( const CountedPtr<Scantable>& in,
1370                                     const std::string & newtype )
1371{
1372  if (in->npol() != 2 && in->npol() != 4)
1373    throw(AipsError("Can only convert two or four polarisations."));
1374  if ( in->getPolType() == newtype )
1375    throw(AipsError("No need to convert."));
1376  if ( ! in->selector_.empty() )
1377    throw(AipsError("Can only convert whole scantable. Unset the selection."));
1378  bool insitu = insitu_;
1379  setInsitu(false);
1380  CountedPtr< Scantable > out = getScantable(in, true);
1381  setInsitu(insitu);
1382  Table& tout = out->table();
1383  tout.rwKeywordSet().define("POLTYPE", String(newtype));
1384
1385  Block<String> cols(4);
1386  cols[0] = "SCANNO";
1387  cols[1] = "CYCLENO";
1388  cols[2] = "BEAMNO";
1389  cols[3] = "IFNO";
1390  TableIterator it(in->originalTable_, cols);
1391  String basetype = in->getPolType();
1392  STPol* stpol = STPol::getPolClass(in->factories_, basetype);
1393  try {
1394    while ( !it.pastEnd() ) {
1395      Table tab = it.table();
1396      uInt row = tab.rowNumbers()[0];
1397      stpol->setSpectra(in->getPolMatrix(row));
1398      Float fang,fhand,parang;
1399      fang = in->focusTable_.getTotalFeedAngle(in->mfocusidCol_(row));
1400      fhand = in->focusTable_.getFeedHand(in->mfocusidCol_(row));
1401      parang = in->paraCol_(row);
1402      /// @todo re-enable this
1403      // disable total feed angle to support paralactifying Caswell style
1404      stpol->setPhaseCorrections(parang, -parang, fhand);
1405      Int npolout = 0;
1406      for (uInt i=0; i<tab.nrow(); ++i) {
1407        Vector<Float> outvec = stpol->getSpectrum(i, newtype);
1408        if ( outvec.nelements() > 0 ) {
1409          tout.addRow();
1410          TableCopy::copyRows(tout, tab, tout.nrow()-1, 0, 1);
1411          ArrayColumn<Float> sCol(tout,"SPECTRA");
1412          ScalarColumn<uInt> pCol(tout,"POLNO");
1413          sCol.put(tout.nrow()-1 ,outvec);
1414          pCol.put(tout.nrow()-1 ,uInt(npolout));
1415          npolout++;
1416       }
1417      }
1418      tout.rwKeywordSet().define("nPol", npolout);
1419      ++it;
1420    }
1421  } catch (AipsError& e) {
1422    delete stpol;
1423    throw(e);
1424  }
1425  delete stpol;
1426  return out;
1427}
1428
1429CountedPtr< Scantable >
1430  asap::STMath::mxExtract( const CountedPtr< Scantable > & in,
1431                           const std::string & scantype )
1432{
1433  bool insitu = insitu_;
1434  setInsitu(false);
1435  CountedPtr< Scantable > out = getScantable(in, true);
1436  setInsitu(insitu);
1437  Table& tout = out->table();
1438  std::string taql = "SELECT FROM $1 WHERE BEAMNO != REFBEAMNO";
1439  if (scantype == "on") {
1440    taql = "SELECT FROM $1 WHERE BEAMNO == REFBEAMNO";
1441  }
1442  Table tab = tableCommand(taql, in->table());
1443  TableCopy::copyRows(tout, tab);
1444  if (scantype == "on") {
1445    // re-index SCANNO to 0
1446    TableVector<uInt> vec(tout, "SCANNO");
1447    vec = 0;
1448  }
1449  return out;
1450}
1451
1452CountedPtr< Scantable >
1453  asap::STMath::lagFlag( const CountedPtr< Scantable > & in,
1454                          double frequency, double width )
1455{
1456  CountedPtr< Scantable > out = getScantable(in, false);
1457  Table& tout = out->table();
1458  TableIterator iter(tout, "FREQ_ID");
1459  FFTServer<Float,Complex> ffts;
1460  while ( !iter.pastEnd() ) {
1461    Table tab = iter.table();
1462    Double rp,rv,inc;
1463    ROTableRow row(tab);
1464    const TableRecord& rec = row.get(0);
1465    uInt freqid = rec.asuInt("FREQ_ID");
1466    out->frequencies().getEntry(rp, rv, inc, freqid);
1467    ArrayColumn<Float> specCol(tab, "SPECTRA");
1468    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1469    for (int i=0; i<int(tab.nrow()); ++i) {
1470      Vector<Float> spec = specCol(i);
1471      Vector<uChar> flag = flagCol(i);
1472      Int lag0 = Int(spec.nelements()*abs(inc)/(frequency+width)+0.5);
1473      Int lag1 = Int(spec.nelements()*abs(inc)/(frequency-width)+0.5);
1474      for (int k=0; k < flag.nelements(); ++k ) {
1475        if (flag[k] > 0) {
1476          spec[k] = 0.0;
1477        }
1478      }
1479      Vector<Complex> lags;
1480      ffts.fft0(lags, spec);
1481      Int start =  max(0, lag0);
1482      Int end =  min(Int(lags.nelements()-1), lag1);
1483      if (start == end) {
1484        lags[start] = Complex(0.0);
1485      } else {
1486        for (int j=start; j <=end ;++j) {
1487          lags[j] = Complex(0.0);
1488        }
1489      }
1490      ffts.fft0(spec, lags);
1491      specCol.put(i, spec);
1492    }
1493    ++iter;
1494  }
1495  return out;
1496}
Note: See TracBrowser for help on using the repository browser.