source: trunk/src/STMath.cpp@ 920

Last change on this file since 920 was 917, checked in by mar637, 18 years ago

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