source: trunk/src/STMath.cpp@ 845

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

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