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

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

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