source: trunk/src/STMath.cpp@ 871

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

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