source: branches/alma/src/STMath.cpp@ 1515

Last change on this file since 1515 was 1514, checked in by Kana Sugimoto, 16 years ago

New Development: No

JIRA Issue: Yes (CAS-1079)

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed:
Added the new functions to get min/max position (channel) of spectra.

  1. std::vector<int> asap::python::stmath::_minmaxpos(ScantableWrapper& in,

std::vector<bool>& mask, std::string& which)

@python_STMath.cpp

  1. std::vector<int> STMathWrapper::minMaxPos(ScantableWrapper& in,

std::vector<bool>& mask, std::string& which)

@STMathWrapper.h

  1. std::vector<int> STMath::minMaxPos(CountedPtr< Scantable > & in,

std::vector< bool > & mask, std::string& which)

@STMath.h & .cpp

  1. IPosition mathutil::minMaxPos(String& which, MaskedArray<Float>& data) @MathUtils.h & .cpp

Test Programs:

Run scantable.stats() with the parameter stat='minpos' or 'maxpos',
and you'll get min/max value with its psition.

Put in Release Notes: No

Module(s): scantable.stats()

Description:

These modifications are to return min/max value with its
position (channel/frequency/velocity) by running scantable.stats().

Diagram:
scantable.stats() ->asap::python::stmath::_minmaxpos
-> STMathWrapper::minMaxPos -> STMath::minMaxPos() -> mathutil::minMaxPos()
-> casa::minMax (@casacore/casa/casa/Arrays/MaskArrMath.tcc)
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 91.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 <casa/Arrays/MaskArrLogi.h>
18#include <casa/Arrays/MaskArrMath.h>
19#include <casa/Arrays/ArrayLogical.h>
20#include <casa/Arrays/ArrayMath.h>
21#include <casa/Arrays/Slice.h>
22#include <casa/Arrays/Slicer.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/TableParse.h>
30#include <tables/Tables/ReadAsciiTable.h>
31#include <tables/Tables/TableIter.h>
32#include <tables/Tables/TableCopy.h>
33#include <scimath/Mathematics/FFTServer.h>
34
35#include <lattices/Lattices/LatticeUtilities.h>
36
37#include <coordinates/Coordinates/SpectralCoordinate.h>
38#include <coordinates/Coordinates/CoordinateSystem.h>
39#include <coordinates/Coordinates/CoordinateUtil.h>
40#include <coordinates/Coordinates/FrequencyAligner.h>
41
42#include <scimath/Mathematics/VectorKernel.h>
43#include <scimath/Mathematics/Convolver.h>
44#include <scimath/Functionals/Polynomial.h>
45
46#include "MathUtils.h"
47#include "RowAccumulator.h"
48#include "STAttr.h"
49#include "STMath.h"
50#include "STSelector.h"
51
52using namespace casa;
53
54using namespace asap;
55
56STMath::STMath(bool insitu) :
57 insitu_(insitu)
58{
59}
60
61
62STMath::~STMath()
63{
64}
65
66CountedPtr<Scantable>
67STMath::average( const std::vector<CountedPtr<Scantable> >& in,
68 const std::vector<bool>& mask,
69 const std::string& weight,
70 const std::string& avmode)
71{
72 if ( avmode == "SCAN" && in.size() != 1 )
73 throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
74 "Use merge first."));
75 WeightType wtype = stringToWeight(weight);
76
77 // output
78 // clone as this is non insitu
79 bool insitu = insitu_;
80 setInsitu(false);
81 CountedPtr< Scantable > out = getScantable(in[0], true);
82 setInsitu(insitu);
83 std::vector<CountedPtr<Scantable> >::const_iterator stit = in.begin();
84 ++stit;
85 while ( stit != in.end() ) {
86 out->appendToHistoryTable((*stit)->history());
87 ++stit;
88 }
89
90 Table& tout = out->table();
91
92 /// @todo check if all scantables are conformant
93
94 ArrayColumn<Float> specColOut(tout,"SPECTRA");
95 ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
96 ArrayColumn<Float> tsysColOut(tout,"TSYS");
97 ScalarColumn<Double> mjdColOut(tout,"TIME");
98 ScalarColumn<Double> intColOut(tout,"INTERVAL");
99 ScalarColumn<uInt> cycColOut(tout,"CYCLENO");
100 ScalarColumn<uInt> scanColOut(tout,"SCANNO");
101
102 // set up the output table rows. These are based on the structure of the
103 // FIRST scantable in the vector
104 const Table& baset = in[0]->table();
105
106 Block<String> cols(3);
107 cols[0] = String("BEAMNO");
108 cols[1] = String("IFNO");
109 cols[2] = String("POLNO");
110 if ( avmode == "SOURCE" ) {
111 cols.resize(4);
112 cols[3] = String("SRCNAME");
113 }
114 if ( avmode == "SCAN" && in.size() == 1) {
115 //cols.resize(4);
116 //cols[3] = String("SCANNO");
117 cols.resize(5);
118 cols[3] = String("SRCNAME");
119 cols[4] = String("SCANNO");
120 }
121 uInt outrowCount = 0;
122 TableIterator iter(baset, cols);
123 while (!iter.pastEnd()) {
124 Table subt = iter.table();
125 // copy the first row of this selection into the new table
126 tout.addRow();
127 TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
128 // re-index to 0
129 if ( avmode != "SCAN" && avmode != "SOURCE" ) {
130 scanColOut.put(outrowCount, uInt(0));
131 }
132 ++outrowCount;
133 ++iter;
134 }
135 RowAccumulator acc(wtype);
136 Vector<Bool> cmask(mask);
137 acc.setUserMask(cmask);
138 ROTableRow row(tout);
139 ROArrayColumn<Float> specCol, tsysCol;
140 ROArrayColumn<uChar> flagCol;
141 ROScalarColumn<Double> mjdCol, intCol;
142 ROScalarColumn<Int> scanIDCol;
143
144 Vector<uInt> rowstodelete;
145
146 for (uInt i=0; i < tout.nrow(); ++i) {
147 for ( int j=0; j < int(in.size()); ++j ) {
148 const Table& tin = in[j]->table();
149 const TableRecord& rec = row.get(i);
150 ROScalarColumn<Double> tmp(tin, "TIME");
151 Double td;tmp.get(0,td);
152 Table basesubt = tin(tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
153 && tin.col("IFNO") == Int(rec.asuInt("IFNO"))
154 && tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
155 Table subt;
156 if ( avmode == "SOURCE") {
157 subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME") );
158 } else if (avmode == "SCAN") {
159 //subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
160 subt = basesubt( basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO"))
161 && basesubt.col("SRCNAME") == rec.asString("SRCNAME") );
162 } else {
163 subt = basesubt;
164 }
165 specCol.attach(subt,"SPECTRA");
166 flagCol.attach(subt,"FLAGTRA");
167 tsysCol.attach(subt,"TSYS");
168 intCol.attach(subt,"INTERVAL");
169 mjdCol.attach(subt,"TIME");
170 Vector<Float> spec,tsys;
171 Vector<uChar> flag;
172 Double inter,time;
173 for (uInt k = 0; k < subt.nrow(); ++k ) {
174 flagCol.get(k, flag);
175 Vector<Bool> bflag(flag.shape());
176 convertArray(bflag, flag);
177 /*
178 if ( allEQ(bflag, True) ) {
179 continue;//don't accumulate
180 }
181 */
182 specCol.get(k, spec);
183 tsysCol.get(k, tsys);
184 intCol.get(k, inter);
185 mjdCol.get(k, time);
186 // spectrum has to be added last to enable weighting by the other values
187 acc.add(spec, !bflag, tsys, inter, time);
188 }
189 }
190 const Vector<Bool>& msk = acc.getMask();
191 if ( allEQ(msk, False) ) {
192 uint n = rowstodelete.nelements();
193 rowstodelete.resize(n+1, True);
194 rowstodelete[n] = i;
195 continue;
196 }
197 //write out
198 if (acc.state()) {
199 Vector<uChar> flg(msk.shape());
200 convertArray(flg, !msk);
201 flagColOut.put(i, flg);
202 specColOut.put(i, acc.getSpectrum());
203 tsysColOut.put(i, acc.getTsys());
204 intColOut.put(i, acc.getInterval());
205 mjdColOut.put(i, acc.getTime());
206 // we should only have one cycle now -> reset it to be 0
207 // frequency switched data has different CYCLENO for different IFNO
208 // which requires resetting this value
209 cycColOut.put(i, uInt(0));
210 } else {
211 ostringstream oss;
212 oss << "For output row="<<i<<", all input rows of data are flagged. no averaging" << endl;
213 pushLog(String(oss));
214 }
215 acc.reset();
216 }
217 if (rowstodelete.nelements() > 0) {
218 cout << rowstodelete << endl;
219 tout.removeRow(rowstodelete);
220 if (tout.nrow() == 0) {
221 throw(AipsError("Can't average fully flagged data."));
222 }
223 }
224 return out;
225}
226
227CountedPtr< Scantable >
228 STMath::averageChannel( const CountedPtr < Scantable > & in,
229 const std::string & mode,
230 const std::string& avmode )
231{
232 // clone as this is non insitu
233 bool insitu = insitu_;
234 setInsitu(false);
235 CountedPtr< Scantable > out = getScantable(in, true);
236 setInsitu(insitu);
237 Table& tout = out->table();
238 ArrayColumn<Float> specColOut(tout,"SPECTRA");
239 ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
240 ArrayColumn<Float> tsysColOut(tout,"TSYS");
241 ScalarColumn<uInt> scanColOut(tout,"SCANNO");
242 ScalarColumn<Double> intColOut(tout, "INTERVAL");
243 Table tmp = in->table().sort("BEAMNO");
244 Block<String> cols(3);
245 cols[0] = String("BEAMNO");
246 cols[1] = String("IFNO");
247 cols[2] = String("POLNO");
248 if ( avmode == "SCAN") {
249 cols.resize(4);
250 cols[3] = String("SCANNO");
251 }
252 uInt outrowCount = 0;
253 uChar userflag = 1 << 7;
254 TableIterator iter(tmp, cols);
255 while (!iter.pastEnd()) {
256 Table subt = iter.table();
257 ROArrayColumn<Float> specCol, tsysCol;
258 ROArrayColumn<uChar> flagCol;
259 ROScalarColumn<Double> intCol(subt, "INTERVAL");
260 specCol.attach(subt,"SPECTRA");
261 flagCol.attach(subt,"FLAGTRA");
262 tsysCol.attach(subt,"TSYS");
263 tout.addRow();
264 TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
265 if ( avmode != "SCAN") {
266 scanColOut.put(outrowCount, uInt(0));
267 }
268 Vector<Float> tmp;
269 specCol.get(0, tmp);
270 uInt nchan = tmp.nelements();
271 // have to do channel by channel here as MaskedArrMath
272 // doesn't have partialMedians
273 Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
274 Vector<Float> outspec(nchan);
275 Vector<uChar> outflag(nchan,0);
276 Vector<Float> outtsys(1);/// @fixme when tsys is channel based
277 for (uInt i=0; i<nchan; ++i) {
278 Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
279 MaskedArray<Float> ma = maskedArray(specs,flags);
280 outspec[i] = median(ma);
281 if ( allEQ(ma.getMask(), False) )
282 outflag[i] = userflag;// flag data
283 }
284 outtsys[0] = median(tsysCol.getColumn());
285 specColOut.put(outrowCount, outspec);
286 flagColOut.put(outrowCount, outflag);
287 tsysColOut.put(outrowCount, outtsys);
288 Double intsum = sum(intCol.getColumn());
289 intColOut.put(outrowCount, intsum);
290 ++outrowCount;
291 ++iter;
292 }
293 return out;
294}
295
296CountedPtr< Scantable > STMath::getScantable(const CountedPtr< Scantable >& in,
297 bool droprows)
298{
299 if (insitu_) return in;
300 else {
301 // clone
302 Scantable* tabp = new Scantable(*in, Bool(droprows));
303 return CountedPtr<Scantable>(tabp);
304 }
305}
306
307CountedPtr< Scantable > STMath::unaryOperate( const CountedPtr< Scantable >& in,
308 float val,
309 const std::string& mode,
310 bool tsys )
311{
312 CountedPtr< Scantable > out = getScantable(in, false);
313 Table& tab = out->table();
314 ArrayColumn<Float> specCol(tab,"SPECTRA");
315 ArrayColumn<Float> tsysCol(tab,"TSYS");
316 for (uInt i=0; i<tab.nrow(); ++i) {
317 Vector<Float> spec;
318 Vector<Float> ts;
319 specCol.get(i, spec);
320 tsysCol.get(i, ts);
321 if (mode == "MUL" || mode == "DIV") {
322 if (mode == "DIV") val = 1.0/val;
323 spec *= val;
324 specCol.put(i, spec);
325 if ( tsys ) {
326 ts *= val;
327 tsysCol.put(i, ts);
328 }
329 } else if ( mode == "ADD" || mode == "SUB") {
330 if (mode == "SUB") val *= -1.0;
331 spec += val;
332 specCol.put(i, spec);
333 if ( tsys ) {
334 ts += val;
335 tsysCol.put(i, ts);
336 }
337 }
338 }
339 return out;
340}
341
342CountedPtr<Scantable> STMath::binaryOperate(const CountedPtr<Scantable>& left,
343 const CountedPtr<Scantable>& right,
344 const std::string& mode)
345{
346 bool insitu = insitu_;
347 if ( ! left->conformant(*right) ) {
348 throw(AipsError("'left' and 'right' scantables are not conformant."));
349 }
350 setInsitu(false);
351 CountedPtr< Scantable > out = getScantable(left, false);
352 setInsitu(insitu);
353 Table& tout = out->table();
354 Block<String> coln(5);
355 coln[0] = "SCANNO"; coln[1] = "CYCLENO"; coln[2] = "BEAMNO";
356 coln[3] = "IFNO"; coln[4] = "POLNO";
357 Table tmpl = tout.sort(coln);
358 Table tmpr = right->table().sort(coln);
359 ArrayColumn<Float> lspecCol(tmpl,"SPECTRA");
360 ROArrayColumn<Float> rspecCol(tmpr,"SPECTRA");
361 ArrayColumn<uChar> lflagCol(tmpl,"FLAGTRA");
362 ROArrayColumn<uChar> rflagCol(tmpr,"FLAGTRA");
363
364 for (uInt i=0; i<tout.nrow(); ++i) {
365 Vector<Float> lspecvec, rspecvec;
366 Vector<uChar> lflagvec, rflagvec;
367 lspecvec = lspecCol(i); rspecvec = rspecCol(i);
368 lflagvec = lflagCol(i); rflagvec = rflagCol(i);
369 MaskedArray<Float> mleft = maskedArray(lspecvec, lflagvec);
370 MaskedArray<Float> mright = maskedArray(rspecvec, rflagvec);
371 if (mode == "ADD") {
372 mleft += mright;
373 } else if ( mode == "SUB") {
374 mleft -= mright;
375 } else if ( mode == "MUL") {
376 mleft *= mright;
377 } else if ( mode == "DIV") {
378 mleft /= mright;
379 } else {
380 throw(AipsError("Illegal binary operator"));
381 }
382 lspecCol.put(i, mleft.getArray());
383 }
384 return out;
385}
386
387
388
389MaskedArray<Float> STMath::maskedArray( const Vector<Float>& s,
390 const Vector<uChar>& f)
391{
392 Vector<Bool> mask;
393 mask.resize(f.shape());
394 convertArray(mask, f);
395 return MaskedArray<Float>(s,!mask);
396}
397
398Vector<uChar> STMath::flagsFromMA(const MaskedArray<Float>& ma)
399{
400 const Vector<Bool>& m = ma.getMask();
401 Vector<uChar> flags(m.shape());
402 convertArray(flags, !m);
403 return flags;
404}
405
406CountedPtr< Scantable > STMath::autoQuotient( const CountedPtr< Scantable >& in,
407 const std::string & mode,
408 bool preserve )
409{
410 /// @todo make other modes available
411 /// modes should be "nearest", "pair"
412 // make this operation non insitu
413 const Table& tin = in->table();
414 Table ons = tin(tin.col("SRCTYPE") == Int(0));
415 Table offs = tin(tin.col("SRCTYPE") == Int(1));
416 if ( offs.nrow() == 0 )
417 throw(AipsError("No 'off' scans present."));
418 // put all "on" scans into output table
419
420 bool insitu = insitu_;
421 setInsitu(false);
422 CountedPtr< Scantable > out = getScantable(in, true);
423 setInsitu(insitu);
424 Table& tout = out->table();
425
426 TableCopy::copyRows(tout, ons);
427 TableRow row(tout);
428 ROScalarColumn<Double> offtimeCol(offs, "TIME");
429 ArrayColumn<Float> outspecCol(tout, "SPECTRA");
430 ROArrayColumn<Float> outtsysCol(tout, "TSYS");
431 ArrayColumn<uChar> outflagCol(tout, "FLAGTRA");
432 for (uInt i=0; i < tout.nrow(); ++i) {
433 const TableRecord& rec = row.get(i);
434 Double ontime = rec.asDouble("TIME");
435 Table presel = offs(offs.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
436 && offs.col("IFNO") == Int(rec.asuInt("IFNO"))
437 && offs.col("POLNO") == Int(rec.asuInt("POLNO")) );
438 ROScalarColumn<Double> offtimeCol(presel, "TIME");
439
440 Double mindeltat = min(abs(offtimeCol.getColumn() - ontime));
441 // Timestamp may vary within a cycle ???!!!
442 // increase this by 0.01 sec in case of rounding errors...
443 // There might be a better way to do this.
444 // fix to this fix. TIME is MJD, so 1.0d not 1.0s
445 mindeltat += 0.01/24./60./60.;
446 Table sel = presel( abs(presel.col("TIME")-ontime) <= mindeltat);
447
448 if ( sel.nrow() < 1 ) {
449 throw(AipsError("No closest in time found... This could be a rounding "
450 "issue. Try quotient instead."));
451 }
452 TableRow offrow(sel);
453 const TableRecord& offrec = offrow.get(0);//should only be one row
454 RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
455 RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
456 RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
457 /// @fixme this assumes tsys is a scalar not vector
458 Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
459 Vector<Float> specon, tsyson;
460 outtsysCol.get(i, tsyson);
461 outspecCol.get(i, specon);
462 Vector<uChar> flagon;
463 outflagCol.get(i, flagon);
464 MaskedArray<Float> mon = maskedArray(specon, flagon);
465 MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
466 MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
467 if (preserve) {
468 quot -= tsysoffscalar;
469 } else {
470 quot -= tsyson[0];
471 }
472 outspecCol.put(i, quot.getArray());
473 outflagCol.put(i, flagsFromMA(quot));
474 }
475 // renumber scanno
476 TableIterator it(tout, "SCANNO");
477 uInt i = 0;
478 while ( !it.pastEnd() ) {
479 Table t = it.table();
480 TableVector<uInt> vec(t, "SCANNO");
481 vec = i;
482 ++i;
483 ++it;
484 }
485 return out;
486}
487
488
489CountedPtr< Scantable > STMath::quotient( const CountedPtr< Scantable > & on,
490 const CountedPtr< Scantable > & off,
491 bool preserve )
492{
493 bool insitu = insitu_;
494 if ( ! on->conformant(*off) ) {
495 throw(AipsError("'on' and 'off' scantables are not conformant."));
496 }
497 setInsitu(false);
498 CountedPtr< Scantable > out = getScantable(on, false);
499 setInsitu(insitu);
500 Table& tout = out->table();
501 const Table& toff = off->table();
502 TableIterator sit(tout, "SCANNO");
503 TableIterator s2it(toff, "SCANNO");
504 while ( !sit.pastEnd() ) {
505 Table ton = sit.table();
506 TableRow row(ton);
507 Table t = s2it.table();
508 ArrayColumn<Float> outspecCol(ton, "SPECTRA");
509 ROArrayColumn<Float> outtsysCol(ton, "TSYS");
510 ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
511 for (uInt i=0; i < ton.nrow(); ++i) {
512 const TableRecord& rec = row.get(i);
513 Table offsel = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
514 && t.col("IFNO") == Int(rec.asuInt("IFNO"))
515 && t.col("POLNO") == Int(rec.asuInt("POLNO")) );
516 if ( offsel.nrow() == 0 )
517 throw AipsError("STMath::quotient: no matching off");
518 TableRow offrow(offsel);
519 const TableRecord& offrec = offrow.get(0);//should be ncycles - take first
520 RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
521 RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
522 RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
523 Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
524 Vector<Float> specon, tsyson;
525 outtsysCol.get(i, tsyson);
526 outspecCol.get(i, specon);
527 Vector<uChar> flagon;
528 outflagCol.get(i, flagon);
529 MaskedArray<Float> mon = maskedArray(specon, flagon);
530 MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
531 MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
532 if (preserve) {
533 quot -= tsysoffscalar;
534 } else {
535 quot -= tsyson[0];
536 }
537 outspecCol.put(i, quot.getArray());
538 outflagCol.put(i, flagsFromMA(quot));
539 }
540 ++sit;
541 ++s2it;
542 // take the first off for each on scan which doesn't have a
543 // matching off scan
544 // non <= noff: matching pairs, non > noff matching pairs then first off
545 if ( s2it.pastEnd() ) s2it.reset();
546 }
547 return out;
548}
549
550// dototalpower (migration of GBTIDL procedure dototalpower.pro)
551// calibrate the CAL on-off pair. It calculate Tsys and average CAL on-off subintegrations
552// do it for each cycles in a specific scan.
553CountedPtr< Scantable > STMath::dototalpower( const CountedPtr< Scantable >& calon,
554 const CountedPtr< Scantable >& caloff, Float tcal )
555{
556if ( ! calon->conformant(*caloff) ) {
557 throw(AipsError("'CAL on' and 'CAL off' scantables are not conformant."));
558 }
559 setInsitu(false);
560 CountedPtr< Scantable > out = getScantable(caloff, false);
561 Table& tout = out->table();
562 const Table& tcon = calon->table();
563 Vector<Float> tcalout;
564 Vector<Float> tcalout2; //debug
565
566 if ( tout.nrow() != tcon.nrow() ) {
567 throw(AipsError("Mismatch in number of rows to form cal on - off pair."));
568 }
569 // iteration by scanno or cycle no.
570 TableIterator sit(tout, "SCANNO");
571 TableIterator s2it(tcon, "SCANNO");
572 while ( !sit.pastEnd() ) {
573 Table toff = sit.table();
574 TableRow row(toff);
575 Table t = s2it.table();
576 ScalarColumn<Double> outintCol(toff, "INTERVAL");
577 ArrayColumn<Float> outspecCol(toff, "SPECTRA");
578 ArrayColumn<Float> outtsysCol(toff, "TSYS");
579 ArrayColumn<uChar> outflagCol(toff, "FLAGTRA");
580 ROScalarColumn<uInt> outtcalIdCol(toff, "TCAL_ID");
581 ROScalarColumn<uInt> outpolCol(toff, "POLNO");
582 ROScalarColumn<Double> onintCol(t, "INTERVAL");
583 ROArrayColumn<Float> onspecCol(t, "SPECTRA");
584 ROArrayColumn<Float> ontsysCol(t, "TSYS");
585 ROArrayColumn<uChar> onflagCol(t, "FLAGTRA");
586 //ROScalarColumn<uInt> ontcalIdCol(t, "TCAL_ID");
587
588 for (uInt i=0; i < toff.nrow(); ++i) {
589 //skip these checks -> assumes the data order are the same between the cal on off pairs
590 //
591 Vector<Float> specCalon, specCaloff;
592 // to store scalar (mean) tsys
593 Vector<Float> tsysout(1);
594 uInt tcalId, polno;
595 Double offint, onint;
596 outpolCol.get(i, polno);
597 outspecCol.get(i, specCaloff);
598 onspecCol.get(i, specCalon);
599 Vector<uChar> flagCaloff, flagCalon;
600 outflagCol.get(i, flagCaloff);
601 onflagCol.get(i, flagCalon);
602 outtcalIdCol.get(i, tcalId);
603 outintCol.get(i, offint);
604 onintCol.get(i, onint);
605 // caluculate mean Tsys
606 uInt nchan = specCaloff.nelements();
607 // percentage of edge cut off
608 uInt pc = 10;
609 uInt bchan = nchan/pc;
610 uInt echan = nchan-bchan;
611
612 Slicer chansl(IPosition(1,bchan-1), IPosition(1,echan-1), IPosition(1,1),Slicer::endIsLast);
613 Vector<Float> testsubsp = specCaloff(chansl);
614 MaskedArray<Float> spoff = maskedArray( specCaloff(chansl),flagCaloff(chansl) );
615 MaskedArray<Float> spon = maskedArray( specCalon(chansl),flagCalon(chansl) );
616 MaskedArray<Float> spdiff = spon-spoff;
617 uInt noff = spoff.nelementsValid();
618 //uInt non = spon.nelementsValid();
619 uInt ndiff = spdiff.nelementsValid();
620 Float meantsys;
621
622/**
623 Double subspec, subdiff;
624 uInt usednchan;
625 subspec = 0;
626 subdiff = 0;
627 usednchan = 0;
628 for(uInt k=(bchan-1); k<echan; k++) {
629 subspec += specCaloff[k];
630 subdiff += static_cast<Double>(specCalon[k]-specCaloff[k]);
631 ++usednchan;
632 }
633**/
634 // get tcal if input tcal <= 0
635 String tcalt;
636 Float tcalUsed;
637 tcalUsed = tcal;
638 if ( tcal <= 0.0 ) {
639 caloff->tcal().getEntry(tcalt, tcalout, tcalId);
640 if (polno<=3) {
641 tcalUsed = tcalout[polno];
642 }
643 else {
644 tcalUsed = tcalout[0];
645 }
646 }
647
648 Float meanoff;
649 Float meandiff;
650 if (noff && ndiff) {
651 //Debug
652 //if(noff!=ndiff) cerr<<"noff and ndiff is not equal"<<endl;
653 meanoff = sum(spoff)/noff;
654 meandiff = sum(spdiff)/ndiff;
655 meantsys= (meanoff/meandiff )*tcalUsed + tcalUsed/2;
656 }
657 else {
658 meantsys=1;
659 }
660
661 tsysout[0] = Float(meantsys);
662 MaskedArray<Float> mcaloff = maskedArray(specCaloff, flagCaloff);
663 MaskedArray<Float> mcalon = maskedArray(specCalon, flagCalon);
664 MaskedArray<Float> sig = Float(0.5) * (mcaloff + mcalon);
665 //uInt ncaloff = mcaloff.nelementsValid();
666 //uInt ncalon = mcalon.nelementsValid();
667
668 outintCol.put(i, offint+onint);
669 outspecCol.put(i, sig.getArray());
670 outflagCol.put(i, flagsFromMA(sig));
671 outtsysCol.put(i, tsysout);
672 }
673 ++sit;
674 ++s2it;
675 }
676 return out;
677}
678
679//dosigref - migrated from GBT IDL's dosigref.pro, do calibration of position switch
680// observatiions.
681// input: sig and ref scantables, and an optional boxcar smoothing width(default width=0,
682// no smoothing).
683// output: resultant scantable [= (sig-ref/ref)*tsys]
684CountedPtr< Scantable > STMath::dosigref( const CountedPtr < Scantable >& sig,
685 const CountedPtr < Scantable >& ref,
686 int smoothref,
687 casa::Float tsysv,
688 casa::Float tau )
689{
690if ( ! ref->conformant(*sig) ) {
691 throw(AipsError("'sig' and 'ref' scantables are not conformant."));
692 }
693 setInsitu(false);
694 CountedPtr< Scantable > out = getScantable(sig, false);
695 CountedPtr< Scantable > smref;
696 if ( smoothref > 1 ) {
697 float fsmoothref = static_cast<float>(smoothref);
698 std::string inkernel = "boxcar";
699 smref = smooth(ref, inkernel, fsmoothref );
700 ostringstream oss;
701 oss<<"Applied smoothing of "<<fsmoothref<<" on the reference."<<endl;
702 pushLog(String(oss));
703 }
704 else {
705 smref = ref;
706 }
707 Table& tout = out->table();
708 const Table& tref = smref->table();
709 if ( tout.nrow() != tref.nrow() ) {
710 throw(AipsError("Mismatch in number of rows to form on-source and reference pair."));
711 }
712 // iteration by scanno? or cycle no.
713 TableIterator sit(tout, "SCANNO");
714 TableIterator s2it(tref, "SCANNO");
715 while ( !sit.pastEnd() ) {
716 Table ton = sit.table();
717 Table t = s2it.table();
718 ScalarColumn<Double> outintCol(ton, "INTERVAL");
719 ArrayColumn<Float> outspecCol(ton, "SPECTRA");
720 ArrayColumn<Float> outtsysCol(ton, "TSYS");
721 ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
722 ArrayColumn<Float> refspecCol(t, "SPECTRA");
723 ROScalarColumn<Double> refintCol(t, "INTERVAL");
724 ROArrayColumn<Float> reftsysCol(t, "TSYS");
725 ArrayColumn<uChar> refflagCol(t, "FLAGTRA");
726 ROScalarColumn<Float> refelevCol(t, "ELEVATION");
727 for (uInt i=0; i < ton.nrow(); ++i) {
728
729 Double onint, refint;
730 Vector<Float> specon, specref;
731 // to store scalar (mean) tsys
732 Vector<Float> tsysref;
733 outintCol.get(i, onint);
734 refintCol.get(i, refint);
735 outspecCol.get(i, specon);
736 refspecCol.get(i, specref);
737 Vector<uChar> flagref, flagon;
738 outflagCol.get(i, flagon);
739 refflagCol.get(i, flagref);
740 reftsysCol.get(i, tsysref);
741
742 Float tsysrefscalar;
743 if ( tsysv > 0.0 ) {
744 ostringstream oss;
745 Float elev;
746 refelevCol.get(i, elev);
747 oss << "user specified Tsys = " << tsysv;
748 // do recalc elevation if EL = 0
749 if ( elev == 0 ) {
750 throw(AipsError("EL=0, elevation data is missing."));
751 } else {
752 if ( tau <= 0.0 ) {
753 throw(AipsError("Valid tau is not supplied."));
754 } else {
755 tsysrefscalar = tsysv * exp(tau/elev);
756 }
757 }
758 oss << ", corrected (for El) tsys= "<<tsysrefscalar;
759 pushLog(String(oss));
760 }
761 else {
762 tsysrefscalar = tsysref[0];
763 }
764 //get quotient spectrum
765 MaskedArray<Float> mref = maskedArray(specref, flagref);
766 MaskedArray<Float> mon = maskedArray(specon, flagon);
767 MaskedArray<Float> specres = tsysrefscalar*((mon - mref)/mref);
768 Double resint = onint*refint*smoothref/(onint+refint*smoothref);
769
770 //Debug
771 //cerr<<"Tsys used="<<tsysrefscalar<<endl;
772 // fill the result, replay signal tsys by reference tsys
773 outintCol.put(i, resint);
774 outspecCol.put(i, specres.getArray());
775 outflagCol.put(i, flagsFromMA(specres));
776 outtsysCol.put(i, tsysref);
777 }
778 ++sit;
779 ++s2it;
780 }
781 return out;
782}
783
784CountedPtr< Scantable > STMath::donod(const casa::CountedPtr<Scantable>& s,
785 const std::vector<int>& scans,
786 int smoothref,
787 casa::Float tsysv,
788 casa::Float tau,
789 casa::Float tcal )
790
791{
792 setInsitu(false);
793 STSelector sel;
794 std::vector<int> scan1, scan2, beams;
795 std::vector< vector<int> > scanpair;
796 std::vector<string> calstate;
797 String msg;
798
799 CountedPtr< Scantable > s1b1on, s1b1off, s1b2on, s1b2off;
800 CountedPtr< Scantable > s2b1on, s2b1off, s2b2on, s2b2off;
801
802 std::vector< CountedPtr< Scantable > > sctables;
803 sctables.push_back(s1b1on);
804 sctables.push_back(s1b1off);
805 sctables.push_back(s1b2on);
806 sctables.push_back(s1b2off);
807 sctables.push_back(s2b1on);
808 sctables.push_back(s2b1off);
809 sctables.push_back(s2b2on);
810 sctables.push_back(s2b2off);
811
812 //check scanlist
813 int n=s->checkScanInfo(scans);
814 if (n==1) {
815 throw(AipsError("Incorrect scan pairs. "));
816 }
817
818 // Assume scans contain only a pair of consecutive scan numbers.
819 // It is assumed that first beam, b1, is on target.
820 // There is no check if the first beam is on or not.
821 if ( scans.size()==1 ) {
822 scan1.push_back(scans[0]);
823 scan2.push_back(scans[0]+1);
824 } else if ( scans.size()==2 ) {
825 scan1.push_back(scans[0]);
826 scan2.push_back(scans[1]);
827 } else {
828 if ( scans.size()%2 == 0 ) {
829 for (uInt i=0; i<scans.size(); i++) {
830 if (i%2 == 0) {
831 scan1.push_back(scans[i]);
832 }
833 else {
834 scan2.push_back(scans[i]);
835 }
836 }
837 } else {
838 throw(AipsError("Odd numbers of scans, cannot form pairs."));
839 }
840 }
841 scanpair.push_back(scan1);
842 scanpair.push_back(scan2);
843 calstate.push_back("*calon");
844 calstate.push_back("*[^calon]");
845 CountedPtr< Scantable > ws = getScantable(s, false);
846 uInt l=0;
847 while ( l < sctables.size() ) {
848 for (uInt i=0; i < 2; i++) {
849 for (uInt j=0; j < 2; j++) {
850 for (uInt k=0; k < 2; k++) {
851 sel.reset();
852 sel.setScans(scanpair[i]);
853 sel.setName(calstate[k]);
854 beams.clear();
855 beams.push_back(j);
856 sel.setBeams(beams);
857 ws->setSelection(sel);
858 sctables[l]= getScantable(ws, false);
859 l++;
860 }
861 }
862 }
863 }
864
865 // replace here by splitData or getData functionality
866 CountedPtr< Scantable > sig1;
867 CountedPtr< Scantable > ref1;
868 CountedPtr< Scantable > sig2;
869 CountedPtr< Scantable > ref2;
870 CountedPtr< Scantable > calb1;
871 CountedPtr< Scantable > calb2;
872
873 msg=String("Processing dototalpower for subset of the data");
874 ostringstream oss1;
875 oss1 << msg << endl;
876 pushLog(String(oss1));
877 // Debug for IRC CS data
878 //float tcal1=7.0;
879 //float tcal2=4.0;
880 sig1 = dototalpower(sctables[0], sctables[1], tcal=tcal);
881 ref1 = dototalpower(sctables[2], sctables[3], tcal=tcal);
882 ref2 = dototalpower(sctables[4], sctables[5], tcal=tcal);
883 sig2 = dototalpower(sctables[6], sctables[7], tcal=tcal);
884
885 // correction of user-specified tsys for elevation here
886
887 // dosigref calibration
888 msg=String("Processing dosigref for subset of the data");
889 ostringstream oss2;
890 oss2 << msg << endl;
891 pushLog(String(oss2));
892 calb1=dosigref(sig1,ref2,smoothref,tsysv,tau);
893 calb2=dosigref(sig2,ref1,smoothref,tsysv,tau);
894
895 // iteration by scanno or cycle no.
896 Table& tcalb1 = calb1->table();
897 Table& tcalb2 = calb2->table();
898 TableIterator sit(tcalb1, "SCANNO");
899 TableIterator s2it(tcalb2, "SCANNO");
900 while ( !sit.pastEnd() ) {
901 Table t1 = sit.table();
902 Table t2= s2it.table();
903 ArrayColumn<Float> outspecCol(t1, "SPECTRA");
904 ArrayColumn<Float> outtsysCol(t1, "TSYS");
905 ArrayColumn<uChar> outflagCol(t1, "FLAGTRA");
906 ScalarColumn<Double> outintCol(t1, "INTERVAL");
907 ArrayColumn<Float> t2specCol(t2, "SPECTRA");
908 ROArrayColumn<Float> t2tsysCol(t2, "TSYS");
909 ArrayColumn<uChar> t2flagCol(t2, "FLAGTRA");
910 ROScalarColumn<Double> t2intCol(t2, "INTERVAL");
911 for (uInt i=0; i < t1.nrow(); ++i) {
912 Vector<Float> spec1, spec2;
913 // to store scalar (mean) tsys
914 Vector<Float> tsys1, tsys2;
915 Vector<uChar> flag1, flag2;
916 Double tint1, tint2;
917 outspecCol.get(i, spec1);
918 t2specCol.get(i, spec2);
919 outflagCol.get(i, flag1);
920 t2flagCol.get(i, flag2);
921 outtsysCol.get(i, tsys1);
922 t2tsysCol.get(i, tsys2);
923 outintCol.get(i, tint1);
924 t2intCol.get(i, tint2);
925 // average
926 // assume scalar tsys for weights
927 Float wt1, wt2, tsyssq1, tsyssq2;
928 tsyssq1 = tsys1[0]*tsys1[0];
929 tsyssq2 = tsys2[0]*tsys2[0];
930 wt1 = Float(tint1)/tsyssq1;
931 wt2 = Float(tint2)/tsyssq2;
932 Float invsumwt=1/(wt1+wt2);
933 MaskedArray<Float> mspec1 = maskedArray(spec1, flag1);
934 MaskedArray<Float> mspec2 = maskedArray(spec2, flag2);
935 MaskedArray<Float> avspec = invsumwt * (wt1*mspec1 + wt2*mspec2);
936 //Array<Float> avtsys = Float(0.5) * (tsys1 + tsys2);
937 // cerr<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<endl;
938 tsys1[0] = sqrt(tsyssq1 + tsyssq2);
939 Array<Float> avtsys = tsys1;
940
941 outspecCol.put(i, avspec.getArray());
942 outflagCol.put(i, flagsFromMA(avspec));
943 outtsysCol.put(i, avtsys);
944 }
945 ++sit;
946 ++s2it;
947 }
948 return calb1;
949}
950
951//GBTIDL version of frequency switched data calibration
952CountedPtr< Scantable > STMath::dofs( const CountedPtr< Scantable >& s,
953 const std::vector<int>& scans,
954 int smoothref,
955 casa::Float tsysv,
956 casa::Float tau,
957 casa::Float tcal )
958{
959
960
961 STSelector sel;
962 CountedPtr< Scantable > ws = getScantable(s, false);
963 CountedPtr< Scantable > sig, sigwcal, ref, refwcal;
964 CountedPtr< Scantable > calsig, calref, out, out1, out2;
965 Bool nofold=False;
966
967 //split the data
968 sel.setName("*_fs");
969 ws->setSelection(sel);
970 sig = getScantable(ws,false);
971 sel.reset();
972 sel.setName("*_fs_calon");
973 ws->setSelection(sel);
974 sigwcal = getScantable(ws,false);
975 sel.reset();
976 sel.setName("*_fsr");
977 ws->setSelection(sel);
978 ref = getScantable(ws,false);
979 sel.reset();
980 sel.setName("*_fsr_calon");
981 ws->setSelection(sel);
982 refwcal = getScantable(ws,false);
983
984 calsig = dototalpower(sigwcal, sig, tcal=tcal);
985 calref = dototalpower(refwcal, ref, tcal=tcal);
986
987 out1=dosigref(calsig,calref,smoothref,tsysv,tau);
988 out2=dosigref(calref,calsig,smoothref,tsysv,tau);
989
990 Table& tabout1=out1->table();
991 Table& tabout2=out2->table();
992 ROScalarColumn<uInt> freqidCol1(tabout1, "FREQ_ID");
993 ROScalarColumn<uInt> freqidCol2(tabout2, "FREQ_ID");
994 ROArrayColumn<Float> specCol(tabout2, "SPECTRA");
995 Vector<Float> spec; specCol.get(0, spec);
996 uInt nchan = spec.nelements();
997 uInt freqid1; freqidCol1.get(0,freqid1);
998 uInt freqid2; freqidCol2.get(0,freqid2);
999 Double rp1, rp2, rv1, rv2, inc1, inc2;
1000 out1->frequencies().getEntry(rp1, rv1, inc1, freqid1);
1001 out2->frequencies().getEntry(rp2, rv2, inc2, freqid2);
1002 if (rp1==rp2) {
1003 Double foffset = rv1 - rv2;
1004 uInt choffset = static_cast<uInt>(foffset/abs(inc2));
1005 if (choffset >= nchan) {
1006 cerr<<"out-band frequency switching, no folding"<<endl;
1007 nofold = True;
1008 }
1009 }
1010
1011 if (nofold) {
1012 std::vector< CountedPtr< Scantable > > tabs;
1013 tabs.push_back(out1);
1014 tabs.push_back(out2);
1015 out = merge(tabs);
1016 }
1017 else { //folding is not implemented yet
1018 out = out1;
1019 }
1020
1021 return out;
1022}
1023
1024CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
1025{
1026 // make copy or reference
1027 CountedPtr< Scantable > out = getScantable(in, false);
1028 Table& tout = out->table();
1029 Block<String> cols(4);
1030 cols[0] = String("SCANNO");
1031 cols[1] = String("CYCLENO");
1032 cols[2] = String("BEAMNO");
1033 cols[3] = String("POLNO");
1034 TableIterator iter(tout, cols);
1035 while (!iter.pastEnd()) {
1036 Table subt = iter.table();
1037 // this should leave us with two rows for the two IFs....if not ignore
1038 if (subt.nrow() != 2 ) {
1039 continue;
1040 }
1041 ArrayColumn<Float> specCol(subt, "SPECTRA");
1042 ArrayColumn<Float> tsysCol(subt, "TSYS");
1043 ArrayColumn<uChar> flagCol(subt, "FLAGTRA");
1044 Vector<Float> onspec,offspec, ontsys, offtsys;
1045 Vector<uChar> onflag, offflag;
1046 tsysCol.get(0, ontsys); tsysCol.get(1, offtsys);
1047 specCol.get(0, onspec); specCol.get(1, offspec);
1048 flagCol.get(0, onflag); flagCol.get(1, offflag);
1049 MaskedArray<Float> on = maskedArray(onspec, onflag);
1050 MaskedArray<Float> off = maskedArray(offspec, offflag);
1051 MaskedArray<Float> oncopy = on.copy();
1052
1053 on /= off; on -= 1.0f;
1054 on *= ontsys[0];
1055 off /= oncopy; off -= 1.0f;
1056 off *= offtsys[0];
1057 specCol.put(0, on.getArray());
1058 const Vector<Bool>& m0 = on.getMask();
1059 Vector<uChar> flags0(m0.shape());
1060 convertArray(flags0, !m0);
1061 flagCol.put(0, flags0);
1062
1063 specCol.put(1, off.getArray());
1064 const Vector<Bool>& m1 = off.getMask();
1065 Vector<uChar> flags1(m1.shape());
1066 convertArray(flags1, !m1);
1067 flagCol.put(1, flags1);
1068 ++iter;
1069 }
1070
1071 return out;
1072}
1073
1074std::vector< float > STMath::statistic( const CountedPtr< Scantable > & in,
1075 const std::vector< bool > & mask,
1076 const std::string& which )
1077{
1078
1079 Vector<Bool> m(mask);
1080 const Table& tab = in->table();
1081 ROArrayColumn<Float> specCol(tab, "SPECTRA");
1082 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1083 std::vector<float> out;
1084 for (uInt i=0; i < tab.nrow(); ++i ) {
1085 Vector<Float> spec; specCol.get(i, spec);
1086 Vector<uChar> flag; flagCol.get(i, flag);
1087 MaskedArray<Float> ma = maskedArray(spec, flag);
1088 float outstat = 0.0;
1089 if ( spec.nelements() == m.nelements() ) {
1090 outstat = mathutil::statistics(which, ma(m));
1091 } else {
1092 outstat = mathutil::statistics(which, ma);
1093 }
1094 out.push_back(outstat);
1095 }
1096 return out;
1097}
1098
1099std::vector< int > STMath::minMaxPos( const CountedPtr< Scantable > & in,
1100 const std::vector< bool > & mask,
1101 const std::string& which )
1102{
1103
1104 Vector<Bool> m(mask);
1105 const Table& tab = in->table();
1106 ROArrayColumn<Float> specCol(tab, "SPECTRA");
1107 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1108 std::vector<int> out;
1109 for (uInt i=0; i < tab.nrow(); ++i ) {
1110 Vector<Float> spec; specCol.get(i, spec);
1111 Vector<uChar> flag; flagCol.get(i, flag);
1112 MaskedArray<Float> ma = maskedArray(spec, flag);
1113 if (ma.ndim() != 1) {
1114 throw (ArrayError(
1115 "std::vector<int> STMath::minMaxPos("
1116 "ContedPtr<Scantable> &in, std::vector<bool> &mask, "
1117 " std::string &which)"
1118 " - MaskedArray is not 1D"));
1119 }
1120 IPosition outpos(1,0);
1121 if ( spec.nelements() == m.nelements() ) {
1122 outpos = mathutil::minMaxPos(which, ma(m));
1123 } else {
1124 outpos = mathutil::minMaxPos(which, ma);
1125 }
1126 out.push_back(outpos[0]);
1127 }
1128 return out;
1129}
1130
1131CountedPtr< Scantable > STMath::bin( const CountedPtr< Scantable > & in,
1132 int width )
1133{
1134 if ( !in->getSelection().empty() ) throw(AipsError("Can't bin subset of the data."));
1135 CountedPtr< Scantable > out = getScantable(in, false);
1136 Table& tout = out->table();
1137 out->frequencies().rescale(width, "BIN");
1138 ArrayColumn<Float> specCol(tout, "SPECTRA");
1139 ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
1140 for (uInt i=0; i < tout.nrow(); ++i ) {
1141 MaskedArray<Float> main = maskedArray(specCol(i), flagCol(i));
1142 MaskedArray<Float> maout;
1143 LatticeUtilities::bin(maout, main, 0, Int(width));
1144 /// @todo implement channel based tsys binning
1145 specCol.put(i, maout.getArray());
1146 flagCol.put(i, flagsFromMA(maout));
1147 // take only the first binned spectrum's length for the deprecated
1148 // global header item nChan
1149 if (i==0) tout.rwKeywordSet().define(String("nChan"),
1150 Int(maout.getArray().nelements()));
1151 }
1152 return out;
1153}
1154
1155CountedPtr< Scantable > STMath::resample( const CountedPtr< Scantable >& in,
1156 const std::string& method,
1157 float width )
1158//
1159// Should add the possibility of width being specified in km/s. This means
1160// that for each freqID (SpectralCoordinate) we will need to convert to an
1161// average channel width (say at the reference pixel). Then we would need
1162// to be careful to make sure each spectrum (of different freqID)
1163// is the same length.
1164//
1165{
1166 //InterpolateArray1D<Double,Float>::InterpolationMethod interp;
1167 Int interpMethod(stringToIMethod(method));
1168
1169 CountedPtr< Scantable > out = getScantable(in, false);
1170 Table& tout = out->table();
1171
1172// Resample SpectralCoordinates (one per freqID)
1173 out->frequencies().rescale(width, "RESAMPLE");
1174 TableIterator iter(tout, "IFNO");
1175 TableRow row(tout);
1176 while ( !iter.pastEnd() ) {
1177 Table tab = iter.table();
1178 ArrayColumn<Float> specCol(tab, "SPECTRA");
1179 //ArrayColumn<Float> tsysCol(tout, "TSYS");
1180 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1181 Vector<Float> spec;
1182 Vector<uChar> flag;
1183 specCol.get(0,spec); // the number of channels should be constant per IF
1184 uInt nChanIn = spec.nelements();
1185 Vector<Float> xIn(nChanIn); indgen(xIn);
1186 Int fac = Int(nChanIn/width);
1187 Vector<Float> xOut(fac+10); // 10 to be safe - resize later
1188 uInt k = 0;
1189 Float x = 0.0;
1190 while (x < Float(nChanIn) ) {
1191 xOut(k) = x;
1192 k++;
1193 x += width;
1194 }
1195 uInt nChanOut = k;
1196 xOut.resize(nChanOut, True);
1197 // process all rows for this IFNO
1198 Vector<Float> specOut;
1199 Vector<Bool> maskOut;
1200 Vector<uChar> flagOut;
1201 for (uInt i=0; i < tab.nrow(); ++i) {
1202 specCol.get(i, spec);
1203 flagCol.get(i, flag);
1204 Vector<Bool> mask(flag.nelements());
1205 convertArray(mask, flag);
1206
1207 IPosition shapeIn(spec.shape());
1208 //sh.nchan = nChanOut;
1209 InterpolateArray1D<Float,Float>::interpolate(specOut, maskOut, xOut,
1210 xIn, spec, mask,
1211 interpMethod, True, True);
1212 /// @todo do the same for channel based Tsys
1213 flagOut.resize(maskOut.nelements());
1214 convertArray(flagOut, maskOut);
1215 specCol.put(i, specOut);
1216 flagCol.put(i, flagOut);
1217 }
1218 ++iter;
1219 }
1220
1221 return out;
1222}
1223
1224STMath::imethod STMath::stringToIMethod(const std::string& in)
1225{
1226 static STMath::imap lookup;
1227
1228 // initialize the lookup table if necessary
1229 if ( lookup.empty() ) {
1230 lookup["nearest"] = InterpolateArray1D<Double,Float>::nearestNeighbour;
1231 lookup["linear"] = InterpolateArray1D<Double,Float>::linear;
1232 lookup["cubic"] = InterpolateArray1D<Double,Float>::cubic;
1233 lookup["spline"] = InterpolateArray1D<Double,Float>::spline;
1234 }
1235
1236 STMath::imap::const_iterator iter = lookup.find(in);
1237
1238 if ( lookup.end() == iter ) {
1239 std::string message = in;
1240 message += " is not a valid interpolation mode";
1241 throw(AipsError(message));
1242 }
1243 return iter->second;
1244}
1245
1246WeightType STMath::stringToWeight(const std::string& in)
1247{
1248 static std::map<std::string, WeightType> lookup;
1249
1250 // initialize the lookup table if necessary
1251 if ( lookup.empty() ) {
1252 lookup["NONE"] = asap::NONE;
1253 lookup["TINT"] = asap::TINT;
1254 lookup["TINTSYS"] = asap::TINTSYS;
1255 lookup["TSYS"] = asap::TSYS;
1256 lookup["VAR"] = asap::VAR;
1257 }
1258
1259 std::map<std::string, WeightType>::const_iterator iter = lookup.find(in);
1260
1261 if ( lookup.end() == iter ) {
1262 std::string message = in;
1263 message += " is not a valid weighting mode";
1264 throw(AipsError(message));
1265 }
1266 return iter->second;
1267}
1268
1269CountedPtr< Scantable > STMath::gainElevation( const CountedPtr< Scantable >& in,
1270 const vector< float > & coeff,
1271 const std::string & filename,
1272 const std::string& method)
1273{
1274 // Get elevation data from Scantable and convert to degrees
1275 CountedPtr< Scantable > out = getScantable(in, false);
1276 Table& tab = out->table();
1277 ROScalarColumn<Float> elev(tab, "ELEVATION");
1278 Vector<Float> x = elev.getColumn();
1279 x *= Float(180 / C::pi); // Degrees
1280
1281 Vector<Float> coeffs(coeff);
1282 const uInt nc = coeffs.nelements();
1283 if ( filename.length() > 0 && nc > 0 ) {
1284 throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
1285 }
1286
1287 // Correct
1288 if ( nc > 0 || filename.length() == 0 ) {
1289 // Find instrument
1290 Bool throwit = True;
1291 Instrument inst =
1292 STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
1293 throwit);
1294
1295 // Set polynomial
1296 Polynomial<Float>* ppoly = 0;
1297 Vector<Float> coeff;
1298 String msg;
1299 if ( nc > 0 ) {
1300 ppoly = new Polynomial<Float>(nc);
1301 coeff = coeffs;
1302 msg = String("user");
1303 } else {
1304 STAttr sdAttr;
1305 coeff = sdAttr.gainElevationPoly(inst);
1306 ppoly = new Polynomial<Float>(3);
1307 msg = String("built in");
1308 }
1309
1310 if ( coeff.nelements() > 0 ) {
1311 ppoly->setCoefficients(coeff);
1312 } else {
1313 delete ppoly;
1314 throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
1315 }
1316 ostringstream oss;
1317 oss << "Making polynomial correction with " << msg << " coefficients:" << endl;
1318 oss << " " << coeff;
1319 pushLog(String(oss));
1320 const uInt nrow = tab.nrow();
1321 Vector<Float> factor(nrow);
1322 for ( uInt i=0; i < nrow; ++i ) {
1323 factor[i] = 1.0 / (*ppoly)(x[i]);
1324 }
1325 delete ppoly;
1326 scaleByVector(tab, factor, true);
1327
1328 } else {
1329 // Read and correct
1330 pushLog("Making correction from ascii Table");
1331 scaleFromAsciiTable(tab, filename, method, x, true);
1332 }
1333 return out;
1334}
1335
1336void STMath::scaleFromAsciiTable(Table& in, const std::string& filename,
1337 const std::string& method,
1338 const Vector<Float>& xout, bool dotsys)
1339{
1340
1341// Read gain-elevation ascii file data into a Table.
1342
1343 String formatString;
1344 Table tbl = readAsciiTable(formatString, Table::Memory, filename, "", "", False);
1345 scaleFromTable(in, tbl, method, xout, dotsys);
1346}
1347
1348void STMath::scaleFromTable(Table& in,
1349 const Table& table,
1350 const std::string& method,
1351 const Vector<Float>& xout, bool dotsys)
1352{
1353
1354 ROScalarColumn<Float> geElCol(table, "ELEVATION");
1355 ROScalarColumn<Float> geFacCol(table, "FACTOR");
1356 Vector<Float> xin = geElCol.getColumn();
1357 Vector<Float> yin = geFacCol.getColumn();
1358 Vector<Bool> maskin(xin.nelements(),True);
1359
1360 // Interpolate (and extrapolate) with desired method
1361
1362 InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
1363
1364 Vector<Float> yout;
1365 Vector<Bool> maskout;
1366 InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout,
1367 xin, yin, maskin, interp,
1368 True, True);
1369
1370 scaleByVector(in, Float(1.0)/yout, dotsys);
1371}
1372
1373void STMath::scaleByVector( Table& in,
1374 const Vector< Float >& factor,
1375 bool dotsys )
1376{
1377 uInt nrow = in.nrow();
1378 if ( factor.nelements() != nrow ) {
1379 throw(AipsError("factors.nelements() != table.nelements()"));
1380 }
1381 ArrayColumn<Float> specCol(in, "SPECTRA");
1382 ArrayColumn<uChar> flagCol(in, "FLAGTRA");
1383 ArrayColumn<Float> tsysCol(in, "TSYS");
1384 for (uInt i=0; i < nrow; ++i) {
1385 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
1386 ma *= factor[i];
1387 specCol.put(i, ma.getArray());
1388 flagCol.put(i, flagsFromMA(ma));
1389 if ( dotsys ) {
1390 Vector<Float> tsys = tsysCol(i);
1391 tsys *= factor[i];
1392 tsysCol.put(i,tsys);
1393 }
1394 }
1395}
1396
1397CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in,
1398 float d, float etaap,
1399 float jyperk )
1400{
1401 CountedPtr< Scantable > out = getScantable(in, false);
1402 Table& tab = in->table();
1403 Unit fluxUnit(tab.keywordSet().asString("FluxUnit"));
1404 Unit K(String("K"));
1405 Unit JY(String("Jy"));
1406
1407 bool tokelvin = true;
1408 Double cfac = 1.0;
1409
1410 if ( fluxUnit == JY ) {
1411 pushLog("Converting to K");
1412 Quantum<Double> t(1.0,fluxUnit);
1413 Quantum<Double> t2 = t.get(JY);
1414 cfac = (t2 / t).getValue(); // value to Jy
1415
1416 tokelvin = true;
1417 out->setFluxUnit("K");
1418 } else if ( fluxUnit == K ) {
1419 pushLog("Converting to Jy");
1420 Quantum<Double> t(1.0,fluxUnit);
1421 Quantum<Double> t2 = t.get(K);
1422 cfac = (t2 / t).getValue(); // value to K
1423
1424 tokelvin = false;
1425 out->setFluxUnit("Jy");
1426 } else {
1427 throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
1428 }
1429 // Make sure input values are converted to either Jy or K first...
1430 Float factor = cfac;
1431
1432 // Select method
1433 if (jyperk > 0.0) {
1434 factor *= jyperk;
1435 if ( tokelvin ) factor = 1.0 / jyperk;
1436 ostringstream oss;
1437 oss << "Jy/K = " << jyperk;
1438 pushLog(String(oss));
1439 Vector<Float> factors(tab.nrow(), factor);
1440 scaleByVector(tab,factors, false);
1441 } else if ( etaap > 0.0) {
1442 if (d < 0) {
1443 Instrument inst =
1444 STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
1445 True);
1446 STAttr sda;
1447 d = sda.diameter(inst);
1448 }
1449 jyperk = STAttr::findJyPerK(etaap, d);
1450 ostringstream oss;
1451 oss << "Jy/K = " << jyperk;
1452 pushLog(String(oss));
1453 factor *= jyperk;
1454 if ( tokelvin ) {
1455 factor = 1.0 / factor;
1456 }
1457 Vector<Float> factors(tab.nrow(), factor);
1458 scaleByVector(tab, factors, False);
1459 } else {
1460
1461 // OK now we must deal with automatic look up of values.
1462 // We must also deal with the fact that the factors need
1463 // to be computed per IF and may be different and may
1464 // change per integration.
1465
1466 pushLog("Looking up conversion factors");
1467 convertBrightnessUnits(out, tokelvin, cfac);
1468 }
1469
1470 return out;
1471}
1472
1473void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in,
1474 bool tokelvin, float cfac )
1475{
1476 Table& table = in->table();
1477 Instrument inst =
1478 STAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True);
1479 TableIterator iter(table, "FREQ_ID");
1480 STFrequencies stfreqs = in->frequencies();
1481 STAttr sdAtt;
1482 while (!iter.pastEnd()) {
1483 Table tab = iter.table();
1484 ArrayColumn<Float> specCol(tab, "SPECTRA");
1485 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1486 ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID");
1487 MEpoch::ROScalarColumn timeCol(tab, "TIME");
1488
1489 uInt freqid; freqidCol.get(0, freqid);
1490 Vector<Float> tmpspec; specCol.get(0, tmpspec);
1491 // STAttr.JyPerK has a Vector interface... change sometime.
1492 Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements()));
1493 for ( uInt i=0; i<tab.nrow(); ++i) {
1494 Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0];
1495 Float factor = cfac * jyperk;
1496 if ( tokelvin ) factor = Float(1.0) / factor;
1497 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
1498 ma *= factor;
1499 specCol.put(i, ma.getArray());
1500 flagCol.put(i, flagsFromMA(ma));
1501 }
1502 ++iter;
1503 }
1504}
1505
1506CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
1507 float tau )
1508{
1509 CountedPtr< Scantable > out = getScantable(in, false);
1510
1511 Table tab = out->table();
1512 ROScalarColumn<Float> elev(tab, "ELEVATION");
1513 ArrayColumn<Float> specCol(tab, "SPECTRA");
1514 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1515 for ( uInt i=0; i<tab.nrow(); ++i) {
1516 Float zdist = Float(C::pi_2) - elev(i);
1517 Float factor = exp(tau/cos(zdist));
1518 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
1519 ma *= factor;
1520 specCol.put(i, ma.getArray());
1521 flagCol.put(i, flagsFromMA(ma));
1522 }
1523 return out;
1524}
1525
1526CountedPtr< Scantable > STMath::smoothOther( const CountedPtr< Scantable >& in,
1527 const std::string& kernel,
1528 float width )
1529{
1530 CountedPtr< Scantable > out = getScantable(in, false);
1531 Table& table = out->table();
1532 ArrayColumn<Float> specCol(table, "SPECTRA");
1533 ArrayColumn<uChar> flagCol(table, "FLAGTRA");
1534 Vector<Float> spec;
1535 Vector<uChar> flag;
1536 for ( uInt i=0; i<table.nrow(); ++i) {
1537 specCol.get(i, spec);
1538 flagCol.get(i, flag);
1539 Vector<Bool> mask(flag.nelements());
1540 convertArray(mask, flag);
1541 Vector<Float> specout;
1542 Vector<Bool> maskout;
1543 if ( kernel == "hanning" ) {
1544 mathutil::hanning(specout, maskout, spec , !mask);
1545 convertArray(flag, !maskout);
1546 } else if ( kernel == "rmedian" ) {
1547 mathutil::runningMedian(specout, maskout, spec , mask, width);
1548 convertArray(flag, maskout);
1549 }
1550 flagCol.put(i, flag);
1551 specCol.put(i, specout);
1552 }
1553 return out;
1554}
1555
1556CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
1557 const std::string& kernel, float width )
1558{
1559 if (kernel == "rmedian" || kernel == "hanning") {
1560 return smoothOther(in, kernel, width);
1561 }
1562 CountedPtr< Scantable > out = getScantable(in, false);
1563 Table& table = out->table();
1564 VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
1565 // same IFNO should have same no of channels
1566 // this saves overhead
1567 TableIterator iter(table, "IFNO");
1568 while (!iter.pastEnd()) {
1569 Table tab = iter.table();
1570 ArrayColumn<Float> specCol(tab, "SPECTRA");
1571 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1572 Vector<Float> tmpspec; specCol.get(0, tmpspec);
1573 uInt nchan = tmpspec.nelements();
1574 Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
1575 Convolver<Float> conv(kvec, IPosition(1,nchan));
1576 Vector<Float> spec;
1577 Vector<uChar> flag;
1578 for ( uInt i=0; i<tab.nrow(); ++i) {
1579 specCol.get(i, spec);
1580 flagCol.get(i, flag);
1581 Vector<Bool> mask(flag.nelements());
1582 convertArray(mask, flag);
1583 Vector<Float> specout;
1584 mathutil::replaceMaskByZero(specout, mask);
1585 conv.linearConv(specout, spec);
1586 specCol.put(i, specout);
1587 }
1588 ++iter;
1589 }
1590 return out;
1591}
1592
1593CountedPtr< Scantable >
1594 STMath::merge( const std::vector< CountedPtr < Scantable > >& in )
1595{
1596 if ( in.size() < 2 ) {
1597 throw(AipsError("Need at least two scantables to perform a merge."));
1598 }
1599 std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
1600 bool insitu = insitu_;
1601 setInsitu(false);
1602 CountedPtr< Scantable > out = getScantable(*it, false);
1603 setInsitu(insitu);
1604 Table& tout = out->table();
1605 ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
1606 ScalarColumn<uInt> scannocol(tout,"SCANNO"), focusidcol(tout,"FOCUS_ID");
1607 // Renumber SCANNO to be 0-based
1608 Vector<uInt> scannos = scannocol.getColumn();
1609 uInt offset = min(scannos);
1610 scannos -= offset;
1611 scannocol.putColumn(scannos);
1612 uInt newscanno = max(scannos)+1;
1613 ++it;
1614 while ( it != in.end() ){
1615 if ( ! (*it)->conformant(*out) ) {
1616 // log message: "ignoring scantable i, as it isn't
1617 // conformant with the other(s)"
1618 cerr << "oh oh" << endl;
1619 ++it;
1620 continue;
1621 }
1622 out->appendToHistoryTable((*it)->history());
1623 const Table& tab = (*it)->table();
1624 TableIterator scanit(tab, "SCANNO");
1625 while (!scanit.pastEnd()) {
1626 TableIterator freqit(scanit.table(), "FREQ_ID");
1627 while ( !freqit.pastEnd() ) {
1628 Table thetab = freqit.table();
1629 uInt nrow = tout.nrow();
1630 //tout.addRow(thetab.nrow());
1631 TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
1632 ROTableRow row(thetab);
1633 for ( uInt i=0; i<thetab.nrow(); ++i) {
1634 uInt k = nrow+i;
1635 scannocol.put(k, newscanno);
1636 const TableRecord& rec = row.get(i);
1637 Double rv,rp,inc;
1638 (*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
1639 uInt id;
1640 id = out->frequencies().addEntry(rp, rv, inc);
1641 freqidcol.put(k,id);
1642 //String name,fname;Double rf;
1643 Vector<String> name,fname;Vector<Double> rf;
1644 (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
1645 id = out->molecules().addEntry(rf, name, fname);
1646 molidcol.put(k, id);
1647 Float frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
1648 (*it)->focus().getEntry(fax, ftan, frot, fhand,
1649 fmount,fuser, fxy, fxyp,
1650 rec.asuInt("FOCUS_ID"));
1651 id = out->focus().addEntry(fax, ftan, frot, fhand,
1652 fmount,fuser, fxy, fxyp);
1653 focusidcol.put(k, id);
1654 }
1655 ++freqit;
1656 }
1657 ++newscanno;
1658 ++scanit;
1659 }
1660 ++it;
1661 }
1662 return out;
1663}
1664
1665CountedPtr< Scantable >
1666 STMath::invertPhase( const CountedPtr < Scantable >& in )
1667{
1668 return applyToPol(in, &STPol::invertPhase, Float(0.0));
1669}
1670
1671CountedPtr< Scantable >
1672 STMath::rotateXYPhase( const CountedPtr < Scantable >& in, float phase )
1673{
1674 return applyToPol(in, &STPol::rotatePhase, Float(phase));
1675}
1676
1677CountedPtr< Scantable >
1678 STMath::rotateLinPolPhase( const CountedPtr < Scantable >& in, float phase )
1679{
1680 return applyToPol(in, &STPol::rotateLinPolPhase, Float(phase));
1681}
1682
1683CountedPtr< Scantable > STMath::applyToPol( const CountedPtr<Scantable>& in,
1684 STPol::polOperation fptr,
1685 Float phase )
1686{
1687 CountedPtr< Scantable > out = getScantable(in, false);
1688 Table& tout = out->table();
1689 Block<String> cols(4);
1690 cols[0] = String("SCANNO");
1691 cols[1] = String("BEAMNO");
1692 cols[2] = String("IFNO");
1693 cols[3] = String("CYCLENO");
1694 TableIterator iter(tout, cols);
1695 CountedPtr<STPol> stpol = STPol::getPolClass(out->factories_,
1696 out->getPolType() );
1697 while (!iter.pastEnd()) {
1698 Table t = iter.table();
1699 ArrayColumn<Float> speccol(t, "SPECTRA");
1700 ScalarColumn<uInt> focidcol(t, "FOCUS_ID");
1701 ScalarColumn<Float> parancol(t, "PARANGLE");
1702 Matrix<Float> pols(speccol.getColumn());
1703 try {
1704 stpol->setSpectra(pols);
1705 Float fang,fhand,parang;
1706 fang = in->focusTable_.getTotalFeedAngle(focidcol(0));
1707 fhand = in->focusTable_.getFeedHand(focidcol(0));
1708 parang = parancol(0);
1709 /// @todo re-enable this
1710 // disable total feed angle to support paralactifying Caswell style
1711 stpol->setPhaseCorrections(parang, -parang, fhand);
1712 // use a member function pointer in STPol. This only works on
1713 // the STPol pointer itself, not the Counted Pointer so
1714 // derefernce it.
1715 (&(*(stpol))->*fptr)(phase);
1716 speccol.putColumn(stpol->getSpectra());
1717 } catch (AipsError& e) {
1718 //delete stpol;stpol=0;
1719 throw(e);
1720 }
1721 ++iter;
1722 }
1723 //delete stpol;stpol=0;
1724 return out;
1725}
1726
1727CountedPtr< Scantable >
1728 STMath::swapPolarisations( const CountedPtr< Scantable > & in )
1729{
1730 CountedPtr< Scantable > out = getScantable(in, false);
1731 Table& tout = out->table();
1732 Table t0 = tout(tout.col("POLNO") == 0);
1733 Table t1 = tout(tout.col("POLNO") == 1);
1734 if ( t0.nrow() != t1.nrow() )
1735 throw(AipsError("Inconsistent number of polarisations"));
1736 ArrayColumn<Float> speccol0(t0, "SPECTRA");
1737 ArrayColumn<uChar> flagcol0(t0, "FLAGTRA");
1738 ArrayColumn<Float> speccol1(t1, "SPECTRA");
1739 ArrayColumn<uChar> flagcol1(t1, "FLAGTRA");
1740 Matrix<Float> s0 = speccol0.getColumn();
1741 Matrix<uChar> f0 = flagcol0.getColumn();
1742 speccol0.putColumn(speccol1.getColumn());
1743 flagcol0.putColumn(flagcol1.getColumn());
1744 speccol1.putColumn(s0);
1745 flagcol1.putColumn(f0);
1746 return out;
1747}
1748
1749CountedPtr< Scantable >
1750 STMath::averagePolarisations( const CountedPtr< Scantable > & in,
1751 const std::vector<bool>& mask,
1752 const std::string& weight )
1753{
1754 if (in->npol() < 2 )
1755 throw(AipsError("averagePolarisations can only be applied to two or more"
1756 "polarisations"));
1757 bool insitu = insitu_;
1758 setInsitu(false);
1759 CountedPtr< Scantable > pols = getScantable(in, true);
1760 setInsitu(insitu);
1761 Table& tout = pols->table();
1762 std::string taql = "SELECT FROM $1 WHERE POLNO IN [0,1]";
1763 Table tab = tableCommand(taql, in->table());
1764 if (tab.nrow() == 0 )
1765 throw(AipsError("Could not find any rows with POLNO==0 and POLNO==1"));
1766 TableCopy::copyRows(tout, tab);
1767 TableVector<uInt> vec(tout, "POLNO");
1768 vec = 0;
1769 pols->table_.rwKeywordSet().define("nPol", Int(1));
1770 //pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
1771 pols->table_.rwKeywordSet().define("POLTYPE", in->getPolType());
1772 std::vector<CountedPtr<Scantable> > vpols;
1773 vpols.push_back(pols);
1774 CountedPtr< Scantable > out = average(vpols, mask, weight, "SCAN");
1775 return out;
1776}
1777
1778CountedPtr< Scantable >
1779 STMath::averageBeams( const CountedPtr< Scantable > & in,
1780 const std::vector<bool>& mask,
1781 const std::string& weight )
1782{
1783 bool insitu = insitu_;
1784 setInsitu(false);
1785 CountedPtr< Scantable > beams = getScantable(in, false);
1786 setInsitu(insitu);
1787 Table& tout = beams->table();
1788 // give all rows the same BEAMNO
1789 TableVector<uInt> vec(tout, "BEAMNO");
1790 vec = 0;
1791 beams->table_.rwKeywordSet().define("nBeam", Int(1));
1792 std::vector<CountedPtr<Scantable> > vbeams;
1793 vbeams.push_back(beams);
1794 CountedPtr< Scantable > out = average(vbeams, mask, weight, "SCAN");
1795 return out;
1796}
1797
1798
1799CountedPtr< Scantable >
1800 asap::STMath::frequencyAlign( const CountedPtr< Scantable > & in,
1801 const std::string & refTime,
1802 const std::string & method)
1803{
1804 // clone as this is not working insitu
1805 bool insitu = insitu_;
1806 setInsitu(false);
1807 CountedPtr< Scantable > out = getScantable(in, false);
1808 setInsitu(insitu);
1809 Table& tout = out->table();
1810 // Get reference Epoch to time of first row or given String
1811 Unit DAY(String("d"));
1812 MEpoch::Ref epochRef(in->getTimeReference());
1813 MEpoch refEpoch;
1814 if (refTime.length()>0) {
1815 Quantum<Double> qt;
1816 if (MVTime::read(qt,refTime)) {
1817 MVEpoch mv(qt);
1818 refEpoch = MEpoch(mv, epochRef);
1819 } else {
1820 throw(AipsError("Invalid format for Epoch string"));
1821 }
1822 } else {
1823 refEpoch = in->timeCol_(0);
1824 }
1825 MPosition refPos = in->getAntennaPosition();
1826
1827 InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
1828 // test if user frame is different to base frame
1829 if ( in->frequencies().getFrameString(true)
1830 == in->frequencies().getFrameString(false) ) {
1831 throw(AipsError("Can't convert as no output frame has been set"
1832 " (use set_freqframe) or it is aligned already."));
1833 }
1834 MFrequency::Types system = in->frequencies().getFrame();
1835 MVTime mvt(refEpoch.getValue());
1836 String epochout = mvt.string(MVTime::YMD) + String(" (") + refEpoch.getRefString() + String(")");
1837 ostringstream oss;
1838 oss << "Aligned at reference Epoch " << epochout
1839 << " in frame " << MFrequency::showType(system);
1840 pushLog(String(oss));
1841 // set up the iterator
1842 Block<String> cols(4);
1843 // select by constant direction
1844 cols[0] = String("SRCNAME");
1845 cols[1] = String("BEAMNO");
1846 // select by IF ( no of channels varies over this )
1847 cols[2] = String("IFNO");
1848 // select by restfrequency
1849 cols[3] = String("MOLECULE_ID");
1850 TableIterator iter(tout, cols);
1851 while ( !iter.pastEnd() ) {
1852 Table t = iter.table();
1853 MDirection::ROScalarColumn dirCol(t, "DIRECTION");
1854 TableIterator fiter(t, "FREQ_ID");
1855 // determine nchan from the first row. This should work as
1856 // we are iterating over BEAMNO and IFNO // we should have constant direction
1857
1858 ROArrayColumn<Float> sCol(t, "SPECTRA");
1859 MDirection direction = dirCol(0);
1860 uInt nchan = sCol(0).nelements();
1861 while ( !fiter.pastEnd() ) {
1862 Table ftab = fiter.table();
1863 ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
1864 // get the SpectralCoordinate for the freqid, which we are iterating over
1865 SpectralCoordinate sC = in->frequencies().getSpectralCoordinate(freqidCol(0));
1866 FrequencyAligner<Float> fa( sC, nchan, refEpoch,
1867 direction, refPos, system );
1868 // realign the SpectralCoordinate and put into the output Scantable
1869 Vector<String> units(1);
1870 units = String("Hz");
1871 Bool linear=True;
1872 SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
1873 sc2.setWorldAxisUnits(units);
1874 uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
1875 sc2.referenceValue()[0],
1876 sc2.increment()[0]);
1877 TableVector<uInt> tvec(ftab, "FREQ_ID");
1878 tvec = id;
1879 // create the "global" abcissa for alignment with same FREQ_ID
1880 Vector<Double> abc(nchan);
1881 Double w;
1882 for (uInt i=0; i<nchan; i++) {
1883 sC.toWorld(w,Double(i));
1884 abc[i] = w;
1885 }
1886 // cache abcissa for same time stamps, so iterate over those
1887 TableIterator timeiter(ftab, "TIME");
1888 while ( !timeiter.pastEnd() ) {
1889 Table tab = timeiter.table();
1890 ArrayColumn<Float> specCol(tab, "SPECTRA");
1891 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1892 MEpoch::ROScalarColumn timeCol(tab, "TIME");
1893 // use align abcissa cache after the first row
1894 bool first = true;
1895 // these rows should be just be POLNO
1896 for (int i=0; i<int(tab.nrow()); ++i) {
1897 // input values
1898 Vector<uChar> flag = flagCol(i);
1899 Vector<Bool> mask(flag.shape());
1900 Vector<Float> specOut, spec;
1901 spec = specCol(i);
1902 Vector<Bool> maskOut;Vector<uChar> flagOut;
1903 convertArray(mask, flag);
1904 // alignment
1905 Bool ok = fa.align(specOut, maskOut, abc, spec,
1906 mask, timeCol(i), !first,
1907 interp, False);
1908 // back into scantable
1909 flagOut.resize(maskOut.nelements());
1910 convertArray(flagOut, maskOut);
1911 flagCol.put(i, flagOut);
1912 specCol.put(i, specOut);
1913 // start abcissa caching
1914 first = false;
1915 }
1916 // next timestamp
1917 ++timeiter;
1918 }
1919 // next FREQ_ID
1920 ++fiter;
1921 }
1922 // next aligner
1923 ++iter;
1924 }
1925 // set this afterwards to ensure we are doing insitu correctly.
1926 out->frequencies().setFrame(system, true);
1927 return out;
1928}
1929
1930CountedPtr<Scantable>
1931 asap::STMath::convertPolarisation( const CountedPtr<Scantable>& in,
1932 const std::string & newtype )
1933{
1934 if (in->npol() != 2 && in->npol() != 4)
1935 throw(AipsError("Can only convert two or four polarisations."));
1936 if ( in->getPolType() == newtype )
1937 throw(AipsError("No need to convert."));
1938 if ( ! in->selector_.empty() )
1939 throw(AipsError("Can only convert whole scantable. Unset the selection."));
1940 bool insitu = insitu_;
1941 setInsitu(false);
1942 CountedPtr< Scantable > out = getScantable(in, true);
1943 setInsitu(insitu);
1944 Table& tout = out->table();
1945 tout.rwKeywordSet().define("POLTYPE", String(newtype));
1946
1947 Block<String> cols(4);
1948 cols[0] = "SCANNO";
1949 cols[1] = "CYCLENO";
1950 cols[2] = "BEAMNO";
1951 cols[3] = "IFNO";
1952 TableIterator it(in->originalTable_, cols);
1953 String basetype = in->getPolType();
1954 STPol* stpol = STPol::getPolClass(in->factories_, basetype);
1955 try {
1956 while ( !it.pastEnd() ) {
1957 Table tab = it.table();
1958 uInt row = tab.rowNumbers()[0];
1959 stpol->setSpectra(in->getPolMatrix(row));
1960 Float fang,fhand,parang;
1961 fang = in->focusTable_.getTotalFeedAngle(in->mfocusidCol_(row));
1962 fhand = in->focusTable_.getFeedHand(in->mfocusidCol_(row));
1963 parang = in->paraCol_(row);
1964 /// @todo re-enable this
1965 // disable total feed angle to support paralactifying Caswell style
1966 stpol->setPhaseCorrections(parang, -parang, fhand);
1967 Int npolout = 0;
1968 for (uInt i=0; i<tab.nrow(); ++i) {
1969 Vector<Float> outvec = stpol->getSpectrum(i, newtype);
1970 if ( outvec.nelements() > 0 ) {
1971 tout.addRow();
1972 TableCopy::copyRows(tout, tab, tout.nrow()-1, 0, 1);
1973 ArrayColumn<Float> sCol(tout,"SPECTRA");
1974 ScalarColumn<uInt> pCol(tout,"POLNO");
1975 sCol.put(tout.nrow()-1 ,outvec);
1976 pCol.put(tout.nrow()-1 ,uInt(npolout));
1977 npolout++;
1978 }
1979 }
1980 tout.rwKeywordSet().define("nPol", npolout);
1981 ++it;
1982 }
1983 } catch (AipsError& e) {
1984 delete stpol;
1985 throw(e);
1986 }
1987 delete stpol;
1988 return out;
1989}
1990
1991CountedPtr< Scantable >
1992 asap::STMath::mxExtract( const CountedPtr< Scantable > & in,
1993 const std::string & scantype )
1994{
1995 bool insitu = insitu_;
1996 setInsitu(false);
1997 CountedPtr< Scantable > out = getScantable(in, true);
1998 setInsitu(insitu);
1999 Table& tout = out->table();
2000 std::string taql = "SELECT FROM $1 WHERE BEAMNO != REFBEAMNO";
2001 if (scantype == "on") {
2002 taql = "SELECT FROM $1 WHERE BEAMNO == REFBEAMNO";
2003 }
2004 Table tab = tableCommand(taql, in->table());
2005 TableCopy::copyRows(tout, tab);
2006 if (scantype == "on") {
2007 // re-index SCANNO to 0
2008 TableVector<uInt> vec(tout, "SCANNO");
2009 vec = 0;
2010 }
2011 return out;
2012}
2013
2014CountedPtr< Scantable >
2015 asap::STMath::lagFlag( const CountedPtr< Scantable > & in,
2016 double frequency, double width )
2017{
2018 CountedPtr< Scantable > out = getScantable(in, false);
2019 Table& tout = out->table();
2020 TableIterator iter(tout, "FREQ_ID");
2021 FFTServer<Float,Complex> ffts;
2022 while ( !iter.pastEnd() ) {
2023 Table tab = iter.table();
2024 Double rp,rv,inc;
2025 ROTableRow row(tab);
2026 const TableRecord& rec = row.get(0);
2027 uInt freqid = rec.asuInt("FREQ_ID");
2028 out->frequencies().getEntry(rp, rv, inc, freqid);
2029 ArrayColumn<Float> specCol(tab, "SPECTRA");
2030 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2031 for (int i=0; i<int(tab.nrow()); ++i) {
2032 Vector<Float> spec = specCol(i);
2033 Vector<uChar> flag = flagCol(i);
2034 Int lag0 = Int(spec.nelements()*abs(inc)/(frequency+width)+0.5);
2035 Int lag1 = Int(spec.nelements()*abs(inc)/(frequency-width)+0.5);
2036 for (int k=0; k < flag.nelements(); ++k ) {
2037 if (flag[k] > 0) {
2038 spec[k] = 0.0;
2039 }
2040 }
2041 Vector<Complex> lags;
2042 ffts.fft0(lags, spec);
2043 Int start = max(0, lag0);
2044 Int end = min(Int(lags.nelements()-1), lag1);
2045 if (start == end) {
2046 lags[start] = Complex(0.0);
2047 } else {
2048 for (int j=start; j <=end ;++j) {
2049 lags[j] = Complex(0.0);
2050 }
2051 }
2052 ffts.fft0(spec, lags);
2053 specCol.put(i, spec);
2054 }
2055 ++iter;
2056 }
2057 return out;
2058}
2059
2060// Averaging spectra with different channel/resolution
2061CountedPtr<Scantable>
2062STMath::new_average( const std::vector<CountedPtr<Scantable> >& in,
2063 const bool& compel,
2064 const std::vector<bool>& mask,
2065 const std::string& weight,
2066 const std::string& avmode )
2067 throw ( casa::AipsError )
2068{
2069 if ( avmode == "SCAN" && in.size() != 1 )
2070 throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
2071 "Use merge first."));
2072
2073 CountedPtr<Scantable> out ; // processed result
2074 if ( compel ) {
2075 std::vector< CountedPtr<Scantable> > newin ; // input for average process
2076 uInt insize = in.size() ; // number of scantables
2077
2078 // TEST: do normal average in each table before IF grouping
2079 cout << "Do preliminary averaging" << endl ;
2080 vector< CountedPtr<Scantable> > tmpin( insize ) ;
2081 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2082 vector< CountedPtr<Scantable> > v( 1, in[itable] ) ;
2083 tmpin[itable] = average( v, mask, weight, avmode ) ;
2084 }
2085
2086 // warning
2087 cout << "Average spectra with different spectral resolution" << endl ;
2088 cout << endl ;
2089
2090 // temporarily set coordinfo
2091 vector<string> oldinfo( insize ) ;
2092 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2093 vector<string> coordinfo = in[itable]->getCoordInfo() ;
2094 oldinfo[itable] = coordinfo[0] ;
2095 coordinfo[0] = "Hz" ;
2096 tmpin[itable]->setCoordInfo( coordinfo ) ;
2097 }
2098
2099 // columns
2100 ScalarColumn<uInt> freqIDCol ;
2101 ScalarColumn<uInt> ifnoCol ;
2102 ScalarColumn<uInt> scannoCol ;
2103
2104
2105 // check IF frequency coverage
2106 //cout << "Check IF settings in each table" << endl ;
2107 vector< vector<uInt> > freqid( insize );
2108 vector< vector<double> > iffreq( insize ) ;
2109 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2110 uInt rows = tmpin[itable]->nrow() ;
2111 uInt freqnrows = tmpin[itable]->frequencies().table().nrow() ;
2112 for ( uInt irow = 0 ; irow < rows ; irow++ ) {
2113 if ( freqid[itable].size() == freqnrows ) {
2114 break ;
2115 }
2116 else {
2117 freqIDCol.attach( tmpin[itable]->table(), "FREQ_ID" ) ;
2118 ifnoCol.attach( tmpin[itable]->table(), "IFNO" ) ;
2119 uInt id = freqIDCol( irow ) ;
2120 if ( freqid[itable].size() == 0 || count( freqid[itable].begin(), freqid[itable].end(), id ) == 0 ) {
2121 //cout << "itable = " << itable << ": IF " << id << " is included in the list" << endl ;
2122 vector<double> abcissa = tmpin[itable]->getAbcissa( irow ) ;
2123 freqid[itable].push_back( id ) ;
2124 iffreq[itable].push_back( abcissa[0] - 0.5 * ( abcissa[1] - abcissa[0] ) ) ;
2125 iffreq[itable].push_back( abcissa[abcissa.size()-1] + 0.5 * ( abcissa[1] - abcissa[0] ) ) ;
2126 }
2127 }
2128 }
2129 }
2130
2131 // debug
2132 //cout << "IF settings summary:" << endl ;
2133 //for ( uInt i = 0 ; i < freqid.size() ; i++ ) {
2134 //cout << " Table" << i << endl ;
2135 //for ( uInt j = 0 ; j < freqid[i].size() ; j++ ) {
2136 //cout << " id = " << freqid[i][j] << " (min,max) = (" << iffreq[i][2*j] << "," << iffreq[i][2*j+1] << ")" << endl ;
2137 //}
2138 //}
2139 //cout << endl ;
2140
2141 // IF grouping based on their frequency coverage
2142 //cout << "IF grouping based on their frequency coverage" << endl ;
2143
2144 // IF group
2145 // ifgrp[numgrp][nummember*2]
2146 // ifgrp = [table00, freqrow00, table01, freqrow01, ...],
2147 // [table11, freqrow11, table11, freqrow11, ...],
2148 // ...
2149 // ifgfreq[numgrp*2]
2150 // ifgfreq = [min0, max0, min1, max1, ...]
2151 vector< vector<uInt> > ifgrp ;
2152 vector<double> ifgfreq ;
2153
2154 // parameter for IF grouping
2155 // groupmode = OR retrieve all region
2156 // AND only retrieve overlaped region
2157 //string groupmode = "AND" ;
2158 string groupmode = "OR" ;
2159 uInt sizecr = 0 ;
2160 if ( groupmode == "AND" )
2161 sizecr = 2 ;
2162 else if ( groupmode == "OR" )
2163 sizecr = 0 ;
2164
2165 vector<double> sortedfreq ;
2166 for ( uInt i = 0 ; i < iffreq.size() ; i++ ) {
2167 for ( uInt j = 0 ; j < iffreq[i].size() ; j++ ) {
2168 if ( count( sortedfreq.begin(), sortedfreq.end(), iffreq[i][j] ) == 0 )
2169 sortedfreq.push_back( iffreq[i][j] ) ;
2170 }
2171 }
2172 sort( sortedfreq.begin(), sortedfreq.end() ) ;
2173 for ( vector<double>::iterator i = sortedfreq.begin() ; i != sortedfreq.end()-1 ; i++ ) {
2174 ifgfreq.push_back( *i ) ;
2175 ifgfreq.push_back( *(i+1) ) ;
2176 }
2177 ifgrp.resize( ifgfreq.size()/2 ) ;
2178 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2179 for ( uInt iif = 0 ; iif < freqid[itable].size() ; iif++ ) {
2180 double range0 = iffreq[itable][2*iif] ;
2181 double range1 = iffreq[itable][2*iif+1] ;
2182 for ( uInt j = 0 ; j < ifgrp.size() ; j++ ) {
2183 double fmin = max( range0, ifgfreq[2*j] ) ;
2184 double fmax = min( range1, ifgfreq[2*j+1] ) ;
2185 if ( fmin < fmax ) {
2186 ifgrp[j].push_back( itable ) ;
2187 ifgrp[j].push_back( freqid[itable][iif] ) ;
2188 }
2189 }
2190 }
2191 }
2192 vector< vector<uInt> >::iterator fiter = ifgrp.begin() ;
2193 vector<double>::iterator giter = ifgfreq.begin() ;
2194 while( fiter != ifgrp.end() ) {
2195 if ( fiter->size() <= sizecr ) {
2196 fiter = ifgrp.erase( fiter ) ;
2197 giter = ifgfreq.erase( giter ) ;
2198 giter = ifgfreq.erase( giter ) ;
2199 }
2200 else {
2201 fiter++ ;
2202 advance( giter, 2 ) ;
2203 }
2204 }
2205
2206 // freqgrp[numgrp][nummember]
2207 // freqgrp = [ifgrp00, ifgrp01, ifgrp02, ...],
2208 // [ifgrp10, ifgrp11, ifgrp12, ...],
2209 // ...
2210 // freqrange[numgrp*2]
2211 // freqrange = [min0, max0, min1, max1, ...]
2212 vector< vector<uInt> > freqgrp ;
2213 double freqrange = 0.0 ;
2214 uInt grpnum = 0 ;
2215 for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
2216 // Assumed that ifgfreq was sorted
2217 if ( grpnum != 0 && freqrange == ifgfreq[2*i] ) {
2218 freqgrp[grpnum-1].push_back( i ) ;
2219 }
2220 else {
2221 vector<uInt> grp0( 1, i ) ;
2222 freqgrp.push_back( grp0 ) ;
2223 grpnum++ ;
2224 }
2225 freqrange = ifgfreq[2*i+1] ;
2226 }
2227
2228
2229 // print IF groups
2230 cout << "IF Group summary: " << endl ;
2231 cout << " GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
2232 for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
2233 cout << " GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
2234 for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
2235 cout << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ;
2236 }
2237 cout << endl ;
2238 }
2239 cout << endl ;
2240
2241 // print frequency group
2242 cout << "Frequency Group summary: " << endl ;
2243 cout << " GROUP_ID [FREQ_MIN, FREQ_MAX]: IF_GROUP_ID" << endl ;
2244 for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
2245 cout << " GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*freqgrp[i][0]] << "," << ifgfreq[2*freqgrp[i][freqgrp[i].size()-1]+1] << "]: " ;
2246 for ( uInt j = 0 ; j < freqgrp[i].size() ; j++ ) {
2247 cout << freqgrp[i][j] << " " ;
2248 }
2249 cout << endl ;
2250 }
2251 cout << endl ;
2252
2253 // membership check
2254 // groups[numtable][numIF][nummembership]
2255 vector< vector< vector<uInt> > > groups( insize ) ;
2256 for ( uInt i = 0 ; i < insize ; i++ ) {
2257 groups[i].resize( freqid[i].size() ) ;
2258 }
2259 for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
2260 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
2261 uInt tableid = ifgrp[igrp][2*imem] ;
2262 vector<uInt>::iterator iter = find( freqid[tableid].begin(), freqid[tableid].end(), ifgrp[igrp][2*imem+1] ) ;
2263 if ( iter != freqid[tableid].end() ) {
2264 uInt rowid = distance( freqid[tableid].begin(), iter ) ;
2265 groups[tableid][rowid].push_back( igrp ) ;
2266 }
2267 }
2268 }
2269
2270 // print membership
2271 //for ( uInt i = 0 ; i < insize ; i++ ) {
2272 //cout << "Table " << i << endl ;
2273 //for ( uInt j = 0 ; j < groups[i].size() ; j++ ) {
2274 //cout << " FREQ_ID " << setw( 2 ) << freqid[i][j] << ": " ;
2275 //for ( uInt k = 0 ; k < groups[i][j].size() ; k++ ) {
2276 //cout << setw( 2 ) << groups[i][j][k] << " " ;
2277 //}
2278 //cout << endl ;
2279 //}
2280 //}
2281
2282 // set back coordinfo
2283 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2284 vector<string> coordinfo = tmpin[itable]->getCoordInfo() ;
2285 coordinfo[0] = oldinfo[itable] ;
2286 tmpin[itable]->setCoordInfo( coordinfo ) ;
2287 }
2288
2289 // Create additional table if needed
2290 bool oldInsitu = insitu_ ;
2291 setInsitu( false ) ;
2292 vector< vector<uInt> > addrow( insize ) ;
2293 vector<uInt> addtable( insize, 0 ) ;
2294 vector<uInt> newtableids( insize ) ;
2295 vector<uInt> newifids( insize, 0 ) ;
2296 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2297 //cout << "Table " << setw(2) << itable << ": " ;
2298 for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
2299 addrow[itable].push_back( groups[itable][ifrow].size()-1 ) ;
2300 //cout << addrow[itable][ifrow] << " " ;
2301 }
2302 addtable[itable] = *max_element( addrow[itable].begin(), addrow[itable].end() ) ;
2303 //cout << "(" << addtable[itable] << ")" << endl ;
2304 }
2305 newin.resize( insize ) ;
2306 copy( tmpin.begin(), tmpin.end(), newin.begin() ) ;
2307 for ( uInt i = 0 ; i < insize ; i++ ) {
2308 newtableids[i] = i ;
2309 }
2310 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2311 for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
2312 CountedPtr<Scantable> add = getScantable( newin[itable], false ) ;
2313 vector<int> freqidlist ;
2314 for ( uInt i = 0 ; i < groups[itable].size() ; i++ ) {
2315 if ( groups[itable][i].size() > iadd + 1 ) {
2316 freqidlist.push_back( freqid[itable][i] ) ;
2317 }
2318 }
2319 stringstream taqlstream ;
2320 taqlstream << "SELECT FROM $1 WHERE FREQ_ID IN [" ;
2321 for ( uInt i = 0 ; i < freqidlist.size() ; i++ ) {
2322 taqlstream << i ;
2323 if ( i < freqidlist.size() - 1 )
2324 taqlstream << "," ;
2325 else
2326 taqlstream << "]" ;
2327 }
2328 string taql = taqlstream.str() ;
2329 //cout << "taql = " << taql << endl ;
2330 STSelector selector = STSelector() ;
2331 selector.setTaQL( taql ) ;
2332 add->setSelection( selector ) ;
2333 newin.push_back( add ) ;
2334 newtableids.push_back( itable ) ;
2335 newifids.push_back( iadd + 1 ) ;
2336 }
2337 }
2338
2339 // udpate ifgrp
2340 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2341 for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
2342 for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
2343 if ( groups[itable][ifrow].size() > iadd + 1 ) {
2344 uInt igrp = groups[itable][ifrow][iadd+1] ;
2345 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
2346 if ( ifgrp[igrp][2*imem] == newtableids[iadd+insize] && ifgrp[igrp][2*imem+1] == freqid[newtableids[iadd+insize]][ifrow] ) {
2347 ifgrp[igrp][2*imem] = insize + iadd ;
2348 }
2349 }
2350 }
2351 }
2352 }
2353 }
2354
2355 // print IF groups again for debug
2356 //cout << "IF Group summary: " << endl ;
2357 //cout << " GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
2358 //for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
2359 //cout << " GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
2360 //for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
2361 //cout << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ;
2362 //}
2363 //cout << endl ;
2364 //}
2365 //cout << endl ;
2366
2367 // reset SCANNO and IF: IF number is reset by the result of sortation
2368 cout << "All scan number is set to 0" << endl ;
2369 //cout << "All IF number is set to IF group index" << endl ;
2370 cout << endl ;
2371 insize = newin.size() ;
2372 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2373 uInt rows = newin[itable]->nrow() ;
2374 Table &tmpt = newin[itable]->table() ;
2375 freqIDCol.attach( tmpt, "FREQ_ID" ) ;
2376 scannoCol.attach( tmpt, "SCANNO" ) ;
2377 ifnoCol.attach( tmpt, "IFNO" ) ;
2378 for ( uInt irow=0 ; irow < rows ; irow++ ) {
2379 scannoCol.put( irow, 0 ) ;
2380 uInt freqID = freqIDCol( irow ) ;
2381 vector<uInt>::iterator iter = find( freqid[newtableids[itable]].begin(), freqid[newtableids[itable]].end(), freqID ) ;
2382 if ( iter != freqid[newtableids[itable]].end() ) {
2383 uInt index = distance( freqid[newtableids[itable]].begin(), iter ) ;
2384 ifnoCol.put( irow, groups[newtableids[itable]][index][newifids[itable]] ) ;
2385 }
2386 else {
2387 throw(AipsError("IF grouping was wrong in additional tables.")) ;
2388 }
2389 }
2390 }
2391 oldinfo.resize( insize ) ;
2392 setInsitu( oldInsitu ) ;
2393
2394 // temporarily set coordinfo
2395 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2396 vector<string> coordinfo = newin[itable]->getCoordInfo() ;
2397 oldinfo[itable] = coordinfo[0] ;
2398 coordinfo[0] = "Hz" ;
2399 newin[itable]->setCoordInfo( coordinfo ) ;
2400 }
2401
2402 // save column values in the vector
2403 vector< vector<uInt> > freqTableIdVec( insize ) ;
2404 vector< vector<uInt> > freqIdVec( insize ) ;
2405 vector< vector<uInt> > ifNoVec( insize ) ;
2406 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2407 ScalarColumn<uInt> freqIDs ;
2408 freqIDs.attach( newin[itable]->frequencies().table(), "ID" ) ;
2409 ifnoCol.attach( newin[itable]->table(), "IFNO" ) ;
2410 freqIDCol.attach( newin[itable]->table(), "FREQ_ID" ) ;
2411 for ( uInt irow = 0 ; irow < newin[itable]->frequencies().table().nrow() ; irow++ ) {
2412 freqTableIdVec[itable].push_back( freqIDs( irow ) ) ;
2413 }
2414 for ( uInt irow = 0 ; irow < newin[itable]->table().nrow() ; irow++ ) {
2415 freqIdVec[itable].push_back( freqIDCol( irow ) ) ;
2416 ifNoVec[itable].push_back( ifnoCol( irow ) ) ;
2417 }
2418 }
2419
2420 // reset spectra and flagtra: pick up common part of frequency coverage
2421 //cout << "Pick common frequency range and align resolution" << endl ;
2422 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2423 uInt rows = newin[itable]->nrow() ;
2424 int nminchan = -1 ;
2425 int nmaxchan = -1 ;
2426 vector<uInt> freqIdUpdate ;
2427 for ( uInt irow = 0 ; irow < rows ; irow++ ) {
2428 uInt ifno = ifNoVec[itable][irow] ; // IFNO is reset by group index
2429 double minfreq = ifgfreq[2*ifno] ;
2430 double maxfreq = ifgfreq[2*ifno+1] ;
2431 //cout << "frequency range: [" << minfreq << "," << maxfreq << "]" << endl ;
2432 vector<double> abcissa = newin[itable]->getAbcissa( irow ) ;
2433 int nchan = abcissa.size() ;
2434 double resol = abcissa[1] - abcissa[0] ;
2435 //cout << "abcissa range : [" << abcissa[0] << "," << abcissa[nchan-1] << "]" << endl ;
2436 if ( minfreq <= abcissa[0] )
2437 nminchan = 0 ;
2438 else {
2439 //double cfreq = ( minfreq - abcissa[0] ) / resol ;
2440 double cfreq = ( minfreq - abcissa[0] + 0.5 * resol ) / resol ;
2441 nminchan = int(cfreq) + ( ( cfreq - int(cfreq) <= 0.5 ) ? 0 : 1 ) ;
2442 }
2443 if ( maxfreq >= abcissa[abcissa.size()-1] )
2444 nmaxchan = abcissa.size() - 1 ;
2445 else {
2446 //double cfreq = ( abcissa[abcissa.size()-1] - maxfreq ) / resol ;
2447 double cfreq = ( abcissa[abcissa.size()-1] - maxfreq + 0.5 * resol ) / resol ;
2448 nmaxchan = abcissa.size() - 1 - int(cfreq) - ( ( cfreq - int(cfreq) >= 0.5 ) ? 1 : 0 ) ;
2449 }
2450 //cout << "channel range (" << irow << "): [" << nminchan << "," << nmaxchan << "]" << endl ;
2451 if ( nmaxchan > nminchan ) {
2452 newin[itable]->reshapeSpectrum( nminchan, nmaxchan, irow ) ;
2453 int newchan = nmaxchan - nminchan + 1 ;
2454 if ( count( freqIdUpdate.begin(), freqIdUpdate.end(), freqIdVec[itable][irow] ) == 0 && newchan < nchan )
2455 freqIdUpdate.push_back( freqIdVec[itable][irow] ) ;
2456 }
2457 else {
2458 throw(AipsError("Failed to pick up common part of frequency range.")) ;
2459 }
2460 }
2461 for ( uInt i = 0 ; i < freqIdUpdate.size() ; i++ ) {
2462 uInt freqId = freqIdUpdate[i] ;
2463 Double refpix ;
2464 Double refval ;
2465 Double increment ;
2466
2467 // update row
2468 newin[itable]->frequencies().getEntry( refpix, refval, increment, freqId ) ;
2469 refval = refval - ( refpix - nminchan ) * increment ;
2470 refpix = 0 ;
2471 newin[itable]->frequencies().setEntry( refpix, refval, increment, freqId ) ;
2472 }
2473 }
2474
2475
2476 // reset spectra and flagtra: align spectral resolution
2477 //cout << "Align spectral resolution" << endl ;
2478 vector<double> gmaxdnu( freqgrp.size(), 0.0 ) ;
2479 vector<uInt> gmemid( freqgrp.size(), 0 ) ;
2480 for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
2481 double maxdnu = 0.0 ; // maximum (coarsest) frequency resolution
2482 int minchan = INT_MAX ; // minimum channel number
2483 Double refpixref = -1 ; // reference of 'reference pixel'
2484 Double refvalref = -1 ; // reference of 'reference frequency'
2485 Double refinc = -1 ; // reference frequency resolution
2486 uInt refreqid ;
2487 uInt reftable = INT_MAX;
2488 // process only if group member > 1
2489 if ( ifgrp[igrp].size() > 2 ) {
2490 // find minchan and maxdnu in each group
2491 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
2492 uInt tableid = ifgrp[igrp][2*imem] ;
2493 uInt rowid = ifgrp[igrp][2*imem+1] ;
2494 vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
2495 if ( iter != freqIdVec[tableid].end() ) {
2496 uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
2497 vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
2498 int nchan = abcissa.size() ;
2499 double dnu = abcissa[1] - abcissa[0] ;
2500 //cout << "GROUP " << igrp << " (" << tableid << "," << rowid << "): nchan = " << nchan << " (minchan = " << minchan << ")" << endl ;
2501 if ( nchan < minchan ) {
2502 minchan = nchan ;
2503 maxdnu = dnu ;
2504 newin[tableid]->frequencies().getEntry( refpixref, refvalref, refinc, rowid ) ;
2505 refreqid = rowid ;
2506 reftable = tableid ;
2507 }
2508 }
2509 }
2510 // regrid spectra in each group
2511 cout << "GROUP " << igrp << endl ;
2512 cout << " Channel number is adjusted to " << minchan << endl ;
2513 cout << " Corresponding frequency resolution is " << maxdnu << "Hz" << endl ;
2514 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
2515 uInt tableid = ifgrp[igrp][2*imem] ;
2516 uInt rowid = ifgrp[igrp][2*imem+1] ;
2517 freqIDCol.attach( newin[tableid]->table(), "FREQ_ID" ) ;
2518 //cout << "tableid = " << tableid << " rowid = " << rowid << ": " << endl ;
2519 //cout << " regridChannel applied to " ;
2520 if ( tableid != reftable )
2521 refreqid = newin[tableid]->frequencies().addEntry( refpixref, refvalref, refinc ) ;
2522 for ( uInt irow = 0 ; irow < newin[tableid]->table().nrow() ; irow++ ) {
2523 uInt tfreqid = freqIdVec[tableid][irow] ;
2524 if ( tfreqid == rowid ) {
2525 //cout << irow << " " ;
2526 newin[tableid]->regridChannel( minchan, maxdnu, irow ) ;
2527 freqIDCol.put( irow, refreqid ) ;
2528 freqIdVec[tableid][irow] = refreqid ;
2529 }
2530 }
2531 //cout << endl ;
2532 }
2533 }
2534 else {
2535 uInt tableid = ifgrp[igrp][0] ;
2536 uInt rowid = ifgrp[igrp][1] ;
2537 vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
2538 if ( iter != freqIdVec[tableid].end() ) {
2539 uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
2540 vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
2541 minchan = abcissa.size() ;
2542 maxdnu = abcissa[1] - abcissa[0] ;
2543 }
2544 }
2545 for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
2546 if ( count( freqgrp[i].begin(), freqgrp[i].end(), igrp ) > 0 ) {
2547 if ( maxdnu > gmaxdnu[i] ) {
2548 gmaxdnu[i] = maxdnu ;
2549 gmemid[i] = igrp ;
2550 }
2551 break ;
2552 }
2553 }
2554 }
2555 cout << endl ;
2556
2557 // set back coordinfo
2558 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2559 vector<string> coordinfo = newin[itable]->getCoordInfo() ;
2560 coordinfo[0] = oldinfo[itable] ;
2561 newin[itable]->setCoordInfo( coordinfo ) ;
2562 }
2563
2564 // accumulate all rows into the first table
2565 // NOTE: assumed in.size() = 1
2566 vector< CountedPtr<Scantable> > tmp( 1 ) ;
2567 if ( newin.size() == 1 )
2568 tmp[0] = newin[0] ;
2569 else
2570 tmp[0] = merge( newin ) ;
2571
2572 //return tmp[0] ;
2573
2574 // average
2575 CountedPtr<Scantable> tmpout = average( tmp, mask, weight, avmode ) ;
2576
2577 //return tmpout ;
2578
2579 // combine frequency group
2580 cout << "Combine spectra based on frequency grouping" << endl ;
2581 cout << "IFNO is renumbered as frequency group ID (see above)" << endl ;
2582 vector<string> coordinfo = tmpout->getCoordInfo() ;
2583 oldinfo[0] = coordinfo[0] ;
2584 coordinfo[0] = "Hz" ;
2585 tmpout->setCoordInfo( coordinfo ) ;
2586 // create proformas of output table
2587 stringstream taqlstream ;
2588 taqlstream << "SELECT FROM $1 WHERE IFNO IN [" ;
2589 for ( uInt i = 0 ; i < gmemid.size() ; i++ ) {
2590 taqlstream << gmemid[i] ;
2591 if ( i < gmemid.size() - 1 )
2592 taqlstream << "," ;
2593 else
2594 taqlstream << "]" ;
2595 }
2596 string taql = taqlstream.str() ;
2597 //cout << "taql = " << taql << endl ;
2598 STSelector selector = STSelector() ;
2599 selector.setTaQL( taql ) ;
2600 oldInsitu = insitu_ ;
2601 setInsitu( false ) ;
2602 out = getScantable( tmpout, false ) ;
2603 setInsitu( oldInsitu ) ;
2604 out->setSelection( selector ) ;
2605 // regrid rows
2606 ifnoCol.attach( tmpout->table(), "IFNO" ) ;
2607 for ( uInt irow = 0 ; irow < tmpout->table().nrow() ; irow++ ) {
2608 uInt ifno = ifnoCol( irow ) ;
2609 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
2610 if ( count( freqgrp[igrp].begin(), freqgrp[igrp].end(), ifno ) > 0 ) {
2611 vector<double> abcissa = tmpout->getAbcissa( irow ) ;
2612 double bw = ( abcissa[1] - abcissa[0] ) * abcissa.size() ;
2613 int nchan = (int)( bw / gmaxdnu[igrp] ) ;
2614 tmpout->regridChannel( nchan, gmaxdnu[igrp], irow ) ;
2615 break ;
2616 }
2617 }
2618 }
2619 // combine spectra
2620 ArrayColumn<Float> specColOut ;
2621 specColOut.attach( out->table(), "SPECTRA" ) ;
2622 ArrayColumn<uChar> flagColOut ;
2623 flagColOut.attach( out->table(), "FLAGTRA" ) ;
2624 ScalarColumn<uInt> ifnoColOut ;
2625 ifnoColOut.attach( out->table(), "IFNO" ) ;
2626 ScalarColumn<uInt> polnoColOut ;
2627 polnoColOut.attach( out->table(), "POLNO" ) ;
2628 ScalarColumn<uInt> freqidColOut ;
2629 freqidColOut.attach( out->table(), "FREQ_ID" ) ;
2630 Table &tab = tmpout->table() ;
2631 TableIterator iter( tab, "POLNO" ) ;
2632 vector< vector<uInt> > sizes( freqgrp.size() ) ;
2633 while( !iter.pastEnd() ) {
2634 vector< vector<Float> > specout( freqgrp.size() ) ;
2635 vector< vector<uChar> > flagout( freqgrp.size() ) ;
2636 ArrayColumn<Float> specCols ;
2637 specCols.attach( iter.table(), "SPECTRA" ) ;
2638 ArrayColumn<uChar> flagCols ;
2639 flagCols.attach( iter.table(), "FLAGTRA" ) ;
2640 ifnoCol.attach( iter.table(), "IFNO" ) ;
2641 ScalarColumn<uInt> polnos ;
2642 polnos.attach( iter.table(), "POLNO" ) ;
2643 uInt polno = polnos( 0 ) ;
2644 //cout << "POLNO iteration: " << polno << endl ;
2645 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
2646 sizes[igrp].resize( freqgrp[igrp].size() ) ;
2647 for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
2648 for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
2649 uInt ifno = ifnoCol( irow ) ;
2650 if ( ifno == freqgrp[igrp][imem] ) {
2651 Vector<Float> spec = specCols( irow ) ;
2652 Vector<uChar> flag = flagCols( irow ) ;
2653 vector<Float> svec ;
2654 spec.tovector( svec ) ;
2655 vector<uChar> fvec ;
2656 flag.tovector( fvec ) ;
2657 //cout << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << endl ;
2658 specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ;
2659 flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ;
2660 //cout << "specout[" << igrp << "].size() = " << specout[igrp].size() << endl ;
2661 sizes[igrp][imem] = spec.nelements() ;
2662 }
2663 }
2664 }
2665 for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
2666 uInt ifout = ifnoColOut( irow ) ;
2667 uInt polout = polnoColOut( irow ) ;
2668 if ( ifout == gmemid[igrp] && polout == polno ) {
2669 // set SPECTRA and FRAGTRA
2670 Vector<Float> newspec( specout[igrp] ) ;
2671 Vector<uChar> newflag( flagout[igrp] ) ;
2672 specColOut.put( irow, newspec ) ;
2673 flagColOut.put( irow, newflag ) ;
2674 // IFNO renumbering
2675 ifnoColOut.put( irow, igrp ) ;
2676 }
2677 }
2678 }
2679 iter++ ;
2680 }
2681 // update FREQUENCIES subtable
2682 vector<bool> updated( freqgrp.size(), false ) ;
2683 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
2684 uInt index = 0 ;
2685 uInt pixShift = 0 ;
2686 while ( freqgrp[igrp][index] != gmemid[igrp] ) {
2687 pixShift += sizes[igrp][index++] ;
2688 }
2689 for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
2690 if ( ifnoColOut( irow ) == gmemid[igrp] && !updated[igrp] ) {
2691 uInt freqidOut = freqidColOut( irow ) ;
2692 //cout << "freqgrp " << igrp << " freqidOut = " << freqidOut << endl ;
2693 double refpix ;
2694 double refval ;
2695 double increm ;
2696 out->frequencies().getEntry( refpix, refval, increm, freqidOut ) ;
2697 refpix += pixShift ;
2698 out->frequencies().setEntry( refpix, refval, increm, freqidOut ) ;
2699 updated[igrp] = true ;
2700 }
2701 }
2702 }
2703
2704 //out = tmpout ;
2705
2706 coordinfo = tmpout->getCoordInfo() ;
2707 coordinfo[0] = oldinfo[0] ;
2708 tmpout->setCoordInfo( coordinfo ) ;
2709 }
2710 else {
2711 // simple average
2712 out = average( in, mask, weight, avmode ) ;
2713 }
2714
2715 return out ;
2716}
Note: See TracBrowser for help on using the repository browser.