source: tags/Release2.1.0b/src/STMath.cpp @ 1227

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

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