source: trunk/src/STMath.cpp @ 1104

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

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