source: trunk/src/STMath.cpp @ 858

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

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