source: trunk/src/STMath.cpp @ 805

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

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