source: trunk/src/STMath.cpp@ 1008

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

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