source: trunk/src/STMath.cpp @ 926

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

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