source: trunk/src/STMath.cpp @ 934

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

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