source: trunk/src/STMath.cpp@ 1014

Last change on this file since 1014 was 1010, checked in by mar637, 19 years ago

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