source: trunk/src/STMath.cpp @ 985

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

excpetion message update for alignFrequencies
  • 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("Can't convert as no output frame has been set"
1033                    " (use set_freqframe) or it is aligned already."));
1034  }
1035  MFrequency::Types system = in->frequencies().getFrame();
1036  MVTime mvt(refEpoch.getValue());
1037  String epochout = mvt.string(MVTime::YMD) + String(" (") + refEpoch.getRefString() + String(")");
1038  ostringstream oss;
1039  oss << "Aligned at reference Epoch " << epochout
1040      << " in frame " << MFrequency::showType(system);
1041  pushLog(String(oss));
1042  // set up the iterator
1043  Block<String> cols(4);
1044  // select by constant direction
1045  cols[0] = String("SRCNAME");
1046  cols[1] = String("BEAMNO");
1047  // select by IF ( no of channels varies over this )
1048  cols[2] = String("IFNO");
1049  // select by restfrequency
1050  cols[3] = String("MOLECULE_ID");
1051  TableIterator iter(tout, cols);
1052  while ( !iter.pastEnd() ) {
1053    Table t = iter.table();
1054    MDirection::ROScalarColumn dirCol(t, "DIRECTION");
1055    TableIterator fiter(t, "FREQ_ID");
1056    // determine nchan from the first row. This should work as
1057    // we are iterating over BEAMNO and IFNO    // we should have constant direction
1058
1059    ROArrayColumn<Float> sCol(t, "SPECTRA");
1060    MDirection direction = dirCol(0);
1061    uInt nchan = sCol(0).nelements();
1062    while ( !fiter.pastEnd() ) {
1063      Table ftab = fiter.table();
1064      ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
1065      // get the SpectralCoordinate for the freqid, which we are iterating over
1066      SpectralCoordinate sC = in->frequencies().getSpectralCoordinate(freqidCol(0));
1067      FrequencyAligner<Float> fa( sC, nchan, refEpoch,
1068                                  direction, refPos, system );
1069      // realign the SpectralCoordinate and put into the output Scantable
1070      Vector<String> units(1);
1071      units = String("Hz");
1072      Bool linear=True;
1073      SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
1074      sc2.setWorldAxisUnits(units);
1075      uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
1076                                            sc2.referenceValue()[0],
1077                                            sc2.increment()[0]);
1078      TableVector<uInt> tvec(ftab, "FREQ_ID");
1079      tvec = id;
1080      // create the "global" abcissa for alignment with same FREQ_ID
1081      Vector<Double> abc(nchan);
1082      Double w;
1083      for (uInt i=0; i<nchan; i++) {
1084        sC.toWorld(w,Double(i));
1085        abc[i] = w;
1086      }
1087      // cache abcissa for same time stamps, so iterate over those
1088      TableIterator timeiter(ftab, "TIME");
1089      while ( !timeiter.pastEnd() ) {
1090        Table tab = timeiter.table();
1091        ArrayColumn<Float> specCol(tab, "SPECTRA");
1092        ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1093        MEpoch::ROScalarColumn timeCol(tab, "TIME");
1094        // use align abcissa cache after the first row
1095        bool first = true;
1096        // these rows should be just be POLNO
1097        for (int i=0; i<tab.nrow(); ++i) {
1098          // input values
1099          Vector<uChar> flag = flagCol(i);
1100          Vector<Bool> mask(flag.shape());
1101          Vector<Float> specOut, spec;
1102          spec  = specCol(i);
1103          Vector<Bool> maskOut;Vector<uChar> flagOut;
1104          convertArray(mask, flag);
1105          // alignment
1106          Bool ok = fa.align(specOut, maskOut, abc, spec,
1107                             mask, timeCol(i), !first,
1108                             interp, False);
1109          // back into scantable
1110          flagOut.resize(maskOut.nelements());
1111          convertArray(flagOut, maskOut);
1112          flagCol.put(i, flagOut);
1113          specCol.put(i, specOut);
1114          // start abcissa caching
1115          first = false;
1116        }
1117        // next timestamp
1118        ++timeiter;
1119      }
1120      // next FREQ_ID
1121      ++fiter;
1122    }
1123    // next aligner
1124    ++iter;
1125  }
1126  // set this afterwards to ensure we are doing insitu correctly.
1127  out->frequencies().setFrame(system, true);
1128  return out;
1129}
Note: See TracBrowser for help on using the repository browser.