source: trunk/src/STMath.cpp @ 961

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

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