source: trunk/src/STMath.cpp@ 1386

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

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