source: trunk/src/STMath.cpp@ 956

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

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