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

Last change on this file since 1607 was 1607, checked in by Takeshi Nakazato, 16 years ago

New Development: No

JIRA Issue: Yes CAS-1423

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: run sdaverage with OTF data

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

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