source: trunk/src/STMath.cpp@ 1329

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

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