source: branches/Release2.1.2/src/STMath.cpp @ 1310

Last change on this file since 1310 was 1310, checked in by mar637, 17 years ago

merge from trunk, to get binary operator changes
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 49.0 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 <casa/Arrays/MaskArrLogi.h>
18#include <casa/Arrays/MaskArrMath.h>
19#include <casa/Arrays/ArrayLogical.h>
20#include <casa/Arrays/ArrayMath.h>
21#include <casa/Arrays/Slice.h>
22#include <casa/Arrays/Slicer.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/TableParse.h>
30#include <tables/Tables/ReadAsciiTable.h>
31#include <tables/Tables/TableIter.h>
32#include <tables/Tables/TableCopy.h>
33#include <scimath/Mathematics/FFTServer.h>
34
35#include <lattices/Lattices/LatticeUtilities.h>
36
37#include <coordinates/Coordinates/SpectralCoordinate.h>
38#include <coordinates/Coordinates/CoordinateSystem.h>
39#include <coordinates/Coordinates/CoordinateUtil.h>
40#include <coordinates/Coordinates/FrequencyAligner.h>
41
42#include <scimath/Mathematics/VectorKernel.h>
43#include <scimath/Mathematics/Convolver.h>
44#include <scimath/Functionals/Polynomial.h>
45
46#include "MathUtils.h"
47#include "RowAccumulator.h"
48#include "STAttr.h"
49#include "STMath.h"
50
51using namespace casa;
52
53using namespace asap;
54
55STMath::STMath(bool insitu) :
56  insitu_(insitu)
57{
58}
59
60
61STMath::~STMath()
62{
63}
64
65CountedPtr<Scantable>
66STMath::average( const std::vector<CountedPtr<Scantable> >& in,
67                 const std::vector<bool>& mask,
68                 const std::string& weight,
69                 const std::string& avmode)
70{
71  if ( avmode == "SCAN" && in.size() != 1 )
72    throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
73                    "Use merge first."));
74  WeightType wtype = stringToWeight(weight);
75
76  // output
77  // clone as this is non insitu
78  bool insitu = insitu_;
79  setInsitu(false);
80  CountedPtr< Scantable > out = getScantable(in[0], true);
81  setInsitu(insitu);
82  std::vector<CountedPtr<Scantable> >::const_iterator stit = in.begin();
83  ++stit;
84  while ( stit != in.end() ) {
85    out->appendToHistoryTable((*stit)->history());
86    ++stit;
87  }
88
89  Table& tout = out->table();
90
91  /// @todo check if all scantables are conformant
92
93  ArrayColumn<Float> specColOut(tout,"SPECTRA");
94  ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
95  ArrayColumn<Float> tsysColOut(tout,"TSYS");
96  ScalarColumn<Double> mjdColOut(tout,"TIME");
97  ScalarColumn<Double> intColOut(tout,"INTERVAL");
98  ScalarColumn<uInt> cycColOut(tout,"CYCLENO");
99  ScalarColumn<uInt> scanColOut(tout,"SCANNO");
100
101  // set up the output table rows. These are based on the structure of the
102  // FIRST scantable in the vector
103  const Table& baset = in[0]->table();
104
105  Block<String> cols(3);
106  cols[0] = String("BEAMNO");
107  cols[1] = String("IFNO");
108  cols[2] = String("POLNO");
109  if ( avmode == "SOURCE" ) {
110    cols.resize(4);
111    cols[3] = String("SRCNAME");
112  }
113  if ( avmode == "SCAN"  && in.size() == 1) {
114    cols.resize(4);
115    cols[3] = String("SCANNO");
116  }
117  uInt outrowCount = 0;
118  TableIterator iter(baset, cols);
119  while (!iter.pastEnd()) {
120    Table subt = iter.table();
121    // copy the first row of this selection into the new table
122    tout.addRow();
123    TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
124    // re-index to 0
125    if ( avmode != "SCAN" && avmode != "SOURCE" ) {
126      scanColOut.put(outrowCount, uInt(0));
127    }
128    ++outrowCount;
129    ++iter;
130  }
131  RowAccumulator acc(wtype);
132  Vector<Bool> cmask(mask);
133  acc.setUserMask(cmask);
134  ROTableRow row(tout);
135  ROArrayColumn<Float> specCol, tsysCol;
136  ROArrayColumn<uChar> flagCol;
137  ROScalarColumn<Double> mjdCol, intCol;
138  ROScalarColumn<Int> scanIDCol;
139
140  for (uInt i=0; i < tout.nrow(); ++i) {
141    for ( int j=0; j < int(in.size()); ++j ) {
142      const Table& tin = in[j]->table();
143      const TableRecord& rec = row.get(i);
144      ROScalarColumn<Double> tmp(tin, "TIME");
145      Double td;tmp.get(0,td);
146      Table basesubt = tin(tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
147                       && tin.col("IFNO") == Int(rec.asuInt("IFNO"))
148                       && tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
149      Table subt;
150      if ( avmode == "SOURCE") {
151        subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME") );
152      } else if (avmode == "SCAN") {
153        subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
154      } else {
155        subt = basesubt;
156      }
157      specCol.attach(subt,"SPECTRA");
158      flagCol.attach(subt,"FLAGTRA");
159      tsysCol.attach(subt,"TSYS");
160      intCol.attach(subt,"INTERVAL");
161      mjdCol.attach(subt,"TIME");
162      Vector<Float> spec,tsys;
163      Vector<uChar> flag;
164      Double inter,time;
165      for (uInt k = 0; k < subt.nrow(); ++k ) {
166        flagCol.get(k, flag);
167        Vector<Bool> bflag(flag.shape());
168        convertArray(bflag, flag);
169        if ( allEQ(bflag, True) ) {
170          continue;//don't accumulate
171        }
172        specCol.get(k, spec);
173        tsysCol.get(k, tsys);
174        intCol.get(k, inter);
175        mjdCol.get(k, time);
176        // spectrum has to be added last to enable weighting by the other values
177        acc.add(spec, !bflag, tsys, inter, time);
178      }
179    }
180    //write out
181    specColOut.put(i, acc.getSpectrum());
182    const Vector<Bool>& msk = acc.getMask();
183    Vector<uChar> flg(msk.shape());
184    convertArray(flg, !msk);
185    flagColOut.put(i, flg);
186    tsysColOut.put(i, acc.getTsys());
187    intColOut.put(i, acc.getInterval());
188    mjdColOut.put(i, acc.getTime());
189    // we should only have one cycle now -> reset it to be 0
190    // frequency switched data has different CYCLENO for different IFNO
191    // which requires resetting this value
192    cycColOut.put(i, uInt(0));
193    acc.reset();
194  }
195  return out;
196}
197
198CountedPtr< Scantable >
199  STMath::averageChannel( const CountedPtr < Scantable > & in,
200                          const std::string & mode,
201                          const std::string& avmode )
202{
203  // clone as this is non insitu
204  bool insitu = insitu_;
205  setInsitu(false);
206  CountedPtr< Scantable > out = getScantable(in, true);
207  setInsitu(insitu);
208  Table& tout = out->table();
209  ArrayColumn<Float> specColOut(tout,"SPECTRA");
210  ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
211  ArrayColumn<Float> tsysColOut(tout,"TSYS");
212  ScalarColumn<uInt> scanColOut(tout,"SCANNO");
213  ScalarColumn<Double> intColOut(tout, "INTERVAL");
214  Table tmp = in->table().sort("BEAMNO");
215  Block<String> cols(3);
216  cols[0] = String("BEAMNO");
217  cols[1] = String("IFNO");
218  cols[2] = String("POLNO");
219  if ( avmode == "SCAN") {
220    cols.resize(4);
221    cols[3] = String("SCANNO");
222  }
223  uInt outrowCount = 0;
224  uChar userflag = 1 << 7;
225  TableIterator iter(tmp, cols);
226  while (!iter.pastEnd()) {
227    Table subt = iter.table();
228    ROArrayColumn<Float> specCol, tsysCol;
229    ROArrayColumn<uChar> flagCol;
230    ROScalarColumn<Double> intCol(subt, "INTERVAL");
231    specCol.attach(subt,"SPECTRA");
232    flagCol.attach(subt,"FLAGTRA");
233    tsysCol.attach(subt,"TSYS");
234    tout.addRow();
235    TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
236    if ( avmode != "SCAN") {
237      scanColOut.put(outrowCount, uInt(0));
238    }
239    Vector<Float> tmp;
240    specCol.get(0, tmp);
241    uInt nchan = tmp.nelements();
242    // have to do channel by channel here as MaskedArrMath
243    // doesn't have partialMedians
244    Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
245    Vector<Float> outspec(nchan);
246    Vector<uChar> outflag(nchan,0);
247    Vector<Float> outtsys(1);/// @fixme when tsys is channel based
248    for (uInt i=0; i<nchan; ++i) {
249      Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
250      MaskedArray<Float> ma = maskedArray(specs,flags);
251      outspec[i] = median(ma);
252      if ( allEQ(ma.getMask(), False) )
253        outflag[i] = userflag;// flag data
254    }
255    outtsys[0] = median(tsysCol.getColumn());
256    specColOut.put(outrowCount, outspec);
257    flagColOut.put(outrowCount, outflag);
258    tsysColOut.put(outrowCount, outtsys);
259    Double intsum = sum(intCol.getColumn());
260    intColOut.put(outrowCount, intsum);
261    ++outrowCount;
262    ++iter;
263  }
264  return out;
265}
266
267CountedPtr< Scantable > STMath::getScantable(const CountedPtr< Scantable >& in,
268                                             bool droprows)
269{
270  if (insitu_) return in;
271  else {
272    // clone
273    Scantable* tabp = new Scantable(*in, Bool(droprows));
274    return CountedPtr<Scantable>(tabp);
275  }
276}
277
278CountedPtr< Scantable > STMath::unaryOperate( const CountedPtr< Scantable >& in,
279                                              float val,
280                                              const std::string& mode,
281                                              bool tsys )
282{
283  CountedPtr< Scantable > out = getScantable(in, false);
284  Table& tab = out->table();
285  ArrayColumn<Float> specCol(tab,"SPECTRA");
286  ArrayColumn<Float> tsysCol(tab,"TSYS");
287  for (uInt i=0; i<tab.nrow(); ++i) {
288    Vector<Float> spec;
289    Vector<Float> ts;
290    specCol.get(i, spec);
291    tsysCol.get(i, ts);
292    if (mode == "MUL" || mode == "DIV") {
293      if (mode == "DIV") val = 1.0/val;
294      spec *= val;
295      specCol.put(i, spec);
296      if ( tsys ) {
297        ts *= val;
298        tsysCol.put(i, ts);
299      }
300    } else if ( mode == "ADD"  || mode == "SUB") {
301      if (mode == "SUB") val *= -1.0;
302      spec += val;
303      specCol.put(i, spec);
304      if ( tsys ) {
305        ts += val;
306        tsysCol.put(i, ts);
307      }
308    }
309  }
310  return out;
311}
312
313CountedPtr<Scantable> STMath::binaryOperate(const CountedPtr<Scantable>& left,
314                                            const CountedPtr<Scantable>& right,
315                                            const std::string& mode)
316{
317  bool insitu = insitu_;
318  if ( ! left->conformant(*right) ) {
319    throw(AipsError("'left' and 'right' scantables are not conformant."));
320  }
321  setInsitu(false);
322  CountedPtr< Scantable > out = getScantable(left, false);
323  setInsitu(insitu);
324  Table& tout = out->table();
325  Block<String> coln(5);
326  coln[0] = "SCANNO";  coln[1] = "CYCLENO";  coln[2] = "BEAMNO";
327  coln[3] = "IFNO";  coln[4] = "POLNO";
328  Table tmpl = tout.sort(coln);
329  Table tmpr = right->table().sort(coln);
330  ArrayColumn<Float> lspecCol(tmpl,"SPECTRA");
331  ROArrayColumn<Float> rspecCol(tmpr,"SPECTRA");
332  ArrayColumn<uChar> lflagCol(tmpl,"FLAGTRA");
333  ROArrayColumn<uChar> rflagCol(tmpr,"FLAGTRA");
334
335  for (uInt i=0; i<tout.nrow(); ++i) {
336    Vector<Float> lspecvec, rspecvec;
337    Vector<uChar> lflagvec, rflagvec;
338    lspecvec = lspecCol(i);    rspecvec = rspecCol(i);
339    lflagvec = lflagCol(i);    rflagvec = rflagCol(i);
340    MaskedArray<Float> mleft = maskedArray(lspecvec, lflagvec);
341    MaskedArray<Float> mright = maskedArray(rspecvec, rflagvec);
342    if (mode == "ADD") {
343      mleft += mright;
344    } else if ( mode == "SUB") {
345      mleft -= mright;
346    } else if ( mode == "MUL") {
347      mleft *= mright;
348    } else if ( mode == "DIV") {
349      mleft /= mright;
350    } else {
351      throw(AipsError("Illegal binary operator"));
352    }
353    lspecCol.put(i, mleft.getArray());
354  }
355  return out;
356}
357
358
359
360MaskedArray<Float> STMath::maskedArray( const Vector<Float>& s,
361                                        const Vector<uChar>& f)
362{
363  Vector<Bool> mask;
364  mask.resize(f.shape());
365  convertArray(mask, f);
366  return MaskedArray<Float>(s,!mask);
367}
368
369Vector<uChar> STMath::flagsFromMA(const MaskedArray<Float>& ma)
370{
371  const Vector<Bool>& m = ma.getMask();
372  Vector<uChar> flags(m.shape());
373  convertArray(flags, !m);
374  return flags;
375}
376
377CountedPtr< Scantable > STMath::autoQuotient( const CountedPtr< Scantable >& in,
378                                              const std::string & mode,
379                                              bool preserve )
380{
381  /// @todo make other modes available
382  /// modes should be "nearest", "pair"
383  // make this operation non insitu
384  const Table& tin = in->table();
385  Table ons = tin(tin.col("SRCTYPE") == Int(0));
386  Table offs = tin(tin.col("SRCTYPE") == Int(1));
387  if ( offs.nrow() == 0 )
388    throw(AipsError("No 'off' scans present."));
389  // put all "on" scans into output table
390
391  bool insitu = insitu_;
392  setInsitu(false);
393  CountedPtr< Scantable > out = getScantable(in, true);
394  setInsitu(insitu);
395  Table& tout = out->table();
396
397  TableCopy::copyRows(tout, ons);
398  TableRow row(tout);
399  ROScalarColumn<Double> offtimeCol(offs, "TIME");
400  ArrayColumn<Float> outspecCol(tout, "SPECTRA");
401  ROArrayColumn<Float> outtsysCol(tout, "TSYS");
402  ArrayColumn<uChar> outflagCol(tout, "FLAGTRA");
403  for (uInt i=0; i < tout.nrow(); ++i) {
404    const TableRecord& rec = row.get(i);
405    Double ontime = rec.asDouble("TIME");
406    ROScalarColumn<Double> offtimeCol(offs, "TIME");
407    Double mindeltat = min(abs(offtimeCol.getColumn() - ontime));
408    // Timestamp may vary within a cycle ???!!!
409    // increase this by 0.5 sec in case of rounding errors...
410    // There might be a better way to do this.
411    mindeltat += 0.5;
412    Table sel = offs( abs(offs.col("TIME")-ontime) <= (mindeltat+0.5)
413                       && offs.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
414                       && offs.col("IFNO") == Int(rec.asuInt("IFNO"))
415                       && offs.col("POLNO") == Int(rec.asuInt("POLNO")) );
416
417    if ( sel.nrow() < 1 )  {
418      throw(AipsError("No closest in time found... This could be a rounding "
419                      "issue. Try quotient instead."));
420    }
421    TableRow offrow(sel);
422    const TableRecord& offrec = offrow.get(0);//should only be one row
423    RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
424    RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
425    RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
426    /// @fixme this assumes tsys is a scalar not vector
427    Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
428    Vector<Float> specon, tsyson;
429    outtsysCol.get(i, tsyson);
430    outspecCol.get(i, specon);
431    Vector<uChar> flagon;
432    outflagCol.get(i, flagon);
433    MaskedArray<Float> mon = maskedArray(specon, flagon);
434    MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
435    MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
436    if (preserve) {
437      quot -= tsysoffscalar;
438    } else {
439      quot -= tsyson[0];
440    }
441    outspecCol.put(i, quot.getArray());
442    outflagCol.put(i, flagsFromMA(quot));
443  }
444  // renumber scanno
445  TableIterator it(tout, "SCANNO");
446  uInt i = 0;
447  while ( !it.pastEnd() ) {
448    Table t = it.table();
449    TableVector<uInt> vec(t, "SCANNO");
450    vec = i;
451    ++i;
452    ++it;
453  }
454  return out;
455}
456
457
458CountedPtr< Scantable > STMath::quotient( const CountedPtr< Scantable > & on,
459                                          const CountedPtr< Scantable > & off,
460                                          bool preserve )
461{
462  bool insitu = insitu_;
463  if ( ! on->conformant(*off) ) {
464    throw(AipsError("'on' and 'off' scantables are not conformant."));
465  }
466  setInsitu(false);
467  CountedPtr< Scantable > out = getScantable(on, false);
468  setInsitu(insitu);
469  Table& tout = out->table();
470  const Table& toff = off->table();
471  TableIterator sit(tout, "SCANNO");
472  TableIterator s2it(toff, "SCANNO");
473  while ( !sit.pastEnd() ) {
474    Table ton = sit.table();
475    TableRow row(ton);
476    Table t = s2it.table();
477    ArrayColumn<Float> outspecCol(ton, "SPECTRA");
478    ROArrayColumn<Float> outtsysCol(ton, "TSYS");
479    ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
480    for (uInt i=0; i < ton.nrow(); ++i) {
481      const TableRecord& rec = row.get(i);
482      Table offsel = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
483                          && t.col("IFNO") == Int(rec.asuInt("IFNO"))
484                          && t.col("POLNO") == Int(rec.asuInt("POLNO")) );
485      if ( offsel.nrow() == 0 )
486        throw AipsError("STMath::quotient: no matching off");
487      TableRow offrow(offsel);
488      const TableRecord& offrec = offrow.get(0);//should be ncycles - take first
489      RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
490      RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
491      RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
492      Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
493      Vector<Float> specon, tsyson;
494      outtsysCol.get(i, tsyson);
495      outspecCol.get(i, specon);
496      Vector<uChar> flagon;
497      outflagCol.get(i, flagon);
498      MaskedArray<Float> mon = maskedArray(specon, flagon);
499      MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
500      MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
501      if (preserve) {
502        quot -= tsysoffscalar;
503      } else {
504        quot -= tsyson[0];
505      }
506      outspecCol.put(i, quot.getArray());
507      outflagCol.put(i, flagsFromMA(quot));
508    }
509    ++sit;
510    ++s2it;
511    // take the first off for each on scan which doesn't have a
512    // matching off scan
513    // non <= noff:  matching pairs, non > noff matching pairs then first off
514    if ( s2it.pastEnd() ) s2it.reset();
515  }
516  return out;
517}
518
519
520CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
521{
522  // make copy or reference
523  CountedPtr< Scantable > out = getScantable(in, false);
524  Table& tout = out->table();
525  Block<String> cols(4);
526  cols[0] = String("SCANNO");
527  cols[1] = String("CYCLENO");
528  cols[2] = String("BEAMNO");
529  cols[3] = String("POLNO");
530  TableIterator iter(tout, cols);
531  while (!iter.pastEnd()) {
532    Table subt = iter.table();
533    // this should leave us with two rows for the two IFs....if not ignore
534    if (subt.nrow() != 2 ) {
535      continue;
536    }
537    ArrayColumn<Float> specCol(subt, "SPECTRA");
538    ArrayColumn<Float> tsysCol(subt, "TSYS");
539    ArrayColumn<uChar> flagCol(subt, "FLAGTRA");
540    Vector<Float> onspec,offspec, ontsys, offtsys;
541    Vector<uChar> onflag, offflag;
542    tsysCol.get(0, ontsys);   tsysCol.get(1, offtsys);
543    specCol.get(0, onspec);   specCol.get(1, offspec);
544    flagCol.get(0, onflag);   flagCol.get(1, offflag);
545    MaskedArray<Float> on  = maskedArray(onspec, onflag);
546    MaskedArray<Float> off = maskedArray(offspec, offflag);
547    MaskedArray<Float> oncopy = on.copy();
548
549    on /= off; on -= 1.0f;
550    on *= ontsys[0];
551    off /= oncopy; off -= 1.0f;
552    off *= offtsys[0];
553    specCol.put(0, on.getArray());
554    const Vector<Bool>& m0 = on.getMask();
555    Vector<uChar> flags0(m0.shape());
556    convertArray(flags0, !m0);
557    flagCol.put(0, flags0);
558
559    specCol.put(1, off.getArray());
560    const Vector<Bool>& m1 = off.getMask();
561    Vector<uChar> flags1(m1.shape());
562    convertArray(flags1, !m1);
563    flagCol.put(1, flags1);
564    ++iter;
565  }
566
567  return out;
568}
569
570std::vector< float > STMath::statistic( const CountedPtr< Scantable > & in,
571                                        const std::vector< bool > & mask,
572                                        const std::string& which )
573{
574
575  Vector<Bool> m(mask);
576  const Table& tab = in->table();
577  ROArrayColumn<Float> specCol(tab, "SPECTRA");
578  ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
579  std::vector<float> out;
580  for (uInt i=0; i < tab.nrow(); ++i ) {
581    Vector<Float> spec; specCol.get(i, spec);
582    Vector<uChar> flag; flagCol.get(i, flag);
583    MaskedArray<Float> ma  = maskedArray(spec, flag);
584    float outstat = 0.0;
585    if ( spec.nelements() == m.nelements() ) {
586      outstat = mathutil::statistics(which, ma(m));
587    } else {
588      outstat = mathutil::statistics(which, ma);
589    }
590    out.push_back(outstat);
591  }
592  return out;
593}
594
595CountedPtr< Scantable > STMath::bin( const CountedPtr< Scantable > & in,
596                                     int width )
597{
598  if ( !in->getSelection().empty() ) throw(AipsError("Can't bin subset of the data."));
599  CountedPtr< Scantable > out = getScantable(in, false);
600  Table& tout = out->table();
601  out->frequencies().rescale(width, "BIN");
602  ArrayColumn<Float> specCol(tout, "SPECTRA");
603  ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
604  for (uInt i=0; i < tout.nrow(); ++i ) {
605    MaskedArray<Float> main  = maskedArray(specCol(i), flagCol(i));
606    MaskedArray<Float> maout;
607    LatticeUtilities::bin(maout, main, 0, Int(width));
608    /// @todo implement channel based tsys binning
609    specCol.put(i, maout.getArray());
610    flagCol.put(i, flagsFromMA(maout));
611    // take only the first binned spectrum's length for the deprecated
612    // global header item nChan
613    if (i==0) tout.rwKeywordSet().define(String("nChan"),
614                                       Int(maout.getArray().nelements()));
615  }
616  return out;
617}
618
619CountedPtr< Scantable > STMath::resample( const CountedPtr< Scantable >& in,
620                                          const std::string& method,
621                                          float width )
622//
623// Should add the possibility of width being specified in km/s. This means
624// that for each freqID (SpectralCoordinate) we will need to convert to an
625// average channel width (say at the reference pixel).  Then we would need
626// to be careful to make sure each spectrum (of different freqID)
627// is the same length.
628//
629{
630  //InterpolateArray1D<Double,Float>::InterpolationMethod interp;
631  Int interpMethod(stringToIMethod(method));
632
633  CountedPtr< Scantable > out = getScantable(in, false);
634  Table& tout = out->table();
635
636// Resample SpectralCoordinates (one per freqID)
637  out->frequencies().rescale(width, "RESAMPLE");
638  TableIterator iter(tout, "IFNO");
639  TableRow row(tout);
640  while ( !iter.pastEnd() ) {
641    Table tab = iter.table();
642    ArrayColumn<Float> specCol(tab, "SPECTRA");
643    //ArrayColumn<Float> tsysCol(tout, "TSYS");
644    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
645    Vector<Float> spec;
646    Vector<uChar> flag;
647    specCol.get(0,spec); // the number of channels should be constant per IF
648    uInt nChanIn = spec.nelements();
649    Vector<Float> xIn(nChanIn); indgen(xIn);
650    Int fac =  Int(nChanIn/width);
651    Vector<Float> xOut(fac+10); // 10 to be safe - resize later
652    uInt k = 0;
653    Float x = 0.0;
654    while (x < Float(nChanIn) ) {
655      xOut(k) = x;
656      k++;
657      x += width;
658    }
659    uInt nChanOut = k;
660    xOut.resize(nChanOut, True);
661    // process all rows for this IFNO
662    Vector<Float> specOut;
663    Vector<Bool> maskOut;
664    Vector<uChar> flagOut;
665    for (uInt i=0; i < tab.nrow(); ++i) {
666      specCol.get(i, spec);
667      flagCol.get(i, flag);
668      Vector<Bool> mask(flag.nelements());
669      convertArray(mask, flag);
670
671      IPosition shapeIn(spec.shape());
672      //sh.nchan = nChanOut;
673      InterpolateArray1D<Float,Float>::interpolate(specOut, maskOut, xOut,
674                                                   xIn, spec, mask,
675                                                   interpMethod, True, True);
676      /// @todo do the same for channel based Tsys
677      flagOut.resize(maskOut.nelements());
678      convertArray(flagOut, maskOut);
679      specCol.put(i, specOut);
680      flagCol.put(i, flagOut);
681    }
682    ++iter;
683  }
684
685  return out;
686}
687
688STMath::imethod STMath::stringToIMethod(const std::string& in)
689{
690  static STMath::imap lookup;
691
692  // initialize the lookup table if necessary
693  if ( lookup.empty() ) {
694    lookup["nearest"]   = InterpolateArray1D<Double,Float>::nearestNeighbour;
695    lookup["linear"] = InterpolateArray1D<Double,Float>::linear;
696    lookup["cubic"]  = InterpolateArray1D<Double,Float>::cubic;
697    lookup["spline"]  = InterpolateArray1D<Double,Float>::spline;
698  }
699
700  STMath::imap::const_iterator iter = lookup.find(in);
701
702  if ( lookup.end() == iter ) {
703    std::string message = in;
704    message += " is not a valid interpolation mode";
705    throw(AipsError(message));
706  }
707  return iter->second;
708}
709
710WeightType STMath::stringToWeight(const std::string& in)
711{
712  static std::map<std::string, WeightType> lookup;
713
714  // initialize the lookup table if necessary
715  if ( lookup.empty() ) {
716    lookup["NONE"]   = asap::NONE;
717    lookup["TINT"] = asap::TINT;
718    lookup["TINTSYS"]  = asap::TINTSYS;
719    lookup["TSYS"]  = asap::TSYS;
720    lookup["VAR"]  = asap::VAR;
721  }
722
723  std::map<std::string, WeightType>::const_iterator iter = lookup.find(in);
724
725  if ( lookup.end() == iter ) {
726    std::string message = in;
727    message += " is not a valid weighting mode";
728    throw(AipsError(message));
729  }
730  return iter->second;
731}
732
733CountedPtr< Scantable > STMath::gainElevation( const CountedPtr< Scantable >& in,
734                                               const vector< float > & coeff,
735                                               const std::string & filename,
736                                               const std::string& method)
737{
738  // Get elevation data from Scantable and convert to degrees
739  CountedPtr< Scantable > out = getScantable(in, false);
740  Table& tab = out->table();
741  ROScalarColumn<Float> elev(tab, "ELEVATION");
742  Vector<Float> x = elev.getColumn();
743  x *= Float(180 / C::pi);                        // Degrees
744
745  Vector<Float> coeffs(coeff);
746  const uInt nc = coeffs.nelements();
747  if ( filename.length() > 0 && nc > 0 ) {
748    throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
749  }
750
751  // Correct
752  if ( nc > 0 || filename.length() == 0 ) {
753    // Find instrument
754    Bool throwit = True;
755    Instrument inst =
756      STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
757                                throwit);
758
759    // Set polynomial
760    Polynomial<Float>* ppoly = 0;
761    Vector<Float> coeff;
762    String msg;
763    if ( nc > 0 ) {
764      ppoly = new Polynomial<Float>(nc);
765      coeff = coeffs;
766      msg = String("user");
767    } else {
768      STAttr sdAttr;
769      coeff = sdAttr.gainElevationPoly(inst);
770      ppoly = new Polynomial<Float>(3);
771      msg = String("built in");
772    }
773
774    if ( coeff.nelements() > 0 ) {
775      ppoly->setCoefficients(coeff);
776    } else {
777      delete ppoly;
778      throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
779    }
780    ostringstream oss;
781    oss << "Making polynomial correction with " << msg << " coefficients:" << endl;
782    oss << "   " <<  coeff;
783    pushLog(String(oss));
784    const uInt nrow = tab.nrow();
785    Vector<Float> factor(nrow);
786    for ( uInt i=0; i < nrow; ++i ) {
787      factor[i] = 1.0 / (*ppoly)(x[i]);
788    }
789    delete ppoly;
790    scaleByVector(tab, factor, true);
791
792  } else {
793    // Read and correct
794    pushLog("Making correction from ascii Table");
795    scaleFromAsciiTable(tab, filename, method, x, true);
796  }
797  return out;
798}
799
800void STMath::scaleFromAsciiTable(Table& in, const std::string& filename,
801                                 const std::string& method,
802                                 const Vector<Float>& xout, bool dotsys)
803{
804
805// Read gain-elevation ascii file data into a Table.
806
807  String formatString;
808  Table tbl = readAsciiTable(formatString, Table::Memory, filename, "", "", False);
809  scaleFromTable(in, tbl, method, xout, dotsys);
810}
811
812void STMath::scaleFromTable(Table& in,
813                            const Table& table,
814                            const std::string& method,
815                            const Vector<Float>& xout, bool dotsys)
816{
817
818  ROScalarColumn<Float> geElCol(table, "ELEVATION");
819  ROScalarColumn<Float> geFacCol(table, "FACTOR");
820  Vector<Float> xin = geElCol.getColumn();
821  Vector<Float> yin = geFacCol.getColumn();
822  Vector<Bool> maskin(xin.nelements(),True);
823
824  // Interpolate (and extrapolate) with desired method
825
826  InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
827
828   Vector<Float> yout;
829   Vector<Bool> maskout;
830   InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout,
831                                                xin, yin, maskin, interp,
832                                                True, True);
833
834   scaleByVector(in, Float(1.0)/yout, dotsys);
835}
836
837void STMath::scaleByVector( Table& in,
838                            const Vector< Float >& factor,
839                            bool dotsys )
840{
841  uInt nrow = in.nrow();
842  if ( factor.nelements() != nrow ) {
843    throw(AipsError("factors.nelements() != table.nelements()"));
844  }
845  ArrayColumn<Float> specCol(in, "SPECTRA");
846  ArrayColumn<uChar> flagCol(in, "FLAGTRA");
847  ArrayColumn<Float> tsysCol(in, "TSYS");
848  for (uInt i=0; i < nrow; ++i) {
849    MaskedArray<Float> ma  = maskedArray(specCol(i), flagCol(i));
850    ma *= factor[i];
851    specCol.put(i, ma.getArray());
852    flagCol.put(i, flagsFromMA(ma));
853    if ( dotsys ) {
854      Vector<Float> tsys = tsysCol(i);
855      tsys *= factor[i];
856      tsysCol.put(i,tsys);
857    }
858  }
859}
860
861CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in,
862                                             float d, float etaap,
863                                             float jyperk )
864{
865  CountedPtr< Scantable > out = getScantable(in, false);
866  Table& tab = in->table();
867  Unit fluxUnit(tab.keywordSet().asString("FluxUnit"));
868  Unit K(String("K"));
869  Unit JY(String("Jy"));
870
871  bool tokelvin = true;
872  Double cfac = 1.0;
873
874  if ( fluxUnit == JY ) {
875    pushLog("Converting to K");
876    Quantum<Double> t(1.0,fluxUnit);
877    Quantum<Double> t2 = t.get(JY);
878    cfac = (t2 / t).getValue();               // value to Jy
879
880    tokelvin = true;
881    out->setFluxUnit("K");
882  } else if ( fluxUnit == K ) {
883    pushLog("Converting to Jy");
884    Quantum<Double> t(1.0,fluxUnit);
885    Quantum<Double> t2 = t.get(K);
886    cfac = (t2 / t).getValue();              // value to K
887
888    tokelvin = false;
889    out->setFluxUnit("Jy");
890  } else {
891    throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
892  }
893  // Make sure input values are converted to either Jy or K first...
894  Float factor = cfac;
895
896  // Select method
897  if (jyperk > 0.0) {
898    factor *= jyperk;
899    if ( tokelvin ) factor = 1.0 / jyperk;
900    ostringstream oss;
901    oss << "Jy/K = " << jyperk;
902    pushLog(String(oss));
903    Vector<Float> factors(tab.nrow(), factor);
904    scaleByVector(tab,factors, false);
905  } else if ( etaap > 0.0) {
906    Instrument inst =
907      STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"), True);
908    STAttr sda;
909    if (d < 0) d = sda.diameter(inst);
910    jyperk = STAttr::findJyPerK(etaap, d);
911    ostringstream oss;
912    oss << "Jy/K = " << jyperk;
913    pushLog(String(oss));
914    factor *= jyperk;
915    if ( tokelvin ) {
916      factor = 1.0 / factor;
917    }
918    Vector<Float> factors(tab.nrow(), factor);
919    scaleByVector(tab, factors, False);
920  } else {
921
922    // OK now we must deal with automatic look up of values.
923    // We must also deal with the fact that the factors need
924    // to be computed per IF and may be different and may
925    // change per integration.
926
927    pushLog("Looking up conversion factors");
928    convertBrightnessUnits(out, tokelvin, cfac);
929  }
930
931  return out;
932}
933
934void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in,
935                                     bool tokelvin, float cfac )
936{
937  Table& table = in->table();
938  Instrument inst =
939    STAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True);
940  TableIterator iter(table, "FREQ_ID");
941  STFrequencies stfreqs = in->frequencies();
942  STAttr sdAtt;
943  while (!iter.pastEnd()) {
944    Table tab = iter.table();
945    ArrayColumn<Float> specCol(tab, "SPECTRA");
946    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
947    ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID");
948    MEpoch::ROScalarColumn timeCol(tab, "TIME");
949
950    uInt freqid; freqidCol.get(0, freqid);
951    Vector<Float> tmpspec; specCol.get(0, tmpspec);
952    // STAttr.JyPerK has a Vector interface... change sometime.
953    Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements()));
954    for ( uInt i=0; i<tab.nrow(); ++i) {
955      Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0];
956      Float factor = cfac * jyperk;
957      if ( tokelvin ) factor = Float(1.0) / factor;
958      MaskedArray<Float> ma  = maskedArray(specCol(i), flagCol(i));
959      ma *= factor;
960      specCol.put(i, ma.getArray());
961      flagCol.put(i, flagsFromMA(ma));
962    }
963  ++iter;
964  }
965}
966
967CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
968                                         float tau )
969{
970  CountedPtr< Scantable > out = getScantable(in, false);
971
972  Table tab = out->table();
973  ROScalarColumn<Float> elev(tab, "ELEVATION");
974  ArrayColumn<Float> specCol(tab, "SPECTRA");
975  ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
976  for ( uInt i=0; i<tab.nrow(); ++i) {
977    Float zdist = Float(C::pi_2) - elev(i);
978    Float factor = exp(tau)/cos(zdist);
979    MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
980    ma *= factor;
981    specCol.put(i, ma.getArray());
982    flagCol.put(i, flagsFromMA(ma));
983  }
984  return out;
985}
986
987CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
988                                        const std::string& kernel, float width )
989{
990  CountedPtr< Scantable > out = getScantable(in, false);
991  Table& table = out->table();
992  VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
993  // same IFNO should have same no of channels
994  // this saves overhead
995  TableIterator iter(table, "IFNO");
996  while (!iter.pastEnd()) {
997    Table tab = iter.table();
998    ArrayColumn<Float> specCol(tab, "SPECTRA");
999    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1000    Vector<Float> tmpspec; specCol.get(0, tmpspec);
1001    uInt nchan = tmpspec.nelements();
1002    Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
1003    Convolver<Float> conv(kvec, IPosition(1,nchan));
1004    Vector<Float> spec;
1005    Vector<uChar> flag;
1006    for ( uInt i=0; i<tab.nrow(); ++i) {
1007      specCol.get(i, spec);
1008      flagCol.get(i, flag);
1009      Vector<Bool> mask(flag.nelements());
1010      convertArray(mask, flag);
1011      Vector<Float> specout;
1012      if ( type == VectorKernel::HANNING ) {
1013        Vector<Bool> maskout;
1014        mathutil::hanning(specout, maskout, spec , !mask);
1015        convertArray(flag, !maskout);
1016        flagCol.put(i, flag);
1017        specCol.put(i, specout);
1018     } else {
1019        mathutil::replaceMaskByZero(specout, mask);
1020        conv.linearConv(specout, spec);
1021        specCol.put(i, specout);
1022      }
1023    }
1024    ++iter;
1025  }
1026  return out;
1027}
1028
1029CountedPtr< Scantable >
1030  STMath::merge( const std::vector< CountedPtr < Scantable > >& in )
1031{
1032  if ( in.size() < 2 ) {
1033    throw(AipsError("Need at least two scantables to perform a merge."));
1034  }
1035  std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
1036  bool insitu = insitu_;
1037  setInsitu(false);
1038  CountedPtr< Scantable > out = getScantable(*it, false);
1039  setInsitu(insitu);
1040  Table& tout = out->table();
1041  ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
1042  ScalarColumn<uInt> scannocol(tout,"SCANNO"), focusidcol(tout,"FOCUS_ID");
1043  // Renumber SCANNO to be 0-based
1044  Vector<uInt> scannos = scannocol.getColumn();
1045  uInt offset = min(scannos);
1046  scannos -= offset;
1047  scannocol.putColumn(scannos);
1048  uInt newscanno = max(scannos)+1;
1049  ++it;
1050  while ( it != in.end() ){
1051    if ( ! (*it)->conformant(*out) ) {
1052      // log message: "ignoring scantable i, as it isn't
1053      // conformant with the other(s)"
1054      cerr << "oh oh" << endl;
1055      ++it;
1056      continue;
1057    }
1058    out->appendToHistoryTable((*it)->history());
1059    const Table& tab = (*it)->table();
1060    TableIterator scanit(tab, "SCANNO");
1061    while (!scanit.pastEnd()) {
1062      TableIterator freqit(scanit.table(), "FREQ_ID");
1063      while ( !freqit.pastEnd() ) {
1064        Table thetab = freqit.table();
1065        uInt nrow = tout.nrow();
1066        //tout.addRow(thetab.nrow());
1067        TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
1068        ROTableRow row(thetab);
1069        for ( uInt i=0; i<thetab.nrow(); ++i) {
1070          uInt k = nrow+i;
1071          scannocol.put(k, newscanno);
1072          const TableRecord& rec = row.get(i);
1073          Double rv,rp,inc;
1074          (*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
1075          uInt id;
1076          id = out->frequencies().addEntry(rp, rv, inc);
1077          freqidcol.put(k,id);
1078          String name,fname;Double rf;
1079          (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
1080          id = out->molecules().addEntry(rf, name, fname);
1081          molidcol.put(k, id);
1082          Float frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
1083          (*it)->focus().getEntry(fax, ftan, frot, fhand,
1084                                  fmount,fuser, fxy, fxyp,
1085                                  rec.asuInt("FOCUS_ID"));
1086          id = out->focus().addEntry(fax, ftan, frot, fhand,
1087                                     fmount,fuser, fxy, fxyp);
1088          focusidcol.put(k, id);
1089        }
1090        ++freqit;
1091      }
1092      ++newscanno;
1093      ++scanit;
1094    }
1095    ++it;
1096  }
1097  return out;
1098}
1099
1100CountedPtr< Scantable >
1101  STMath::invertPhase( const CountedPtr < Scantable >& in )
1102{
1103  return applyToPol(in, &STPol::invertPhase, Float(0.0));
1104}
1105
1106CountedPtr< Scantable >
1107  STMath::rotateXYPhase( const CountedPtr < Scantable >& in, float phase )
1108{
1109   return applyToPol(in, &STPol::rotatePhase, Float(phase));
1110}
1111
1112CountedPtr< Scantable >
1113  STMath::rotateLinPolPhase( const CountedPtr < Scantable >& in, float phase )
1114{
1115  return applyToPol(in, &STPol::rotateLinPolPhase, Float(phase));
1116}
1117
1118CountedPtr< Scantable > STMath::applyToPol( const CountedPtr<Scantable>& in,
1119                                             STPol::polOperation fptr,
1120                                             Float phase )
1121{
1122  CountedPtr< Scantable > out = getScantable(in, false);
1123  Table& tout = out->table();
1124  Block<String> cols(4);
1125  cols[0] = String("SCANNO");
1126  cols[1] = String("BEAMNO");
1127  cols[2] = String("IFNO");
1128  cols[3] = String("CYCLENO");
1129  TableIterator iter(tout, cols);
1130  STPol* stpol = STPol::getPolClass(out->factories_, out->getPolType() );
1131  while (!iter.pastEnd()) {
1132    Table t = iter.table();
1133    ArrayColumn<Float> speccol(t, "SPECTRA");
1134    ScalarColumn<uInt> focidcol(t, "FOCUS_ID");
1135    ScalarColumn<Float> parancol(t, "PARANGLE");
1136    Matrix<Float> pols = speccol.getColumn();
1137    try {
1138      stpol->setSpectra(pols);
1139      Float fang,fhand,parang;
1140      fang = in->focusTable_.getTotalFeedAngle(focidcol(0));
1141      fhand = in->focusTable_.getFeedHand(focidcol(0));
1142      parang = parancol(0);
1143      /// @todo re-enable this
1144      // disable total feed angle to support paralactifying Caswell style
1145      stpol->setPhaseCorrections(parang, -parang, fhand);
1146      (stpol->*fptr)(phase);
1147      speccol.putColumn(stpol->getSpectra());
1148      Matrix<Float> tmp = stpol->getSpectra();
1149    } catch (AipsError& e) {
1150      delete stpol;stpol=0;
1151      throw(e);
1152    }
1153    ++iter;
1154  }
1155  delete stpol;stpol=0;
1156  return out;
1157}
1158
1159CountedPtr< Scantable >
1160  STMath::swapPolarisations( const CountedPtr< Scantable > & in )
1161{
1162  CountedPtr< Scantable > out = getScantable(in, false);
1163  Table& tout = out->table();
1164  Table t0 = tout(tout.col("POLNO") == 0);
1165  Table t1 = tout(tout.col("POLNO") == 1);
1166  if ( t0.nrow() != t1.nrow() )
1167    throw(AipsError("Inconsistent number of polarisations"));
1168  ArrayColumn<Float> speccol0(t0, "SPECTRA");
1169  ArrayColumn<uChar> flagcol0(t0, "FLAGTRA");
1170  ArrayColumn<Float> speccol1(t1, "SPECTRA");
1171  ArrayColumn<uChar> flagcol1(t1, "FLAGTRA");
1172  Matrix<Float> s0 = speccol0.getColumn();
1173  Matrix<uChar> f0 = flagcol0.getColumn();
1174  speccol0.putColumn(speccol1.getColumn());
1175  flagcol0.putColumn(flagcol1.getColumn());
1176  speccol1.putColumn(s0);
1177  flagcol1.putColumn(f0);
1178  return out;
1179}
1180
1181CountedPtr< Scantable >
1182  STMath::averagePolarisations( const CountedPtr< Scantable > & in,
1183                                const std::vector<bool>& mask,
1184                                const std::string& weight )
1185{
1186  if (in->npol() < 2 )
1187    throw(AipsError("averagePolarisations can only be applied to two or more"
1188                    "polarisations"));
1189  bool insitu = insitu_;
1190  setInsitu(false);
1191  CountedPtr< Scantable > pols = getScantable(in, true);
1192  setInsitu(insitu);
1193  Table& tout = pols->table();
1194  std::string taql = "SELECT FROM $1 WHERE POLNO IN [0,1]";
1195  Table tab = tableCommand(taql, in->table());
1196  if (tab.nrow() == 0 )
1197    throw(AipsError("Could not find  any rows with POLNO==0 and POLNO==1"));
1198  TableCopy::copyRows(tout, tab);
1199  TableVector<uInt> vec(tout, "POLNO");
1200  vec = 0;
1201  pols->table_.rwKeywordSet().define("nPol", Int(1));
1202  pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
1203  std::vector<CountedPtr<Scantable> > vpols;
1204  vpols.push_back(pols);
1205  CountedPtr< Scantable > out = average(vpols, mask, weight, "SCAN");
1206  return out;
1207}
1208
1209CountedPtr< Scantable >
1210  STMath::averageBeams( const CountedPtr< Scantable > & in,
1211                        const std::vector<bool>& mask,
1212                        const std::string& weight )
1213{
1214  bool insitu = insitu_;
1215  setInsitu(false);
1216  CountedPtr< Scantable > beams = getScantable(in, false);
1217  setInsitu(insitu);
1218  Table& tout = beams->table();
1219  // give all rows the same BEAMNO
1220  TableVector<uInt> vec(tout, "BEAMNO");
1221  vec = 0;
1222  beams->table_.rwKeywordSet().define("nBeam", Int(1));
1223  std::vector<CountedPtr<Scantable> > vbeams;
1224  vbeams.push_back(beams);
1225  CountedPtr< Scantable > out = average(vbeams, mask, weight, "SCAN");
1226  return out;
1227}
1228
1229
1230CountedPtr< Scantable >
1231  asap::STMath::frequencyAlign( const CountedPtr< Scantable > & in,
1232                                const std::string & refTime,
1233                                const std::string & method)
1234{
1235  // clone as this is not working insitu
1236  bool insitu = insitu_;
1237  setInsitu(false);
1238  CountedPtr< Scantable > out = getScantable(in, false);
1239  setInsitu(insitu);
1240  Table& tout = out->table();
1241  // Get reference Epoch to time of first row or given String
1242  Unit DAY(String("d"));
1243  MEpoch::Ref epochRef(in->getTimeReference());
1244  MEpoch refEpoch;
1245  if (refTime.length()>0) {
1246    Quantum<Double> qt;
1247    if (MVTime::read(qt,refTime)) {
1248      MVEpoch mv(qt);
1249      refEpoch = MEpoch(mv, epochRef);
1250   } else {
1251      throw(AipsError("Invalid format for Epoch string"));
1252   }
1253  } else {
1254    refEpoch = in->timeCol_(0);
1255  }
1256  MPosition refPos = in->getAntennaPosition();
1257
1258  InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
1259  // test if user frame is different to base frame
1260  if ( in->frequencies().getFrameString(true)
1261       == in->frequencies().getFrameString(false) ) {
1262    throw(AipsError("Can't convert as no output frame has been set"
1263                    " (use set_freqframe) or it is aligned already."));
1264  }
1265  MFrequency::Types system = in->frequencies().getFrame();
1266  MVTime mvt(refEpoch.getValue());
1267  String epochout = mvt.string(MVTime::YMD) + String(" (") + refEpoch.getRefString() + String(")");
1268  ostringstream oss;
1269  oss << "Aligned at reference Epoch " << epochout
1270      << " in frame " << MFrequency::showType(system);
1271  pushLog(String(oss));
1272  // set up the iterator
1273  Block<String> cols(4);
1274  // select by constant direction
1275  cols[0] = String("SRCNAME");
1276  cols[1] = String("BEAMNO");
1277  // select by IF ( no of channels varies over this )
1278  cols[2] = String("IFNO");
1279  // select by restfrequency
1280  cols[3] = String("MOLECULE_ID");
1281  TableIterator iter(tout, cols);
1282  while ( !iter.pastEnd() ) {
1283    Table t = iter.table();
1284    MDirection::ROScalarColumn dirCol(t, "DIRECTION");
1285    TableIterator fiter(t, "FREQ_ID");
1286    // determine nchan from the first row. This should work as
1287    // we are iterating over BEAMNO and IFNO    // we should have constant direction
1288
1289    ROArrayColumn<Float> sCol(t, "SPECTRA");
1290    MDirection direction = dirCol(0);
1291    uInt nchan = sCol(0).nelements();
1292    while ( !fiter.pastEnd() ) {
1293      Table ftab = fiter.table();
1294      ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
1295      // get the SpectralCoordinate for the freqid, which we are iterating over
1296      SpectralCoordinate sC = in->frequencies().getSpectralCoordinate(freqidCol(0));
1297      FrequencyAligner<Float> fa( sC, nchan, refEpoch,
1298                                  direction, refPos, system );
1299      // realign the SpectralCoordinate and put into the output Scantable
1300      Vector<String> units(1);
1301      units = String("Hz");
1302      Bool linear=True;
1303      SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
1304      sc2.setWorldAxisUnits(units);
1305      uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
1306                                            sc2.referenceValue()[0],
1307                                            sc2.increment()[0]);
1308      TableVector<uInt> tvec(ftab, "FREQ_ID");
1309      tvec = id;
1310      // create the "global" abcissa for alignment with same FREQ_ID
1311      Vector<Double> abc(nchan);
1312      Double w;
1313      for (uInt i=0; i<nchan; i++) {
1314        sC.toWorld(w,Double(i));
1315        abc[i] = w;
1316      }
1317      // cache abcissa for same time stamps, so iterate over those
1318      TableIterator timeiter(ftab, "TIME");
1319      while ( !timeiter.pastEnd() ) {
1320        Table tab = timeiter.table();
1321        ArrayColumn<Float> specCol(tab, "SPECTRA");
1322        ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1323        MEpoch::ROScalarColumn timeCol(tab, "TIME");
1324        // use align abcissa cache after the first row
1325        bool first = true;
1326        // these rows should be just be POLNO
1327        for (int i=0; i<int(tab.nrow()); ++i) {
1328          // input values
1329          Vector<uChar> flag = flagCol(i);
1330          Vector<Bool> mask(flag.shape());
1331          Vector<Float> specOut, spec;
1332          spec  = specCol(i);
1333          Vector<Bool> maskOut;Vector<uChar> flagOut;
1334          convertArray(mask, flag);
1335          // alignment
1336          Bool ok = fa.align(specOut, maskOut, abc, spec,
1337                             mask, timeCol(i), !first,
1338                             interp, False);
1339          // back into scantable
1340          flagOut.resize(maskOut.nelements());
1341          convertArray(flagOut, maskOut);
1342          flagCol.put(i, flagOut);
1343          specCol.put(i, specOut);
1344          // start abcissa caching
1345          first = false;
1346        }
1347        // next timestamp
1348        ++timeiter;
1349      }
1350      // next FREQ_ID
1351      ++fiter;
1352    }
1353    // next aligner
1354    ++iter;
1355  }
1356  // set this afterwards to ensure we are doing insitu correctly.
1357  out->frequencies().setFrame(system, true);
1358  return out;
1359}
1360
1361CountedPtr<Scantable>
1362  asap::STMath::convertPolarisation( const CountedPtr<Scantable>& in,
1363                                     const std::string & newtype )
1364{
1365  if (in->npol() != 2 && in->npol() != 4)
1366    throw(AipsError("Can only convert two or four polarisations."));
1367  if ( in->getPolType() == newtype )
1368    throw(AipsError("No need to convert."));
1369  if ( ! in->selector_.empty() )
1370    throw(AipsError("Can only convert whole scantable. Unset the selection."));
1371  bool insitu = insitu_;
1372  setInsitu(false);
1373  CountedPtr< Scantable > out = getScantable(in, true);
1374  setInsitu(insitu);
1375  Table& tout = out->table();
1376  tout.rwKeywordSet().define("POLTYPE", String(newtype));
1377
1378  Block<String> cols(4);
1379  cols[0] = "SCANNO";
1380  cols[1] = "CYCLENO";
1381  cols[2] = "BEAMNO";
1382  cols[3] = "IFNO";
1383  TableIterator it(in->originalTable_, cols);
1384  String basetype = in->getPolType();
1385  STPol* stpol = STPol::getPolClass(in->factories_, basetype);
1386  try {
1387    while ( !it.pastEnd() ) {
1388      Table tab = it.table();
1389      uInt row = tab.rowNumbers()[0];
1390      stpol->setSpectra(in->getPolMatrix(row));
1391      Float fang,fhand,parang;
1392      fang = in->focusTable_.getTotalFeedAngle(in->mfocusidCol_(row));
1393      fhand = in->focusTable_.getFeedHand(in->mfocusidCol_(row));
1394      parang = in->paraCol_(row);
1395      /// @todo re-enable this
1396      // disable total feed angle to support paralactifying Caswell style
1397      stpol->setPhaseCorrections(parang, -parang, fhand);
1398      Int npolout = 0;
1399      for (uInt i=0; i<tab.nrow(); ++i) {
1400        Vector<Float> outvec = stpol->getSpectrum(i, newtype);
1401        if ( outvec.nelements() > 0 ) {
1402          tout.addRow();
1403          TableCopy::copyRows(tout, tab, tout.nrow()-1, 0, 1);
1404          ArrayColumn<Float> sCol(tout,"SPECTRA");
1405          ScalarColumn<uInt> pCol(tout,"POLNO");
1406          sCol.put(tout.nrow()-1 ,outvec);
1407          pCol.put(tout.nrow()-1 ,uInt(npolout));
1408          npolout++;
1409       }
1410      }
1411      tout.rwKeywordSet().define("nPol", npolout);
1412      ++it;
1413    }
1414  } catch (AipsError& e) {
1415    delete stpol;
1416    throw(e);
1417  }
1418  delete stpol;
1419  return out;
1420}
1421
1422CountedPtr< Scantable >
1423  asap::STMath::mxExtract( const CountedPtr< Scantable > & in,
1424                           const std::string & scantype )
1425{
1426  bool insitu = insitu_;
1427  setInsitu(false);
1428  CountedPtr< Scantable > out = getScantable(in, true);
1429  setInsitu(insitu);
1430  Table& tout = out->table();
1431  std::string taql = "SELECT FROM $1 WHERE BEAMNO != REFBEAMNO";
1432  if (scantype == "on") {
1433    taql = "SELECT FROM $1 WHERE BEAMNO == REFBEAMNO";
1434  }
1435  Table tab = tableCommand(taql, in->table());
1436  TableCopy::copyRows(tout, tab);
1437  if (scantype == "on") {
1438    // re-index SCANNO to 0
1439    TableVector<uInt> vec(tout, "SCANNO");
1440    vec = 0;
1441  }
1442  return out;
1443}
1444
1445CountedPtr< Scantable >
1446  asap::STMath::lagFlag( const CountedPtr< Scantable > & in,
1447                          double frequency, double width )
1448{
1449  CountedPtr< Scantable > out = getScantable(in, false);
1450  Table& tout = out->table();
1451  TableIterator iter(tout, "FREQ_ID");
1452  FFTServer<Float,Complex> ffts;
1453  while ( !iter.pastEnd() ) {
1454    Table tab = iter.table();
1455    Double rp,rv,inc;
1456    ROTableRow row(tab);
1457    const TableRecord& rec = row.get(0);
1458    uInt freqid = rec.asuInt("FREQ_ID");
1459    out->frequencies().getEntry(rp, rv, inc, freqid);
1460    ArrayColumn<Float> specCol(tab, "SPECTRA");
1461    ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1462    for (int i=0; i<int(tab.nrow()); ++i) {
1463      Vector<Float> spec = specCol(i);
1464      Vector<uChar> flag = flagCol(i);
1465      Int lag0 = Int(spec.nelements()*abs(inc)/(frequency+width)+0.5);
1466      Int lag1 = Int(spec.nelements()*abs(inc)/(frequency-width)+0.5);
1467      for (int k=0; k < flag.nelements(); ++k ) {
1468        if (flag[k] > 0) {
1469          spec[k] = 0.0;
1470        }
1471      }
1472      Vector<Complex> lags;
1473      ffts.fft0(lags, spec);
1474      Int start =  max(0, lag0);
1475      Int end =  min(Int(lags.nelements()-1), lag1);
1476      if (start == end) {
1477        lags[start] = Complex(0.0);
1478      } else {
1479        for (int j=start; j <=end ;++j) {
1480          lags[j] = Complex(0.0);
1481        }
1482      }
1483      ffts.fft0(spec, lags);
1484      specCol.put(i, spec);
1485    }
1486    ++iter;
1487  }
1488  return out;
1489}
Note: See TracBrowser for help on using the repository browser.