source: trunk/src/STMath.cpp @ 977

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

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