source: trunk/src/STMath.cpp @ 1145

Last change on this file since 1145 was 1145, checked in by mar637, 18 years ago

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