source: trunk/src/STMath.cpp @ 862

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

added history append to average and merge
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 27.6 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 "SDAttr.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
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    MaskedArray<Float> ma  = maskedArray(spec, flagCol(i));
365    float outstat;
366    if ( spec.nelements() == m.nelements() ) {
367      outstat = mathutil::statistics(which, ma(m));
368    } else {
369      outstat = mathutil::statistics(which, ma);
370    }
371    out.push_back(outstat);
372  }
373  return out;
374}
375
376CountedPtr< Scantable > STMath::bin( const CountedPtr< Scantable > & in,
377                                     int width )
378{
379  if ( !in->getSelection().empty() ) throw(AipsError("Can't bin subset of the data."));
380  CountedPtr< Scantable > out = getScantable(in, false);
381  Table& tout = out->table();
382  out->frequencies().rescale(width, "BIN");
383  ArrayColumn<Float> specCol(tout, "SPECTRA");
384  ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
385  for (uInt i=0; i < tout.nrow(); ++i ) {
386    MaskedArray<Float> main  = maskedArray(specCol(i), flagCol(i));
387    MaskedArray<Float> maout;
388    LatticeUtilities::bin(maout, main, 0, Int(width));
389    /// @todo implement channel based tsys binning
390    specCol.put(i, maout.getArray());
391    flagCol.put(i, flagsFromMA(maout));
392    // take only the first binned spectrum's length for the deprecated
393    // global header item nChan
394    if (i==0) tout.rwKeywordSet().define(String("nChan"),
395                                       Int(maout.getArray().nelements()));
396  }
397  return out;
398}
399
400CountedPtr< Scantable > STMath::resample( const CountedPtr< Scantable >& in,
401                                          const std::string& method,
402                                          float width )
403//
404// Should add the possibility of width being specified in km/s. This means
405// that for each freqID (SpectralCoordinate) we will need to convert to an
406// average channel width (say at the reference pixel).  Then we would need
407// to be careful to make sure each spectrum (of different freqID)
408// is the same length.
409//
410{
411  InterpolateArray1D<Double,Float>::InterpolationMethod interp;
412  Int interpMethod(stringToIMethod(method));
413
414  CountedPtr< Scantable > out = getScantable(in, false);
415  Table& tout = out->table();
416
417// Resample SpectralCoordinates (one per freqID)
418  out->frequencies().rescale(width, "RESAMPLE");
419  TableIterator iter(tout, "IFNO");
420  TableRow row(tout);
421  while ( !iter.pastEnd() ) {
422    Table tab = iter.table();
423    ArrayColumn<Float> specCol(tab, "SPECTRA");
424    //ArrayColumn<Float> tsysCol(tout, "TSYS");
425    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
426    Vector<Float> spec;
427    Vector<uChar> flag;
428    specCol.get(0,spec); // the number of channels should be constant per IF
429    uInt nChanIn = spec.nelements();
430    Vector<Float> xIn(nChanIn); indgen(xIn);
431    Int fac =  Int(nChanIn/width);
432    Vector<Float> xOut(fac+10); // 10 to be safe - resize later
433    uInt k = 0;
434    Float x = 0.0;
435    while (x < Float(nChanIn) ) {
436      xOut(k) = x;
437      k++;
438      x += width;
439    }
440    uInt nChanOut = k;
441    xOut.resize(nChanOut, True);
442    // process all rows for this IFNO
443    Vector<Float> specOut;
444    Vector<Bool> maskOut;
445    Vector<uChar> flagOut;
446    for (uInt i=0; i < tab.nrow(); ++i) {
447      specCol.get(i, spec);
448      flagCol.get(i, flag);
449      Vector<Bool> mask(flag.nelements());
450      convertArray(mask, flag);
451
452      IPosition shapeIn(spec.shape());
453      //sh.nchan = nChanOut;
454      InterpolateArray1D<Float,Float>::interpolate(specOut, maskOut, xOut,
455                                                   xIn, spec, mask,
456                                                   interpMethod, True, True);
457      /// @todo do the same for channel based Tsys
458      flagOut.resize(maskOut.nelements());
459      convertArray(flagOut, maskOut);
460      specCol.put(i, specOut);
461      flagCol.put(i, flagOut);
462    }
463    ++iter;
464  }
465
466  return out;
467}
468
469STMath::imethod STMath::stringToIMethod(const std::string& in)
470{
471  static STMath::imap lookup;
472
473  // initialize the lookup table if necessary
474  if ( lookup.empty() ) {
475    lookup["NEAR"]   = InterpolateArray1D<Double,Float>::nearestNeighbour;
476    lookup["LIN"] = InterpolateArray1D<Double,Float>::linear;
477    lookup["CUB"]  = InterpolateArray1D<Double,Float>::cubic;
478    lookup["SPL"]  = InterpolateArray1D<Double,Float>::spline;
479  }
480
481  STMath::imap::const_iterator iter = lookup.find(in);
482
483  if ( lookup.end() == iter ) {
484    std::string message = in;
485    message += " is not a valid interpolation mode";
486    throw(AipsError(message));
487  }
488  return iter->second;
489}
490
491WeightType STMath::stringToWeight(const std::string& in)
492{
493  static std::map<std::string, WeightType> lookup;
494
495  // initialize the lookup table if necessary
496  if ( lookup.empty() ) {
497    lookup["NONE"]   = asap::NONE;
498    lookup["TINT"] = asap::TINT;
499    lookup["TINTSYS"]  = asap::TINTSYS;
500    lookup["TSYS"]  = asap::TSYS;
501    lookup["VAR"]  = asap::VAR;
502  }
503
504  std::map<std::string, WeightType>::const_iterator iter = lookup.find(in);
505
506  if ( lookup.end() == iter ) {
507    std::string message = in;
508    message += " is not a valid weighting mode";
509    throw(AipsError(message));
510  }
511  return iter->second;
512}
513
514CountedPtr< Scantable > STMath::gainElevation( const CountedPtr< Scantable >& in,
515                                               const Vector< Float > & coeffs,
516                                               const std::string & filename,
517                                               const std::string& method)
518{
519  // Get elevation data from Scantable and convert to degrees
520  CountedPtr< Scantable > out = getScantable(in, false);
521  Table& tab = in->table();
522  ROScalarColumn<Float> elev(tab, "ELEVATION");
523  Vector<Float> x = elev.getColumn();
524  x *= Float(180 / C::pi);                        // Degrees
525
526  const uInt nc = coeffs.nelements();
527  if ( filename.length() > 0 && nc > 0 ) {
528    throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
529  }
530
531  // Correct
532  if ( nc > 0 || filename.length() == 0 ) {
533    // Find instrument
534    Bool throwit = True;
535    Instrument inst =
536      SDAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
537                                throwit);
538
539    // Set polynomial
540    Polynomial<Float>* ppoly = 0;
541    Vector<Float> coeff;
542    String msg;
543    if ( nc > 0 ) {
544      ppoly = new Polynomial<Float>(nc);
545      coeff = coeffs;
546      msg = String("user");
547    } else {
548      SDAttr sdAttr;
549      coeff = sdAttr.gainElevationPoly(inst);
550      ppoly = new Polynomial<Float>(3);
551      msg = String("built in");
552    }
553
554    if ( coeff.nelements() > 0 ) {
555      ppoly->setCoefficients(coeff);
556    } else {
557      delete ppoly;
558      throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
559    }
560    ostringstream oss;
561    oss << "Making polynomial correction with " << msg << " coefficients:" << endl;
562    oss << "   " <<  coeff;
563    pushLog(String(oss));
564    const uInt nrow = tab.nrow();
565    Vector<Float> factor(nrow);
566    for ( uInt i=0; i < nrow; ++i ) {
567      factor[i] = 1.0 / (*ppoly)(x[i]);
568    }
569    delete ppoly;
570    scaleByVector(tab, factor, true);
571
572  } else {
573    // Read and correct
574    pushLog("Making correction from ascii Table");
575    scaleFromAsciiTable(tab, filename, method, x, true);
576  }
577  return out;
578}
579
580void STMath::scaleFromAsciiTable(Table& in, const std::string& filename,
581                                 const std::string& method,
582                                 const Vector<Float>& xout, bool dotsys)
583{
584
585// Read gain-elevation ascii file data into a Table.
586
587  String formatString;
588  Table tbl = readAsciiTable(formatString, Table::Memory, filename, "", "", False);
589  scaleFromTable(in, tbl, method, xout, dotsys);
590}
591
592void STMath::scaleFromTable(Table& in,
593                            const Table& table,
594                            const std::string& method,
595                            const Vector<Float>& xout, bool dotsys)
596{
597
598  ROScalarColumn<Float> geElCol(table, "ELEVATION");
599  ROScalarColumn<Float> geFacCol(table, "FACTOR");
600  Vector<Float> xin = geElCol.getColumn();
601  Vector<Float> yin = geFacCol.getColumn();
602  Vector<Bool> maskin(xin.nelements(),True);
603
604  // Interpolate (and extrapolate) with desired method
605
606   //InterpolateArray1D<Double,Float>::InterpolationMethod method;
607   Int intmethod(stringToIMethod(method));
608
609   Vector<Float> yout;
610   Vector<Bool> maskout;
611   InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout,
612                                                xin, yin, maskin, intmethod,
613                                                True, True);
614
615   scaleByVector(in, Float(1.0)/yout, dotsys);
616}
617
618void STMath::scaleByVector( Table& in,
619                            const Vector< Float >& factor,
620                            bool dotsys )
621{
622  uInt nrow = in.nrow();
623  if ( factor.nelements() != nrow ) {
624    throw(AipsError("factors.nelements() != table.nelements()"));
625  }
626  ArrayColumn<Float> specCol(in, "SPECTRA");
627  ArrayColumn<uChar> flagCol(in, "FLAGTRA");
628  ArrayColumn<Float> tsysCol(in, "TSYS");
629  for (uInt i=0; i < nrow; ++i) {
630    MaskedArray<Float> ma  = maskedArray(specCol(i), flagCol(i));
631    ma *= factor[i];
632    specCol.put(i, ma.getArray());
633    flagCol.put(i, flagsFromMA(ma));
634    if ( dotsys ) {
635      Vector<Float> tsys;
636      tsysCol.get(i, tsys);
637      tsys *= factor[i];
638      specCol.put(i,tsys);
639    }
640  }
641}
642
643CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in,
644                                             float d, float etaap,
645                                             float jyperk )
646{
647  CountedPtr< Scantable > out = getScantable(in, false);
648  Table& tab = in->table();
649  Unit fluxUnit(tab.keywordSet().asString("FluxUnit"));
650  Unit K(String("K"));
651  Unit JY(String("Jy"));
652
653  bool tokelvin = true;
654  Double cfac = 1.0;
655
656  if ( fluxUnit == JY ) {
657    pushLog("Converting to K");
658    Quantum<Double> t(1.0,fluxUnit);
659    Quantum<Double> t2 = t.get(JY);
660    cfac = (t2 / t).getValue();               // value to Jy
661
662    tokelvin = true;
663    out->setFluxUnit("K");
664  } else if ( fluxUnit == K ) {
665    pushLog("Converting to Jy");
666    Quantum<Double> t(1.0,fluxUnit);
667    Quantum<Double> t2 = t.get(K);
668    cfac = (t2 / t).getValue();              // value to K
669
670    tokelvin = false;
671    out->setFluxUnit("Jy");
672  } else {
673    throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
674  }
675  // Make sure input values are converted to either Jy or K first...
676  Float factor = cfac;
677
678  // Select method
679  if (jyperk > 0.0) {
680    factor *= jyperk;
681    if ( tokelvin ) factor = 1.0 / jyperk;
682    ostringstream oss;
683    oss << "Jy/K = " << jyperk;
684    pushLog(String(oss));
685    Vector<Float> factors(tab.nrow(), factor);
686    scaleByVector(tab,factors, false);
687  } else if ( etaap > 0.0) {
688    Instrument inst =
689      SDAttr::convertInstrument(tab.keywordSet().asString("AntennaName"), True);
690    SDAttr sda;
691    if (d < 0) d = sda.diameter(inst);
692    Float jyPerk = SDAttr::findJyPerK(etaap, d);
693    ostringstream oss;
694    oss << "Jy/K = " << jyperk;
695    pushLog(String(oss));
696    factor *= jyperk;
697    if ( tokelvin ) {
698      factor = 1.0 / factor;
699    }
700    Vector<Float> factors(tab.nrow(), factor);
701    scaleByVector(tab, factors, False);
702  } else {
703
704    // OK now we must deal with automatic look up of values.
705    // We must also deal with the fact that the factors need
706    // to be computed per IF and may be different and may
707    // change per integration.
708
709    pushLog("Looking up conversion factors");
710    convertBrightnessUnits(out, tokelvin, cfac);
711  }
712
713  return out;
714}
715
716void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in,
717                                     bool tokelvin, float cfac )
718{
719  Table& table = in->table();
720  Instrument inst =
721    SDAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True);
722  TableIterator iter(table, "FREQ_ID");
723  STFrequencies stfreqs = in->frequencies();
724  SDAttr sdAtt;
725  while (!iter.pastEnd()) {
726    Table tab = iter.table();
727    ArrayColumn<Float> specCol(tab, "SPECTRA");
728    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
729    ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID");
730    MEpoch::ROScalarColumn timeCol(tab, "TIME");
731
732    uInt freqid; freqidCol.get(0, freqid);
733    Vector<Float> tmpspec; specCol.get(0, tmpspec);
734    // SDAttr.JyPerK has a Vector interface... change sometime.
735    Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements()));
736    for ( uInt i=0; i<tab.nrow(); ++i) {
737      Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0];
738      Float factor = cfac * jyperk;
739      if ( tokelvin ) factor = Float(1.0) / factor;
740      MaskedArray<Float> ma  = maskedArray(specCol(i), flagCol(i));
741      ma *= factor;
742      specCol.put(i, ma.getArray());
743      flagCol.put(i, flagsFromMA(ma));
744    }
745  }
746}
747
748CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
749                                         float tau )
750{
751  CountedPtr< Scantable > out = getScantable(in, false);
752  Table& tab = in->table();
753  ROScalarColumn<Float> elev(tab, "ELEVATION");
754  ArrayColumn<Float> specCol(tab, "SPECTRA");
755  ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
756  for ( uInt i=0; i<tab.nrow(); ++i) {
757    Float zdist = Float(C::pi_2) - elev(i);
758    Float factor = exp(tau)/cos(zdist);
759    MaskedArray<Float> ma  = maskedArray(specCol(i), flagCol(i));
760    ma *= factor;
761    specCol.put(i, ma.getArray());
762    flagCol.put(i, flagsFromMA(ma));
763  }
764  return out;
765}
766
767CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
768                                        const std::string& kernel, float width )
769{
770  CountedPtr< Scantable > out = getScantable(in, false);
771  Table& table = in->table();
772  VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
773  // same IFNO should have same no of channels
774  // this saves overhead
775  TableIterator iter(table, "IFNO");
776  while (!iter.pastEnd()) {
777    Table tab = iter.table();
778    ArrayColumn<Float> specCol(tab, "SPECTRA");
779    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
780    Vector<Float> tmpspec; specCol.get(0, tmpspec);
781    uInt nchan = tmpspec.nelements();
782    Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
783    Convolver<Float> conv(kvec, IPosition(1,nchan));
784    Vector<Float> spec;
785    Vector<uChar> flag;
786    for ( uInt i=0; i<tab.nrow(); ++i) {
787      specCol.get(i, spec);
788      flagCol.get(i, flag);
789      Vector<Bool> mask(flag.nelements());
790      convertArray(mask, flag);
791      Vector<Float> specout;
792      if ( type == VectorKernel::HANNING ) {
793        Vector<Bool> maskout;
794        mathutil::hanning(specout, maskout, spec , mask);
795        convertArray(flag, maskout);
796        flagCol.put(i, flag);
797      } else {
798        mathutil::replaceMaskByZero(specout, mask);
799        conv.linearConv(specout, spec);
800      }
801      specCol.put(i, specout);
802    }
803  }
804  return out;
805}
806
807CountedPtr< Scantable >
808  STMath::merge( const std::vector< CountedPtr < Scantable > >& in )
809{
810  if ( in.size() < 2 ) {
811    throw(AipsError("Need at least two scantables to perform a merge."));
812  }
813  std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
814  bool insitu = insitu_;
815  setInsitu(false);
816  CountedPtr< Scantable > out = getScantable(*it, false);
817  setInsitu(insitu);
818  Table& tout = out->table();
819  ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
820  ScalarColumn<uInt> scannocol(tout,"SCANNO"),focusidcol(tout,"FOCUS_ID");
821  uInt newscanno = max(scannocol.getColumn())+1;
822  ++it;
823  while ( it != in.end() ){
824    if ( ! (*it)->conformant(*out) ) {
825      // log message: "ignoring scantable i, as it isn't
826      // conformant with the other(s)"
827      cerr << "oh oh" << endl;
828      ++it;
829      continue;
830    }
831    out->appendToHistoryTable((*it)->history());
832    const Table& tab = (*it)->table();
833    TableIterator scanit(tab, "SCANNO");
834    while (!scanit.pastEnd()) {
835      TableIterator freqit(scanit.table(), "FREQ_ID");
836      while ( !freqit.pastEnd() ) {
837        Table thetab = freqit.table();
838        uInt nrow = tout.nrow();
839        //tout.addRow(thetab.nrow());
840        TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
841        ROTableRow row(thetab);
842        for ( uInt i=0; i<thetab.nrow(); ++i) {
843          uInt k = nrow+i;
844          scannocol.put(k, newscanno);
845          const TableRecord& rec = row.get(i);
846          Double rv,rp,inc;
847          (*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
848          uInt id;
849          id = out->frequencies().addEntry(rp, rv, inc);
850          freqidcol.put(k,id);
851          String name,fname;Double rf;
852          (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
853          id = out->molecules().addEntry(rf, name, fname);
854          molidcol.put(k, id);
855          Float frot,fang,ftan;
856          (*it)->focus().getEntry(frot, fang, ftan, rec.asuInt("FOCUS_ID"));
857          id = out->focus().addEntry(frot, fang, ftan);
858          focusidcol.put(k, id);
859        }
860        ++freqit;
861      }
862      ++newscanno;
863      ++scanit;
864    }
865    ++it;
866  }
867  return out;
868}
Note: See TracBrowser for help on using the repository browser.