source: trunk/src/STMath.cpp@ 934

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

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