source: trunk/src/STMath.cpp @ 1066

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

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