source: trunk/src/STMath.cpp@ 858

Last change on this file since 858 was 858, checked in by mar637, 20 years ago

changed average interface to use std::vector<bool> instead of casa::Vector
added default arguments to average
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 27.4 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/ExprNode.h>
27#include <tables/Tables/TableRecord.h>
28#include <tables/Tables/ReadAsciiTable.h>
29
30#include <lattices/Lattices/LatticeUtilities.h>
31
32#include <scimath/Mathematics/VectorKernel.h>
33#include <scimath/Mathematics/Convolver.h>
34#include <scimath/Functionals/Polynomial.h>
35
36#include "MathUtils.h"
37#include "RowAccumulator.h"
38#include "SDAttr.h"
39#include "STMath.h"
40
41using namespace casa;
42
43using namespace asap;
44
45STMath::STMath(bool insitu) :
46 insitu_(insitu)
47{
48}
49
50
51STMath::~STMath()
52{
53}
54
55CountedPtr<Scantable>
56STMath::average( const std::vector<CountedPtr<Scantable> >& in,
57 const std::vector<bool>& mask,
58 const std::string& weight,
59 const std::string& avmode,
60 bool alignfreq)
61{
62 if ( avmode == "SCAN" && in.size() != 1 )
63 throw(AipsError("Can't perform 'SCAN' averaging on multiple tables"));
64 WeightType wtype = stringToWeight(weight);
65 // output
66 // clone as this is non insitu
67 bool insitu = insitu_;
68 setInsitu(false);
69 CountedPtr< Scantable > out = getScantable(in[0], true);
70 setInsitu(insitu);
71
72 Table& tout = out->table();
73
74 /// @todo check if all scantables are conformant
75
76 ArrayColumn<Float> specColOut(tout,"SPECTRA");
77 ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
78 ArrayColumn<Float> tsysColOut(tout,"TSYS");
79 ScalarColumn<Double> mjdColOut(tout,"TIME");
80 ScalarColumn<Double> intColOut(tout,"INTERVAL");
81
82 // set up the output table rows. These are based on the structure of the
83 // FIRST scantabel in the vector
84 const Table& baset = in[0]->table();
85
86 Block<String> cols(3);
87 cols[0] = String("BEAMNO");
88 cols[1] = String("IFNO");
89 cols[2] = String("POLNO");
90 if ( avmode == "SOURCE" ) {
91 cols.resize(4);
92 cols[3] = String("SRCNAME");
93 }
94 if ( avmode == "SCAN" && in.size() == 1) {
95 cols.resize(4);
96 cols[3] = String("SCANNO");
97 }
98 uInt outrowCount = 0;
99 TableIterator iter(baset, cols);
100 while (!iter.pastEnd()) {
101 Table subt = iter.table();
102 // copy the first row of this selection into the new table
103 tout.addRow();
104 TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
105 ++outrowCount;
106 ++iter;
107 }
108 RowAccumulator acc(wtype);
109 Vector<Bool> cmask(mask);
110 acc.setUserMask(cmask);
111 ROTableRow row(tout);
112 ROArrayColumn<Float> specCol, tsysCol;
113 ROArrayColumn<uChar> flagCol;
114 ROScalarColumn<Double> mjdCol, intCol;
115 ROScalarColumn<Int> scanIDCol;
116
117 for (uInt i=0; i < tout.nrow(); ++i) {
118 for ( int j=0; j < in.size(); ++j ) {
119 const Table& tin = in[j]->table();
120 const TableRecord& rec = row.get(i);
121 ROScalarColumn<Double> tmp(tin, "TIME");
122 Double td;tmp.get(0,td);
123 Table basesubt = tin(tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
124 && tin.col("IFNO") == Int(rec.asuInt("IFNO"))
125 && tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
126 Table subt;
127 if ( avmode == "SOURCE") {
128 subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME") );
129 } else if (avmode == "SCAN") {
130 subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
131 } else {
132 subt = basesubt;
133 }
134 specCol.attach(subt,"SPECTRA");
135 flagCol.attach(subt,"FLAGTRA");
136 tsysCol.attach(subt,"TSYS");
137 intCol.attach(subt,"INTERVAL");
138 mjdCol.attach(subt,"TIME");
139 Vector<Float> spec,tsys;
140 Vector<uChar> flag;
141 Double inter,time;
142 for (uInt k = 0; k < subt.nrow(); ++k ) {
143 flagCol.get(k, flag);
144 Vector<Bool> bflag(flag.shape());
145 convertArray(bflag, flag);
146 if ( allEQ(bflag, True) ) {
147 continue;//don't accumulate
148 }
149 specCol.get(k, spec);
150 tsysCol.get(k, tsys);
151 intCol.get(k, inter);
152 mjdCol.get(k, time);
153 // spectrum has to be added last to enable weighting by the other values
154 acc.add(spec, !bflag, tsys, inter, time);
155 }
156 }
157 //write out
158 specColOut.put(i, acc.getSpectrum());
159 const Vector<Bool>& msk = acc.getMask();
160 Vector<uChar> flg(msk.shape());
161 convertArray(flg, !msk);
162 flagColOut.put(i, flg);
163 tsysColOut.put(i, acc.getTsys());
164 intColOut.put(i, acc.getInterval());
165 mjdColOut.put(i, acc.getTime());
166 acc.reset();
167 }
168 return out;
169}
170
171CountedPtr< Scantable > STMath::getScantable(const CountedPtr< Scantable >& in,
172 bool droprows)
173{
174 if (insitu_) return in;
175 else {
176 // clone
177 Scantable* tabp = new Scantable(*in, Bool(droprows));
178 return CountedPtr<Scantable>(tabp);
179 }
180}
181
182CountedPtr< Scantable > STMath::unaryOperate( const CountedPtr< Scantable >& in,
183 float val,
184 const std::string& mode,
185 bool tsys )
186{
187 // modes are "ADD" and "MUL"
188 CountedPtr< Scantable > out = getScantable(in, false);
189 Table& tab = out->table();
190 ArrayColumn<Float> specCol(tab,"SPECTRA");
191 ArrayColumn<Float> tsysCol(tab,"TSYS");
192 for (uInt i=0; i<tab.nrow(); ++i) {
193 Vector<Float> spec;
194 Vector<Float> ts;
195 specCol.get(i, spec);
196 tsysCol.get(i, ts);
197 if (mode == "MUL") {
198 spec *= val;
199 specCol.put(i, spec);
200 if ( tsys ) {
201 ts *= val;
202 tsysCol.put(i, ts);
203 }
204 } else if ( mode == "ADD" ) {
205 spec += val;
206 specCol.put(i, spec);
207 if ( tsys ) {
208 ts += val;
209 tsysCol.put(i, ts);
210 }
211 }
212 }
213 return out;
214}
215
216MaskedArray<Float> STMath::maskedArray( const Vector<Float>& s,
217 const Vector<uChar>& f)
218{
219 Vector<Bool> mask;
220 mask.resize(f.shape());
221 convertArray(mask, f);
222 return MaskedArray<Float>(s,!mask);
223}
224
225Vector<uChar> STMath::flagsFromMA(const MaskedArray<Float>& ma)
226{
227 const Vector<Bool>& m = ma.getMask();
228 Vector<uChar> flags(m.shape());
229 convertArray(flags, !m);
230 return flags;
231}
232
233CountedPtr< Scantable > STMath::quotient( const CountedPtr< Scantable >& in,
234 const std::string & mode,
235 bool preserve )
236{
237 /// @todo make other modes available
238 /// modes should be "nearest", "pair"
239 // make this operation non insitu
240 const Table& tin = in->table();
241 Table ons = tin(tin.col("SRCTYPE") == Int(0));
242 Table offs = tin(tin.col("SRCTYPE") == Int(1));
243 if ( offs.nrow() == 0 )
244 throw(AipsError("No 'off' scans present."));
245 // put all "on" scans into output table
246
247 bool insitu = insitu_;
248 setInsitu(false);
249 CountedPtr< Scantable > out = getScantable(in, true);
250 setInsitu(insitu);
251 Table& tout = out->table();
252
253 TableCopy::copyRows(tout, ons);
254 TableRow row(tout);
255 ROScalarColumn<Double> offtimeCol(offs, "TIME");
256
257 ArrayColumn<Float> outspecCol(tout, "SPECTRA");
258 ROArrayColumn<Float> outtsysCol(tout, "TSYS");
259 ArrayColumn<uChar> outflagCol(tout, "FLAGTRA");
260
261 for (uInt i=0; i < tout.nrow(); ++i) {
262 const TableRecord& rec = row.get(i);
263 Double ontime = rec.asDouble("TIME");
264 ROScalarColumn<Double> offtimeCol(offs, "TIME");
265 Double mindeltat = min(abs(offtimeCol.getColumn() - ontime));
266 Table sel = offs( abs(offs.col("TIME")-ontime) <= mindeltat
267 && offs.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
268 && offs.col("IFNO") == Int(rec.asuInt("IFNO"))
269 && offs.col("POLNO") == Int(rec.asuInt("POLNO")) );
270
271 TableRow offrow(sel);
272 const TableRecord& offrec = offrow.get(0);//should only be one row
273 RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
274 RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
275 RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
276 /// @fixme this assumes tsys is a scalar not vector
277 Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
278 Vector<Float> specon, tsyson;
279 outtsysCol.get(i, tsyson);
280 outspecCol.get(i, specon);
281 Vector<uChar> flagon;
282 outflagCol.get(i, flagon);
283 MaskedArray<Float> mon = maskedArray(specon, flagon);
284 MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
285 MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
286 if (preserve) {
287 quot -= tsysoffscalar;
288 } else {
289 quot -= tsyson[0];
290 }
291 outspecCol.put(i, quot.getArray());
292 outflagCol.put(i, flagsFromMA(quot));
293 }
294 return out;
295}
296
297CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
298{
299 // make copy or reference
300 CountedPtr< Scantable > out = getScantable(in, false);
301 Table& tout = out->table();
302 Block<String> cols(3);
303 cols[0] = String("SCANNO");
304 cols[1] = String("BEAMNO");
305 cols[2] = String("POLNO");
306 TableIterator iter(tout, cols);
307 while (!iter.pastEnd()) {
308 Table subt = iter.table();
309 // this should leave us with two rows for the two IFs....if not ignore
310 if (subt.nrow() != 2 ) {
311 continue;
312 }
313 ArrayColumn<Float> specCol(tout, "SPECTRA");
314 ArrayColumn<Float> tsysCol(tout, "TSYS");
315 ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
316 Vector<Float> onspec,offspec, ontsys, offtsys;
317 Vector<uChar> onflag, offflag;
318 tsysCol.get(0, ontsys); tsysCol.get(1, offtsys);
319 specCol.get(0, onspec); specCol.get(1, offspec);
320 flagCol.get(0, onflag); flagCol.get(1, offflag);
321 MaskedArray<Float> on = maskedArray(onspec, onflag);
322 MaskedArray<Float> off = maskedArray(offspec, offflag);
323 MaskedArray<Float> oncopy = on.copy();
324
325 on /= off; on -= 1.0f;
326 on *= ontsys[0];
327 off /= oncopy; off -= 1.0f;
328 off *= offtsys[0];
329 specCol.put(0, on.getArray());
330 const Vector<Bool>& m0 = on.getMask();
331 Vector<uChar> flags0(m0.shape());
332 convertArray(flags0, !m0);
333 flagCol.put(0, flags0);
334
335 specCol.put(1, off.getArray());
336 const Vector<Bool>& m1 = off.getMask();
337 Vector<uChar> flags1(m1.shape());
338 convertArray(flags1, !m1);
339 flagCol.put(1, flags1);
340
341 }
342
343 return out;
344}
345
346std::vector< float > STMath::statistic( const CountedPtr< Scantable > & in,
347 const std::vector< bool > & mask,
348 const std::string& which )
349{
350
351 Vector<Bool> m(mask);
352 const Table& tab = in->table();
353 ROArrayColumn<Float> specCol(tab, "SPECTRA");
354 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
355 std::vector<float> out;
356 for (uInt i=0; i < tab.nrow(); ++i ) {
357 Vector<Float> spec; specCol.get(i, spec);
358 MaskedArray<Float> ma = maskedArray(spec, flagCol(i));
359 float outstat;
360 if ( spec.nelements() == m.nelements() ) {
361 outstat = mathutil::statistics(which, ma(m));
362 } else {
363 outstat = mathutil::statistics(which, ma);
364 }
365 out.push_back(outstat);
366 }
367 return out;
368}
369
370CountedPtr< Scantable > STMath::bin( const CountedPtr< Scantable > & in,
371 int width )
372{
373 if ( !in->getSelection().empty() ) throw(AipsError("Can't bin subset of the data."));
374 CountedPtr< Scantable > out = getScantable(in, false);
375 Table& tout = out->table();
376 out->frequencies().rescale(width, "BIN");
377 ArrayColumn<Float> specCol(tout, "SPECTRA");
378 ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
379 for (uInt i=0; i < tout.nrow(); ++i ) {
380 MaskedArray<Float> main = maskedArray(specCol(i), flagCol(i));
381 MaskedArray<Float> maout;
382 LatticeUtilities::bin(maout, main, 0, Int(width));
383 /// @todo implement channel based tsys binning
384 specCol.put(i, maout.getArray());
385 flagCol.put(i, flagsFromMA(maout));
386 // take only the first binned spectrum's length for the deprecated
387 // global header item nChan
388 if (i==0) tout.rwKeywordSet().define(String("nChan"),
389 Int(maout.getArray().nelements()));
390 }
391 return out;
392}
393
394CountedPtr< Scantable > STMath::resample( const CountedPtr< Scantable >& in,
395 const std::string& method,
396 float width )
397//
398// Should add the possibility of width being specified in km/s. This means
399// that for each freqID (SpectralCoordinate) we will need to convert to an
400// average channel width (say at the reference pixel). Then we would need
401// to be careful to make sure each spectrum (of different freqID)
402// is the same length.
403//
404{
405 InterpolateArray1D<Double,Float>::InterpolationMethod interp;
406 Int interpMethod(stringToIMethod(method));
407
408 CountedPtr< Scantable > out = getScantable(in, false);
409 Table& tout = out->table();
410
411// Resample SpectralCoordinates (one per freqID)
412 out->frequencies().rescale(width, "RESAMPLE");
413 TableIterator iter(tout, "IFNO");
414 TableRow row(tout);
415 while ( !iter.pastEnd() ) {
416 Table tab = iter.table();
417 ArrayColumn<Float> specCol(tab, "SPECTRA");
418 //ArrayColumn<Float> tsysCol(tout, "TSYS");
419 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
420 Vector<Float> spec;
421 Vector<uChar> flag;
422 specCol.get(0,spec); // the number of channels should be constant per IF
423 uInt nChanIn = spec.nelements();
424 Vector<Float> xIn(nChanIn); indgen(xIn);
425 Int fac = Int(nChanIn/width);
426 Vector<Float> xOut(fac+10); // 10 to be safe - resize later
427 uInt k = 0;
428 Float x = 0.0;
429 while (x < Float(nChanIn) ) {
430 xOut(k) = x;
431 k++;
432 x += width;
433 }
434 uInt nChanOut = k;
435 xOut.resize(nChanOut, True);
436 // process all rows for this IFNO
437 Vector<Float> specOut;
438 Vector<Bool> maskOut;
439 Vector<uChar> flagOut;
440 for (uInt i=0; i < tab.nrow(); ++i) {
441 specCol.get(i, spec);
442 flagCol.get(i, flag);
443 Vector<Bool> mask(flag.nelements());
444 convertArray(mask, flag);
445
446 IPosition shapeIn(spec.shape());
447 //sh.nchan = nChanOut;
448 InterpolateArray1D<Float,Float>::interpolate(specOut, maskOut, xOut,
449 xIn, spec, mask,
450 interpMethod, True, True);
451 /// @todo do the same for channel based Tsys
452 flagOut.resize(maskOut.nelements());
453 convertArray(flagOut, maskOut);
454 specCol.put(i, specOut);
455 flagCol.put(i, flagOut);
456 }
457 ++iter;
458 }
459
460 return out;
461}
462
463STMath::imethod STMath::stringToIMethod(const std::string& in)
464{
465 static STMath::imap lookup;
466
467 // initialize the lookup table if necessary
468 if ( lookup.empty() ) {
469 lookup["NEAR"] = InterpolateArray1D<Double,Float>::nearestNeighbour;
470 lookup["LIN"] = InterpolateArray1D<Double,Float>::linear;
471 lookup["CUB"] = InterpolateArray1D<Double,Float>::cubic;
472 lookup["SPL"] = InterpolateArray1D<Double,Float>::spline;
473 }
474
475 STMath::imap::const_iterator iter = lookup.find(in);
476
477 if ( lookup.end() == iter ) {
478 std::string message = in;
479 message += " is not a valid interpolation mode";
480 throw(AipsError(message));
481 }
482 return iter->second;
483}
484
485WeightType STMath::stringToWeight(const std::string& in)
486{
487 static std::map<std::string, WeightType> lookup;
488
489 // initialize the lookup table if necessary
490 if ( lookup.empty() ) {
491 lookup["NONE"] = asap::NONE;
492 lookup["TINT"] = asap::TINT;
493 lookup["TINTSYS"] = asap::TINTSYS;
494 lookup["TSYS"] = asap::TSYS;
495 lookup["VAR"] = asap::VAR;
496 }
497
498 std::map<std::string, WeightType>::const_iterator iter = lookup.find(in);
499
500 if ( lookup.end() == iter ) {
501 std::string message = in;
502 message += " is not a valid weighting mode";
503 throw(AipsError(message));
504 }
505 return iter->second;
506}
507
508CountedPtr< Scantable > STMath::gainElevation( const CountedPtr< Scantable >& in,
509 const Vector< Float > & coeffs,
510 const std::string & filename,
511 const std::string& method)
512{
513 // Get elevation data from Scantable and convert to degrees
514 CountedPtr< Scantable > out = getScantable(in, false);
515 Table& tab = in->table();
516 ROScalarColumn<Float> elev(tab, "ELEVATION");
517 Vector<Float> x = elev.getColumn();
518 x *= Float(180 / C::pi); // Degrees
519
520 const uInt nc = coeffs.nelements();
521 if ( filename.length() > 0 && nc > 0 ) {
522 throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
523 }
524
525 // Correct
526 if ( nc > 0 || filename.length() == 0 ) {
527 // Find instrument
528 Bool throwit = True;
529 Instrument inst =
530 SDAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
531 throwit);
532
533 // Set polynomial
534 Polynomial<Float>* ppoly = 0;
535 Vector<Float> coeff;
536 String msg;
537 if ( nc > 0 ) {
538 ppoly = new Polynomial<Float>(nc);
539 coeff = coeffs;
540 msg = String("user");
541 } else {
542 SDAttr sdAttr;
543 coeff = sdAttr.gainElevationPoly(inst);
544 ppoly = new Polynomial<Float>(3);
545 msg = String("built in");
546 }
547
548 if ( coeff.nelements() > 0 ) {
549 ppoly->setCoefficients(coeff);
550 } else {
551 delete ppoly;
552 throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
553 }
554 ostringstream oss;
555 oss << "Making polynomial correction with " << msg << " coefficients:" << endl;
556 oss << " " << coeff;
557 pushLog(String(oss));
558 const uInt nrow = tab.nrow();
559 Vector<Float> factor(nrow);
560 for ( uInt i=0; i < nrow; ++i ) {
561 factor[i] = 1.0 / (*ppoly)(x[i]);
562 }
563 delete ppoly;
564 scaleByVector(tab, factor, true);
565
566 } else {
567 // Read and correct
568 pushLog("Making correction from ascii Table");
569 scaleFromAsciiTable(tab, filename, method, x, true);
570 }
571 return out;
572}
573
574void STMath::scaleFromAsciiTable(Table& in, const std::string& filename,
575 const std::string& method,
576 const Vector<Float>& xout, bool dotsys)
577{
578
579// Read gain-elevation ascii file data into a Table.
580
581 String formatString;
582 Table tbl = readAsciiTable(formatString, Table::Memory, filename, "", "", False);
583 scaleFromTable(in, tbl, method, xout, dotsys);
584}
585
586void STMath::scaleFromTable(Table& in,
587 const Table& table,
588 const std::string& method,
589 const Vector<Float>& xout, bool dotsys)
590{
591
592 ROScalarColumn<Float> geElCol(table, "ELEVATION");
593 ROScalarColumn<Float> geFacCol(table, "FACTOR");
594 Vector<Float> xin = geElCol.getColumn();
595 Vector<Float> yin = geFacCol.getColumn();
596 Vector<Bool> maskin(xin.nelements(),True);
597
598 // Interpolate (and extrapolate) with desired method
599
600 //InterpolateArray1D<Double,Float>::InterpolationMethod method;
601 Int intmethod(stringToIMethod(method));
602
603 Vector<Float> yout;
604 Vector<Bool> maskout;
605 InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout,
606 xin, yin, maskin, intmethod,
607 True, True);
608
609 scaleByVector(in, Float(1.0)/yout, dotsys);
610}
611
612void STMath::scaleByVector( Table& in,
613 const Vector< Float >& factor,
614 bool dotsys )
615{
616 uInt nrow = in.nrow();
617 if ( factor.nelements() != nrow ) {
618 throw(AipsError("factors.nelements() != table.nelements()"));
619 }
620 ArrayColumn<Float> specCol(in, "SPECTRA");
621 ArrayColumn<uChar> flagCol(in, "FLAGTRA");
622 ArrayColumn<Float> tsysCol(in, "TSYS");
623 for (uInt i=0; i < nrow; ++i) {
624 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
625 ma *= factor[i];
626 specCol.put(i, ma.getArray());
627 flagCol.put(i, flagsFromMA(ma));
628 if ( dotsys ) {
629 Vector<Float> tsys;
630 tsysCol.get(i, tsys);
631 tsys *= factor[i];
632 specCol.put(i,tsys);
633 }
634 }
635}
636
637CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in,
638 float d, float etaap,
639 float jyperk )
640{
641 CountedPtr< Scantable > out = getScantable(in, false);
642 Table& tab = in->table();
643 Unit fluxUnit(tab.keywordSet().asString("FluxUnit"));
644 Unit K(String("K"));
645 Unit JY(String("Jy"));
646
647 bool tokelvin = true;
648 Double cfac = 1.0;
649
650 if ( fluxUnit == JY ) {
651 pushLog("Converting to K");
652 Quantum<Double> t(1.0,fluxUnit);
653 Quantum<Double> t2 = t.get(JY);
654 cfac = (t2 / t).getValue(); // value to Jy
655
656 tokelvin = true;
657 out->setFluxUnit("K");
658 } else if ( fluxUnit == K ) {
659 pushLog("Converting to Jy");
660 Quantum<Double> t(1.0,fluxUnit);
661 Quantum<Double> t2 = t.get(K);
662 cfac = (t2 / t).getValue(); // value to K
663
664 tokelvin = false;
665 out->setFluxUnit("Jy");
666 } else {
667 throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
668 }
669 // Make sure input values are converted to either Jy or K first...
670 Float factor = cfac;
671
672 // Select method
673 if (jyperk > 0.0) {
674 factor *= jyperk;
675 if ( tokelvin ) factor = 1.0 / jyperk;
676 ostringstream oss;
677 oss << "Jy/K = " << jyperk;
678 pushLog(String(oss));
679 Vector<Float> factors(tab.nrow(), factor);
680 scaleByVector(tab,factors, false);
681 } else if ( etaap > 0.0) {
682 Instrument inst =
683 SDAttr::convertInstrument(tab.keywordSet().asString("AntennaName"), True);
684 SDAttr sda;
685 if (d < 0) d = sda.diameter(inst);
686 Float jyPerk = SDAttr::findJyPerK(etaap, d);
687 ostringstream oss;
688 oss << "Jy/K = " << jyperk;
689 pushLog(String(oss));
690 factor *= jyperk;
691 if ( tokelvin ) {
692 factor = 1.0 / factor;
693 }
694 Vector<Float> factors(tab.nrow(), factor);
695 scaleByVector(tab, factors, False);
696 } else {
697
698 // OK now we must deal with automatic look up of values.
699 // We must also deal with the fact that the factors need
700 // to be computed per IF and may be different and may
701 // change per integration.
702
703 pushLog("Looking up conversion factors");
704 convertBrightnessUnits(out, tokelvin, cfac);
705 }
706
707 return out;
708}
709
710void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in,
711 bool tokelvin, float cfac )
712{
713 Table& table = in->table();
714 Instrument inst =
715 SDAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True);
716 TableIterator iter(table, "FREQ_ID");
717 STFrequencies stfreqs = in->frequencies();
718 SDAttr sdAtt;
719 while (!iter.pastEnd()) {
720 Table tab = iter.table();
721 ArrayColumn<Float> specCol(tab, "SPECTRA");
722 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
723 ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID");
724 MEpoch::ROScalarColumn timeCol(tab, "TIME");
725
726 uInt freqid; freqidCol.get(0, freqid);
727 Vector<Float> tmpspec; specCol.get(0, tmpspec);
728 // SDAttr.JyPerK has a Vector interface... change sometime.
729 Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements()));
730 for ( uInt i=0; i<tab.nrow(); ++i) {
731 Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0];
732 Float factor = cfac * jyperk;
733 if ( tokelvin ) factor = Float(1.0) / factor;
734 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
735 ma *= factor;
736 specCol.put(i, ma.getArray());
737 flagCol.put(i, flagsFromMA(ma));
738 }
739 }
740}
741
742CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
743 float tau )
744{
745 CountedPtr< Scantable > out = getScantable(in, false);
746 Table& tab = in->table();
747 ROScalarColumn<Float> elev(tab, "ELEVATION");
748 ArrayColumn<Float> specCol(tab, "SPECTRA");
749 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
750 for ( uInt i=0; i<tab.nrow(); ++i) {
751 Float zdist = Float(C::pi_2) - elev(i);
752 Float factor = exp(tau)/cos(zdist);
753 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
754 ma *= factor;
755 specCol.put(i, ma.getArray());
756 flagCol.put(i, flagsFromMA(ma));
757 }
758 return out;
759}
760
761CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
762 const std::string& kernel, float width )
763{
764 CountedPtr< Scantable > out = getScantable(in, false);
765 Table& table = in->table();
766 VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
767 // same IFNO should have same no of channels
768 // this saves overhead
769 TableIterator iter(table, "IFNO");
770 while (!iter.pastEnd()) {
771 Table tab = iter.table();
772 ArrayColumn<Float> specCol(tab, "SPECTRA");
773 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
774 Vector<Float> tmpspec; specCol.get(0, tmpspec);
775 uInt nchan = tmpspec.nelements();
776 Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
777 Convolver<Float> conv(kvec, IPosition(1,nchan));
778 Vector<Float> spec;
779 Vector<uChar> flag;
780 for ( uInt i=0; i<tab.nrow(); ++i) {
781 specCol.get(i, spec);
782 flagCol.get(i, flag);
783 Vector<Bool> mask(flag.nelements());
784 convertArray(mask, flag);
785 Vector<Float> specout;
786 if ( type == VectorKernel::HANNING ) {
787 Vector<Bool> maskout;
788 mathutil::hanning(specout, maskout, spec , mask);
789 convertArray(flag, maskout);
790 flagCol.put(i, flag);
791 } else {
792 mathutil::replaceMaskByZero(specout, mask);
793 conv.linearConv(specout, spec);
794 }
795 specCol.put(i, specout);
796 }
797 }
798 return out;
799}
800
801CountedPtr< Scantable >
802 STMath::merge( const std::vector< CountedPtr < Scantable > >& in )
803{
804 if ( in.size() < 2 ) {
805 throw(AipsError("Need at least two scantables to perform merge."));
806 }
807 std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
808 ++it;
809 bool insitu = insitu_;
810 setInsitu(false);
811 CountedPtr< Scantable > out = getScantable(in[0], false);
812 setInsitu(insitu);
813 Table& tout = out->table();
814 ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
815 ScalarColumn<uInt> scannocol(tout,"SCANNO"),focusidcol(tout,"FOCUS_ID");
816 uInt newscanno = max(scannocol.getColumn())+1;
817 while ( it != in.end() ){
818 if ( ! (*it)->conformant(*out) ) {
819 // log message: "ignoring scantable i, as it isn't
820 // conformant with the other(s)"
821 cerr << "oh oh" << endl;
822 ++it;
823 continue;
824 }
825 const Table& tab = (*it)->table();
826 TableIterator scanit(tab, "SCANNO");
827 while (!scanit.pastEnd()) {
828 TableIterator freqit(scanit.table(), "FREQ_ID");
829 while ( !freqit.pastEnd() ) {
830 Table thetab = freqit.table();
831 uInt nrow = tout.nrow();
832 //tout.addRow(thetab.nrow());
833 TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
834 ROTableRow row(thetab);
835 for ( uInt i=0; i<thetab.nrow(); ++i) {
836 uInt k = nrow+i;
837 scannocol.put(k, newscanno);
838 const TableRecord& rec = row.get(i);
839 Double rv,rp,inc;
840 (*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
841 uInt id;
842 id = out->frequencies().addEntry(rp, rv, inc);
843 freqidcol.put(k,id);
844 String name,fname;Double rf;
845 (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
846 id = out->molecules().addEntry(rf, name, fname);
847 molidcol.put(k, id);
848 Float frot,fang,ftan;
849 (*it)->focus().getEntry(frot, fang, ftan, rec.asuInt("FOCUS_ID"));
850 id = out->focus().addEntry(frot, fang, ftan);
851 focusidcol.put(k, id);
852 }
853 ++freqit;
854 }
855 ++newscanno;
856 ++scanit;
857 }
858 ++it;
859 }
860 return out;
861}
Note: See TracBrowser for help on using the repository browser.