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

Last change on this file since 1624 was 1616, checked in by Takeshi Nakazato, 15 years ago

New Development: No

JIRA Issue: Yes CAS-729, CAS-1147

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

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