source: trunk/src/STMath.cpp @ 1069

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

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