source: trunk/src/STMath.cpp @ 917

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

0-based re-indexing SCANNO on merge. added frequencyAlign. Compiles, but not yet tested.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 34.8 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                 bool alignfreq)
67{
68  if ( avmode == "SCAN" && in.size() != 1 )
69    throw(AipsError("Can't perform 'SCAN' averaging on multiple tables"));
70  WeightType wtype = stringToWeight(weight);
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
273  for (uInt i=0; i < tout.nrow(); ++i) {
274    const TableRecord& rec = row.get(i);
275    Double ontime = rec.asDouble("TIME");
276    ROScalarColumn<Double> offtimeCol(offs, "TIME");
277    Double mindeltat = min(abs(offtimeCol.getColumn() - ontime));
278    Table sel = offs( abs(offs.col("TIME")-ontime) <= mindeltat
279                       && offs.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
280                       && offs.col("IFNO") == Int(rec.asuInt("IFNO"))
281                       && offs.col("POLNO") == Int(rec.asuInt("POLNO")) );
282
283    TableRow offrow(sel);
284    const TableRecord& offrec = offrow.get(0);//should only be one row
285    RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
286    RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
287    RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
288    /// @fixme this assumes tsys is a scalar not vector
289    Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
290    Vector<Float> specon, tsyson;
291    outtsysCol.get(i, tsyson);
292    outspecCol.get(i, specon);
293    Vector<uChar> flagon;
294    outflagCol.get(i, flagon);
295    MaskedArray<Float> mon = maskedArray(specon, flagon);
296    MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
297    MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
298    if (preserve) {
299      quot -= tsysoffscalar;
300    } else {
301      quot -= tsyson[0];
302    }
303    outspecCol.put(i, quot.getArray());
304    outflagCol.put(i, flagsFromMA(quot));
305  }
306  return out;
307}
308
309CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
310{
311  // make copy or reference
312  CountedPtr< Scantable > out = getScantable(in, false);
313  Table& tout = out->table();
314  Block<String> cols(3);
315  cols[0] = String("SCANNO");
316  cols[1] = String("BEAMNO");
317  cols[2] = String("POLNO");
318  TableIterator iter(tout, cols);
319  while (!iter.pastEnd()) {
320    Table subt = iter.table();
321    // this should leave us with two rows for the two IFs....if not ignore
322    if (subt.nrow() != 2 ) {
323      continue;
324    }
325    ArrayColumn<Float> specCol(tout, "SPECTRA");
326    ArrayColumn<Float> tsysCol(tout, "TSYS");
327    ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
328    Vector<Float> onspec,offspec, ontsys, offtsys;
329    Vector<uChar> onflag, offflag;
330    tsysCol.get(0, ontsys);   tsysCol.get(1, offtsys);
331    specCol.get(0, onspec);   specCol.get(1, offspec);
332    flagCol.get(0, onflag);   flagCol.get(1, offflag);
333    MaskedArray<Float> on  = maskedArray(onspec, onflag);
334    MaskedArray<Float> off = maskedArray(offspec, offflag);
335    MaskedArray<Float> oncopy = on.copy();
336
337    on /= off; on -= 1.0f;
338    on *= ontsys[0];
339    off /= oncopy; off -= 1.0f;
340    off *= offtsys[0];
341    specCol.put(0, on.getArray());
342    const Vector<Bool>& m0 = on.getMask();
343    Vector<uChar> flags0(m0.shape());
344    convertArray(flags0, !m0);
345    flagCol.put(0, flags0);
346
347    specCol.put(1, off.getArray());
348    const Vector<Bool>& m1 = off.getMask();
349    Vector<uChar> flags1(m1.shape());
350    convertArray(flags1, !m1);
351    flagCol.put(1, flags1);
352    ++iter;
353  }
354
355  return out;
356}
357
358std::vector< float > STMath::statistic( const CountedPtr< Scantable > & in,
359                                        const std::vector< bool > & mask,
360                                        const std::string& which )
361{
362
363  Vector<Bool> m(mask);
364  const Table& tab = in->table();
365  ROArrayColumn<Float> specCol(tab, "SPECTRA");
366  ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
367  std::vector<float> out;
368  for (uInt i=0; i < tab.nrow(); ++i ) {
369    Vector<Float> spec; specCol.get(i, spec);
370    Vector<uChar> flag; flagCol.get(i, flag);
371    MaskedArray<Float> ma  = maskedArray(spec, flag);
372    float outstat = 0.0;
373    if ( spec.nelements() == m.nelements() ) {
374      outstat = mathutil::statistics(which, ma(m));
375    } else {
376      outstat = mathutil::statistics(which, ma);
377    }
378    out.push_back(outstat);
379  }
380  return out;
381}
382
383CountedPtr< Scantable > STMath::bin( const CountedPtr< Scantable > & in,
384                                     int width )
385{
386  if ( !in->getSelection().empty() ) throw(AipsError("Can't bin subset of the data."));
387  CountedPtr< Scantable > out = getScantable(in, false);
388  Table& tout = out->table();
389  out->frequencies().rescale(width, "BIN");
390  ArrayColumn<Float> specCol(tout, "SPECTRA");
391  ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
392  for (uInt i=0; i < tout.nrow(); ++i ) {
393    MaskedArray<Float> main  = maskedArray(specCol(i), flagCol(i));
394    MaskedArray<Float> maout;
395    LatticeUtilities::bin(maout, main, 0, Int(width));
396    /// @todo implement channel based tsys binning
397    specCol.put(i, maout.getArray());
398    flagCol.put(i, flagsFromMA(maout));
399    // take only the first binned spectrum's length for the deprecated
400    // global header item nChan
401    if (i==0) tout.rwKeywordSet().define(String("nChan"),
402                                       Int(maout.getArray().nelements()));
403  }
404  return out;
405}
406
407CountedPtr< Scantable > STMath::resample( const CountedPtr< Scantable >& in,
408                                          const std::string& method,
409                                          float width )
410//
411// Should add the possibility of width being specified in km/s. This means
412// that for each freqID (SpectralCoordinate) we will need to convert to an
413// average channel width (say at the reference pixel).  Then we would need
414// to be careful to make sure each spectrum (of different freqID)
415// is the same length.
416//
417{
418  InterpolateArray1D<Double,Float>::InterpolationMethod interp;
419  Int interpMethod(stringToIMethod(method));
420
421  CountedPtr< Scantable > out = getScantable(in, false);
422  Table& tout = out->table();
423
424// Resample SpectralCoordinates (one per freqID)
425  out->frequencies().rescale(width, "RESAMPLE");
426  TableIterator iter(tout, "IFNO");
427  TableRow row(tout);
428  while ( !iter.pastEnd() ) {
429    Table tab = iter.table();
430    ArrayColumn<Float> specCol(tab, "SPECTRA");
431    //ArrayColumn<Float> tsysCol(tout, "TSYS");
432    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
433    Vector<Float> spec;
434    Vector<uChar> flag;
435    specCol.get(0,spec); // the number of channels should be constant per IF
436    uInt nChanIn = spec.nelements();
437    Vector<Float> xIn(nChanIn); indgen(xIn);
438    Int fac =  Int(nChanIn/width);
439    Vector<Float> xOut(fac+10); // 10 to be safe - resize later
440    uInt k = 0;
441    Float x = 0.0;
442    while (x < Float(nChanIn) ) {
443      xOut(k) = x;
444      k++;
445      x += width;
446    }
447    uInt nChanOut = k;
448    xOut.resize(nChanOut, True);
449    // process all rows for this IFNO
450    Vector<Float> specOut;
451    Vector<Bool> maskOut;
452    Vector<uChar> flagOut;
453    for (uInt i=0; i < tab.nrow(); ++i) {
454      specCol.get(i, spec);
455      flagCol.get(i, flag);
456      Vector<Bool> mask(flag.nelements());
457      convertArray(mask, flag);
458
459      IPosition shapeIn(spec.shape());
460      //sh.nchan = nChanOut;
461      InterpolateArray1D<Float,Float>::interpolate(specOut, maskOut, xOut,
462                                                   xIn, spec, mask,
463                                                   interpMethod, True, True);
464      /// @todo do the same for channel based Tsys
465      flagOut.resize(maskOut.nelements());
466      convertArray(flagOut, maskOut);
467      specCol.put(i, specOut);
468      flagCol.put(i, flagOut);
469    }
470    ++iter;
471  }
472
473  return out;
474}
475
476STMath::imethod STMath::stringToIMethod(const std::string& in)
477{
478  static STMath::imap lookup;
479
480  // initialize the lookup table if necessary
481  if ( lookup.empty() ) {
482    lookup["NEAR"]   = InterpolateArray1D<Double,Float>::nearestNeighbour;
483    lookup["LIN"] = InterpolateArray1D<Double,Float>::linear;
484    lookup["CUB"]  = InterpolateArray1D<Double,Float>::cubic;
485    lookup["SPL"]  = InterpolateArray1D<Double,Float>::spline;
486  }
487
488  STMath::imap::const_iterator iter = lookup.find(in);
489
490  if ( lookup.end() == iter ) {
491    std::string message = in;
492    message += " is not a valid interpolation mode";
493    throw(AipsError(message));
494  }
495  return iter->second;
496}
497
498WeightType STMath::stringToWeight(const std::string& in)
499{
500  static std::map<std::string, WeightType> lookup;
501
502  // initialize the lookup table if necessary
503  if ( lookup.empty() ) {
504    lookup["NONE"]   = asap::NONE;
505    lookup["TINT"] = asap::TINT;
506    lookup["TINTSYS"]  = asap::TINTSYS;
507    lookup["TSYS"]  = asap::TSYS;
508    lookup["VAR"]  = asap::VAR;
509  }
510
511  std::map<std::string, WeightType>::const_iterator iter = lookup.find(in);
512
513  if ( lookup.end() == iter ) {
514    std::string message = in;
515    message += " is not a valid weighting mode";
516    throw(AipsError(message));
517  }
518  return iter->second;
519}
520
521CountedPtr< Scantable > STMath::gainElevation( const CountedPtr< Scantable >& in,
522                                               const vector< float > & coeff,
523                                               const std::string & filename,
524                                               const std::string& method)
525{
526  // Get elevation data from Scantable and convert to degrees
527  CountedPtr< Scantable > out = getScantable(in, false);
528  Table& tab = in->table();
529  ROScalarColumn<Float> elev(tab, "ELEVATION");
530  Vector<Float> x = elev.getColumn();
531  x *= Float(180 / C::pi);                        // Degrees
532
533  Vector<Float> coeffs(coeff);
534  const uInt nc = coeffs.nelements();
535  if ( filename.length() > 0 && nc > 0 ) {
536    throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
537  }
538
539  // Correct
540  if ( nc > 0 || filename.length() == 0 ) {
541    // Find instrument
542    Bool throwit = True;
543    Instrument inst =
544      STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
545                                throwit);
546
547    // Set polynomial
548    Polynomial<Float>* ppoly = 0;
549    Vector<Float> coeff;
550    String msg;
551    if ( nc > 0 ) {
552      ppoly = new Polynomial<Float>(nc);
553      coeff = coeffs;
554      msg = String("user");
555    } else {
556      STAttr sdAttr;
557      coeff = sdAttr.gainElevationPoly(inst);
558      ppoly = new Polynomial<Float>(3);
559      msg = String("built in");
560    }
561
562    if ( coeff.nelements() > 0 ) {
563      ppoly->setCoefficients(coeff);
564    } else {
565      delete ppoly;
566      throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
567    }
568    ostringstream oss;
569    oss << "Making polynomial correction with " << msg << " coefficients:" << endl;
570    oss << "   " <<  coeff;
571    pushLog(String(oss));
572    const uInt nrow = tab.nrow();
573    Vector<Float> factor(nrow);
574    for ( uInt i=0; i < nrow; ++i ) {
575      factor[i] = 1.0 / (*ppoly)(x[i]);
576    }
577    delete ppoly;
578    scaleByVector(tab, factor, true);
579
580  } else {
581    // Read and correct
582    pushLog("Making correction from ascii Table");
583    scaleFromAsciiTable(tab, filename, method, x, true);
584  }
585  return out;
586}
587
588void STMath::scaleFromAsciiTable(Table& in, const std::string& filename,
589                                 const std::string& method,
590                                 const Vector<Float>& xout, bool dotsys)
591{
592
593// Read gain-elevation ascii file data into a Table.
594
595  String formatString;
596  Table tbl = readAsciiTable(formatString, Table::Memory, filename, "", "", False);
597  scaleFromTable(in, tbl, method, xout, dotsys);
598}
599
600void STMath::scaleFromTable(Table& in,
601                            const Table& table,
602                            const std::string& method,
603                            const Vector<Float>& xout, bool dotsys)
604{
605
606  ROScalarColumn<Float> geElCol(table, "ELEVATION");
607  ROScalarColumn<Float> geFacCol(table, "FACTOR");
608  Vector<Float> xin = geElCol.getColumn();
609  Vector<Float> yin = geFacCol.getColumn();
610  Vector<Bool> maskin(xin.nelements(),True);
611
612  // Interpolate (and extrapolate) with desired method
613
614   //InterpolateArray1D<Double,Float>::InterpolationMethod method;
615   Int intmethod(stringToIMethod(method));
616
617   Vector<Float> yout;
618   Vector<Bool> maskout;
619   InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout,
620                                                xin, yin, maskin, intmethod,
621                                                True, True);
622
623   scaleByVector(in, Float(1.0)/yout, dotsys);
624}
625
626void STMath::scaleByVector( Table& in,
627                            const Vector< Float >& factor,
628                            bool dotsys )
629{
630  uInt nrow = in.nrow();
631  if ( factor.nelements() != nrow ) {
632    throw(AipsError("factors.nelements() != table.nelements()"));
633  }
634  ArrayColumn<Float> specCol(in, "SPECTRA");
635  ArrayColumn<uChar> flagCol(in, "FLAGTRA");
636  ArrayColumn<Float> tsysCol(in, "TSYS");
637  for (uInt i=0; i < nrow; ++i) {
638    MaskedArray<Float> ma  = maskedArray(specCol(i), flagCol(i));
639    ma *= factor[i];
640    specCol.put(i, ma.getArray());
641    flagCol.put(i, flagsFromMA(ma));
642    if ( dotsys ) {
643      Vector<Float> tsys;
644      tsysCol.get(i, tsys);
645      tsys *= factor[i];
646      specCol.put(i,tsys);
647    }
648  }
649}
650
651CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in,
652                                             float d, float etaap,
653                                             float jyperk )
654{
655  CountedPtr< Scantable > out = getScantable(in, false);
656  Table& tab = in->table();
657  Unit fluxUnit(tab.keywordSet().asString("FluxUnit"));
658  Unit K(String("K"));
659  Unit JY(String("Jy"));
660
661  bool tokelvin = true;
662  Double cfac = 1.0;
663
664  if ( fluxUnit == JY ) {
665    pushLog("Converting to K");
666    Quantum<Double> t(1.0,fluxUnit);
667    Quantum<Double> t2 = t.get(JY);
668    cfac = (t2 / t).getValue();               // value to Jy
669
670    tokelvin = true;
671    out->setFluxUnit("K");
672  } else if ( fluxUnit == K ) {
673    pushLog("Converting to Jy");
674    Quantum<Double> t(1.0,fluxUnit);
675    Quantum<Double> t2 = t.get(K);
676    cfac = (t2 / t).getValue();              // value to K
677
678    tokelvin = false;
679    out->setFluxUnit("Jy");
680  } else {
681    throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
682  }
683  // Make sure input values are converted to either Jy or K first...
684  Float factor = cfac;
685
686  // Select method
687  if (jyperk > 0.0) {
688    factor *= jyperk;
689    if ( tokelvin ) factor = 1.0 / jyperk;
690    ostringstream oss;
691    oss << "Jy/K = " << jyperk;
692    pushLog(String(oss));
693    Vector<Float> factors(tab.nrow(), factor);
694    scaleByVector(tab,factors, false);
695  } else if ( etaap > 0.0) {
696    Instrument inst =
697      STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"), True);
698    STAttr sda;
699    if (d < 0) d = sda.diameter(inst);
700    Float jyPerk = STAttr::findJyPerK(etaap, d);
701    ostringstream oss;
702    oss << "Jy/K = " << jyperk;
703    pushLog(String(oss));
704    factor *= jyperk;
705    if ( tokelvin ) {
706      factor = 1.0 / factor;
707    }
708    Vector<Float> factors(tab.nrow(), factor);
709    scaleByVector(tab, factors, False);
710  } else {
711
712    // OK now we must deal with automatic look up of values.
713    // We must also deal with the fact that the factors need
714    // to be computed per IF and may be different and may
715    // change per integration.
716
717    pushLog("Looking up conversion factors");
718    convertBrightnessUnits(out, tokelvin, cfac);
719  }
720
721  return out;
722}
723
724void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in,
725                                     bool tokelvin, float cfac )
726{
727  Table& table = in->table();
728  Instrument inst =
729    STAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True);
730  TableIterator iter(table, "FREQ_ID");
731  STFrequencies stfreqs = in->frequencies();
732  STAttr sdAtt;
733  while (!iter.pastEnd()) {
734    Table tab = iter.table();
735    ArrayColumn<Float> specCol(tab, "SPECTRA");
736    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
737    ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID");
738    MEpoch::ROScalarColumn timeCol(tab, "TIME");
739
740    uInt freqid; freqidCol.get(0, freqid);
741    Vector<Float> tmpspec; specCol.get(0, tmpspec);
742    // STAttr.JyPerK has a Vector interface... change sometime.
743    Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements()));
744    for ( uInt i=0; i<tab.nrow(); ++i) {
745      Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0];
746      Float factor = cfac * jyperk;
747      if ( tokelvin ) factor = Float(1.0) / factor;
748      MaskedArray<Float> ma  = maskedArray(specCol(i), flagCol(i));
749      ma *= factor;
750      specCol.put(i, ma.getArray());
751      flagCol.put(i, flagsFromMA(ma));
752    }
753  ++iter;
754  }
755}
756
757CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
758                                         float tau )
759{
760  CountedPtr< Scantable > out = getScantable(in, false);
761  Table& tab = in->table();
762  ROScalarColumn<Float> elev(tab, "ELEVATION");
763  ArrayColumn<Float> specCol(tab, "SPECTRA");
764  ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
765  for ( uInt i=0; i<tab.nrow(); ++i) {
766    Float zdist = Float(C::pi_2) - elev(i);
767    Float factor = exp(tau)/cos(zdist);
768    MaskedArray<Float> ma  = maskedArray(specCol(i), flagCol(i));
769    ma *= factor;
770    specCol.put(i, ma.getArray());
771    flagCol.put(i, flagsFromMA(ma));
772  }
773  return out;
774}
775
776CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
777                                        const std::string& kernel, float width )
778{
779  CountedPtr< Scantable > out = getScantable(in, false);
780  Table& table = in->table();
781  VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
782  // same IFNO should have same no of channels
783  // this saves overhead
784  TableIterator iter(table, "IFNO");
785  while (!iter.pastEnd()) {
786    Table tab = iter.table();
787    ArrayColumn<Float> specCol(tab, "SPECTRA");
788    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
789    Vector<Float> tmpspec; specCol.get(0, tmpspec);
790    uInt nchan = tmpspec.nelements();
791    Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
792    Convolver<Float> conv(kvec, IPosition(1,nchan));
793    Vector<Float> spec;
794    Vector<uChar> flag;
795    for ( uInt i=0; i<tab.nrow(); ++i) {
796      specCol.get(i, spec);
797      flagCol.get(i, flag);
798      Vector<Bool> mask(flag.nelements());
799      convertArray(mask, flag);
800      Vector<Float> specout;
801      if ( type == VectorKernel::HANNING ) {
802        Vector<Bool> maskout;
803        mathutil::hanning(specout, maskout, spec , mask);
804        convertArray(flag, maskout);
805        flagCol.put(i, flag);
806      } else {
807        mathutil::replaceMaskByZero(specout, mask);
808        conv.linearConv(specout, spec);
809      }
810      specCol.put(i, specout);
811    }
812    ++iter;
813  }
814  return out;
815}
816
817CountedPtr< Scantable >
818  STMath::merge( const std::vector< CountedPtr < Scantable > >& in )
819{
820  if ( in.size() < 2 ) {
821    throw(AipsError("Need at least two scantables to perform a merge."));
822  }
823  std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
824  bool insitu = insitu_;
825  setInsitu(false);
826  CountedPtr< Scantable > out = getScantable(*it, false);
827  setInsitu(insitu);
828  Table& tout = out->table();
829  ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
830  ScalarColumn<uInt> scannocol(tout,"SCANNO"), focusidcol(tout,"FOCUS_ID");
831  // Renumber SCANNO to be 0-based
832  uInt offset = min(scannocol.getColumn());
833  TableVector<uInt> scannos(tout, "SCANNO");
834  scannos -= offset;
835  uInt newscanno = max(scannocol.getColumn())+1;
836  ++it;
837  while ( it != in.end() ){
838    if ( ! (*it)->conformant(*out) ) {
839      // log message: "ignoring scantable i, as it isn't
840      // conformant with the other(s)"
841      cerr << "oh oh" << endl;
842      ++it;
843      continue;
844    }
845    out->appendToHistoryTable((*it)->history());
846    const Table& tab = (*it)->table();
847    TableIterator scanit(tab, "SCANNO");
848    while (!scanit.pastEnd()) {
849      TableIterator freqit(scanit.table(), "FREQ_ID");
850      while ( !freqit.pastEnd() ) {
851        Table thetab = freqit.table();
852        uInt nrow = tout.nrow();
853        //tout.addRow(thetab.nrow());
854        TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
855        ROTableRow row(thetab);
856        for ( uInt i=0; i<thetab.nrow(); ++i) {
857          uInt k = nrow+i;
858          scannocol.put(k, newscanno);
859          const TableRecord& rec = row.get(i);
860          Double rv,rp,inc;
861          (*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
862          uInt id;
863          id = out->frequencies().addEntry(rp, rv, inc);
864          freqidcol.put(k,id);
865          String name,fname;Double rf;
866          (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
867          id = out->molecules().addEntry(rf, name, fname);
868          molidcol.put(k, id);
869          Float frot,fang,ftan;
870          (*it)->focus().getEntry(frot, fang, ftan, rec.asuInt("FOCUS_ID"));
871          id = out->focus().addEntry(frot, fang, ftan);
872          focusidcol.put(k, id);
873        }
874        ++freqit;
875      }
876      ++newscanno;
877      ++scanit;
878    }
879    ++it;
880  }
881  return out;
882}
883
884CountedPtr< Scantable >
885  STMath::invertPhase( const CountedPtr < Scantable >& in )
886{
887  applyToPol(in, &STPol::invertPhase, Float(0.0));
888}
889
890CountedPtr< Scantable >
891  STMath::rotateXYPhase( const CountedPtr < Scantable >& in, float phase )
892{
893   return applyToPol(in, &STPol::rotatePhase, Float(phase));
894}
895
896CountedPtr< Scantable >
897  STMath::rotateLinPolPhase( const CountedPtr < Scantable >& in, float phase )
898{
899  return applyToPol(in, &STPol::rotateLinPolPhase, Float(phase));
900}
901
902CountedPtr< Scantable > STMath::applyToPol( const CountedPtr<Scantable>& in,
903                                             STPol::polOperation fptr,
904                                             Float phase )
905{
906  CountedPtr< Scantable > out = getScantable(in, false);
907  Table& tout = out->table();
908  Block<String> cols(4);
909  cols[0] = String("SCANNO");
910  cols[1] = String("BEAMNO");
911  cols[2] = String("IFNO");
912  cols[3] = String("CYCLENO");
913  TableIterator iter(tout, cols);
914  STPol* stpol = NULL;
915  stpol =STPol::getPolClass(out->factories_, out->getPolType() );
916  while (!iter.pastEnd()) {
917    Table t = iter.table();
918    ArrayColumn<Float> speccol(t, "SPECTRA");
919    Matrix<Float> pols = speccol.getColumn();
920    Matrix<Float> out;
921    try {
922      stpol->setSpectra(pols);
923      (stpol->*fptr)(phase);
924      speccol.putColumn(stpol->getSpectra());
925      delete stpol;stpol=0;
926    } catch (AipsError& e) {
927      delete stpol;stpol=0;
928      throw(e);
929    }
930    ++iter;
931  }
932  return out;
933}
934
935CountedPtr< Scantable >
936  STMath::swapPolarisations( const CountedPtr< Scantable > & in )
937{
938  CountedPtr< Scantable > out = getScantable(in, false);
939  Table& tout = out->table();
940  Table t0 = tout(tout.col("POLNO") == 0);
941  Table t1 = tout(tout.col("POLNO") == 1);
942  if ( t0.nrow() != t1.nrow() )
943    throw(AipsError("Inconsistent number of polarisations"));
944  ArrayColumn<Float> speccol0(t0, "SPECTRA");
945  ArrayColumn<uChar> flagcol0(t0, "FLAGTRA");
946  ArrayColumn<Float> speccol1(t1, "SPECTRA");
947  ArrayColumn<uChar> flagcol1(t1, "FLAGTRA");
948  Matrix<Float> s0 = speccol0.getColumn();
949  Matrix<uChar> f0 = flagcol0.getColumn();
950  speccol0.putColumn(speccol1.getColumn());
951  flagcol0.putColumn(flagcol1.getColumn());
952  speccol1.putColumn(s0);
953  flagcol1.putColumn(f0);
954  return out;
955}
956
957CountedPtr< Scantable >
958  asap::STMath::frequencyAlign( const CountedPtr< Scantable > & in,
959                                const std::string & refTime,
960                                const std::string & method, bool perfreqid )
961{
962  CountedPtr< Scantable > out = getScantable(in, false);
963  Table& tout = out->table();
964  // clear ouput frequency table
965  Table ftable = out->frequencies().table();
966  ftable.removeRow(ftable.rowNumbers());
967  // Get reference Epoch to time of first row or given String
968  Unit DAY(String("d"));
969  MEpoch::Ref epochRef(in->getTimeReference());
970  MEpoch refEpoch;
971  if (refTime.length()>0) {
972    Quantum<Double> qt;
973    if (MVTime::read(qt,refTime)) {
974      MVEpoch mv(qt);
975      refEpoch = MEpoch(mv, epochRef);
976   } else {
977      throw(AipsError("Invalid format for Epoch string"));
978   }
979  } else {
980    refEpoch = in->timeCol_(0);
981  }
982  MPosition refPos = in->getAntennaPosition();
983/*  ostringstream oss;
984  oss << "Aligned at reference Epoch " << formatEpoch(refEpoch)
985      << " in frame " << MFrequency::showType(freqSystem);
986  pushLog(String(oss));
987*/
988  InterpolateArray1D<Double,Float>::InterpolationMethod interp;
989  Int interpMethod(stringToIMethod(method));
990  // test if user frame is different to base frame
991  if ( in->frequencies().getFrameString(true)
992       == in->frequencies().getFrameString(false) ) {
993    throw(AipsError("You have not set a frequency frame different from the initial - use function set_freqframe"));
994  }
995  MFrequency::Types system = in->frequencies().getFrame();
996  // set up the iterator
997  Block<String> cols(3);
998  // select the by constant direction
999  cols[0] = String("SRCNAME");
1000  cols[1] = String("BEAMNO");
1001  // select by IF ( no of channels varies over this )
1002  cols[2] = String("IFNO");
1003  // handle freqid based alignment
1004  if ( perfreqid ) {
1005    cols.resize(4);
1006    cols[3] = String("FREQ_ID");
1007  }
1008  TableIterator iter(tout, cols);
1009  while (!iter.pastEnd()) {
1010    Table t = iter.table();
1011    MDirection::ROScalarColumn dirCol(t, "DIRECTION");
1012    ScalarColumn<uInt> freqidCol(t, "FREQ_ID");
1013    // get the SpectralCoordinate for the freqid, which we are iterating over
1014    SpectralCoordinate sC = in->frequencies().getSpectralCoordinate(freqidCol(0));
1015    // determine nchan from the first row. This should work as
1016    // we are iterting over IFNO
1017    ROArrayColumn<Float> sCol(t, "SPECTRA");
1018    uInt nchan = sCol(0).nelements();
1019    // we should have constant direction
1020    FrequencyAligner<Float> fa(sC, nchan, refEpoch, dirCol(0), refPos, system);
1021    // realign the SpectralCOordinate and put into the output Scantable
1022    Vector<String> units(1);
1023    units = String("Hz");
1024    Bool linear=True;
1025    SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
1026    sc2.setWorldAxisUnits(units);
1027    out->frequencies().addEntry(sc2.referencePixel()[0],
1028                                sc2.referenceValue()[0],
1029                                sc2.increment()[0]);
1030    // create the abcissa for IF based alignment
1031    Vector<Double> abc;
1032    if ( !perfreqid ) {
1033      abc.resize(nchan);
1034      Double w;
1035      for (uInt i=0; i<nchan; i++) {
1036        sC.toWorld(w,Double(i));
1037        abc[i] = w;
1038      }
1039    }
1040    // cache abcissa for same time stamps, so iterate over those
1041    TableIterator timeiter(t, "TIME");
1042    while ( !timeiter.pastEnd() ) {
1043      Table tab = timeiter.table();
1044      ArrayColumn<Float> specCol(tab, "SPECTRA");
1045      ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1046      MEpoch::ROScalarColumn timeCol(tab, "TIME");
1047      // use align abcissa cache after the first row
1048      bool first = true;
1049      // these rows should be just be POLNO
1050      for (int i=0; i<tab.nrow(); ++i) {
1051        // input values
1052        Vector<uChar> flag = flagCol(i);
1053        Vector<Bool> mask(flag.shape());
1054        Vector<Float> specOut;
1055        Vector<Bool> maskOut;Vector<uChar> flagOut;
1056        convertArray(mask, flag);
1057        // alignment
1058        if ( perfreqid ) {
1059          Bool ok = fa.align(specOut, maskOut, specCol(i),
1060                        mask, timeCol(i), !first,
1061                        interp, False);
1062        } else {
1063          Bool ok = fa.align(specOut, maskOut, abc, specCol(i),
1064                        mask, timeCol(i), !first,
1065                        interp, False);
1066        }
1067        // back into scantable
1068        flagOut.resize(maskOut.nelements());
1069        convertArray(flagOut, maskOut);
1070        flagCol.put(i, flagOut);
1071        specCol.put(i, specOut);
1072        // start ancissa caching
1073        first = false;
1074      }
1075      ++timeiter;
1076    }
1077    // next aligner
1078    ++iter;
1079  }
1080  return out;
1081}
Note: See TracBrowser for help on using the repository browser.