source: trunk/src/STMath.cpp @ 1143

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

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