source: trunk/src/STMath.cpp @ 902

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

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