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

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

New Development: No

JIRA Issue: Yes CAS-1448

Ready to Release: No

Interface Changes: Yes/No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

I have defined a function for folding frequency-switch data.
It is still experimental. Test data is needed.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 108.3 KB
Line 
1//
2// C++ Implementation: STMath
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12
13#include <casa/iomanip.h>
14#include <casa/Exceptions/Error.h>
15#include <casa/Containers/Block.h>
16#include <casa/BasicSL/String.h>
17#include <casa/Arrays/MaskArrLogi.h>
18#include <casa/Arrays/MaskArrMath.h>
19#include <casa/Arrays/ArrayLogical.h>
20#include <casa/Arrays/ArrayMath.h>
21#include <casa/Arrays/Slice.h>
22#include <casa/Arrays/Slicer.h>
23#include <casa/Containers/RecordField.h>
24#include <tables/Tables/TableRow.h>
25#include <tables/Tables/TableVector.h>
26#include <tables/Tables/TabVecMath.h>
27#include <tables/Tables/ExprNode.h>
28#include <tables/Tables/TableRecord.h>
29#include <tables/Tables/TableParse.h>
30#include <tables/Tables/ReadAsciiTable.h>
31#include <tables/Tables/TableIter.h>
32#include <tables/Tables/TableCopy.h>
33#include <scimath/Mathematics/FFTServer.h>
34
35#include <lattices/Lattices/LatticeUtilities.h>
36
37#include <coordinates/Coordinates/SpectralCoordinate.h>
38#include <coordinates/Coordinates/CoordinateSystem.h>
39#include <coordinates/Coordinates/CoordinateUtil.h>
40#include <coordinates/Coordinates/FrequencyAligner.h>
41
42#include <scimath/Mathematics/VectorKernel.h>
43#include <scimath/Mathematics/Convolver.h>
44#include <scimath/Functionals/Polynomial.h>
45
46#include "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 ScalarColumn<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 //cerr << out1->frequencies().table().nrow() << " " << out2->frequencies().table().nrow() << endl ;
1204 if (rp1==rp2) {
1205 Double foffset = rv1 - rv2;
1206 uInt choffset = static_cast<uInt>(foffset/abs(inc2));
1207 if (choffset >= nchan) {
1208 cerr<<"out-band frequency switching, no folding"<<endl;
1209 nofold = True;
1210 }
1211 }
1212
1213 if (nofold) {
1214 std::vector< CountedPtr< Scantable > > tabs;
1215 tabs.push_back(out1);
1216 tabs.push_back(out2);
1217 out = merge(tabs);
1218 }
1219 else { //folding is not implemented yet
1220 //out = out1;
1221 Int choffset = static_cast<Int>((rv1-rv2)/inc2) ;
1222 out = dofold( out1, out2, choffset ) ;
1223 }
1224
1225 return out;
1226}
1227
1228CountedPtr<Scantable> STMath::dofold( const CountedPtr<Scantable> &sig,
1229 const CountedPtr<Scantable> &ref,
1230 Int choffset )
1231{
1232 // output scantable
1233 CountedPtr<Scantable> out = getScantable( sig, false ) ;
1234
1235 // get column
1236 ROArrayColumn<Float> specCol1( sig->table(), "SPECTRA" ) ;
1237 ROArrayColumn<Float> specCol2( ref->table(), "SPECTRA" ) ;
1238 ROArrayColumn<Float> tsysCol1( sig->table(), "TSYS" ) ;
1239 ROArrayColumn<Float> tsysCol2( ref->table(), "TSYS" ) ;
1240 ROArrayColumn<uChar> flagCol1( sig->table(), "FLAGTRA" ) ;
1241 ROArrayColumn<uChar> flagCol2( ref->table(), "FLAGTRA" ) ;
1242 ROScalarColumn<Double> mjdCol1( sig->table(), "TIME" ) ;
1243 ROScalarColumn<Double> mjdCol2( ref->table(), "TIME" ) ;
1244 ROScalarColumn<Double> intervalCol1( sig->table(), "INTERVAL" ) ;
1245 ROScalarColumn<Double> intervalCol2( ref->table(), "INTERVAL" ) ;
1246
1247 // check
1248 if ( choffset == 0 ) {
1249 cerr << "channel offset is zero, no folding" << endl ;
1250 return out ;
1251 }
1252 int nchan = ref->nchan() ;
1253 if ( abs(choffset) >= nchan ) {
1254 cerr << "out-band frequency switching, no folding" << endl ;
1255 return out ;
1256 }
1257
1258 // attach column for output scantable
1259 ArrayColumn<Float> specColOut( out->table(), "SPECTRA" ) ;
1260 ArrayColumn<uChar> flagColOut( out->table(), "FLAGTRA" ) ;
1261 ArrayColumn<Float> tsysColOut( out->table(), "TSYS" ) ;
1262 ScalarColumn<Double> mjdColOut( out->table(), "TIME" ) ;
1263 ScalarColumn<Double> intervalColOut( out->table(), "INTERVAL" ) ;
1264
1265 // for each row
1266 // assume that the data order are same between sig and ref
1267 RowAccumulator acc( asap::TINTSYS ) ;
1268 for ( int i = 0 ; i < sig->nrow() ; i++ ) {
1269 // get values
1270 Vector<Float> spsig ;
1271 specCol1.get( i, spsig ) ;
1272 Vector<Float> spref ;
1273 specCol2.get( i, spref ) ;
1274 Vector<Float> tsyssig ;
1275 tsysCol1.get( i, tsyssig ) ;
1276 Vector<Float> tsysref ;
1277 tsysCol2.get( i, tsysref ) ;
1278 Vector<uChar> flagsig ;
1279 flagCol1.get( i, flagsig ) ;
1280 Vector<uChar> flagref ;
1281 flagCol2.get( i, flagref ) ;
1282 Double timesig ;
1283 mjdCol1.get( i, timesig ) ;
1284 Double timeref ;
1285 mjdCol2.get( i, timeref ) ;
1286 Double intsig ;
1287 intervalCol1.get( i, intsig ) ;
1288 Double intref ;
1289 intervalCol2.get( i, intref ) ;
1290
1291 // shift reference spectra
1292 int refchan = spref.nelements() ;
1293 if ( choffset > 0 ) {
1294 for ( int j = 0 ; j < refchan-choffset ; j++ ) {
1295 spref[j] = spref[j+choffset] ;
1296 tsysref[j] = tsysref[j+choffset] ;
1297 flagref[j] = flagref[j+choffset] ;
1298 }
1299 for ( int j = refchan-choffset ; j < refchan ; j++ ) {
1300 spref[j] = spref[j-refchan+choffset] ;
1301 tsysref[j] = tsysref[j-refchan+choffset] ;
1302 flagref[j] = flagref[j-refchan+choffset] ;
1303 }
1304 }
1305 else {
1306 for ( int j = 0 ; j < abs(choffset) ; j++ ) {
1307 spref[j] = spref[refchan+choffset+j] ;
1308 tsysref[j] = tsysref[refchan+choffset+j] ;
1309 flagref[j] = flagref[refchan+choffset+j] ;
1310 }
1311 for ( int j = abs(choffset) ; j < refchan ; j++ ) {
1312 spref[j] = spref[j+choffset] ;
1313 tsysref[j] = tsysref[j+choffset] ;
1314 flagref[j] = flagref[j+choffset] ;
1315 }
1316 }
1317
1318 // folding
1319 acc.add( spsig, !flagsig, tsyssig, intsig, timesig ) ;
1320 acc.add( spref, !flagref, tsysref, intref, timeref ) ;
1321
1322 // put result
1323 specColOut.put( i, acc.getSpectrum() ) ;
1324 const Vector<Bool> &msk = acc.getMask() ;
1325 Vector<uChar> flg( msk.shape() ) ;
1326 convertArray( flg, !msk ) ;
1327 flagColOut.put( i, flg ) ;
1328 tsysColOut.put( i, acc.getTsys() ) ;
1329 intervalColOut.put( i, acc.getInterval() ) ;
1330 mjdColOut.put( i, acc.getTime() ) ;
1331
1332 acc.reset() ;
1333 }
1334
1335 return out ;
1336}
1337
1338
1339CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
1340{
1341 // make copy or reference
1342 CountedPtr< Scantable > out = getScantable(in, false);
1343 Table& tout = out->table();
1344 Block<String> cols(4);
1345 cols[0] = String("SCANNO");
1346 cols[1] = String("CYCLENO");
1347 cols[2] = String("BEAMNO");
1348 cols[3] = String("POLNO");
1349 TableIterator iter(tout, cols);
1350 while (!iter.pastEnd()) {
1351 Table subt = iter.table();
1352 // this should leave us with two rows for the two IFs....if not ignore
1353 if (subt.nrow() != 2 ) {
1354 continue;
1355 }
1356 ArrayColumn<Float> specCol(subt, "SPECTRA");
1357 ArrayColumn<Float> tsysCol(subt, "TSYS");
1358 ArrayColumn<uChar> flagCol(subt, "FLAGTRA");
1359 Vector<Float> onspec,offspec, ontsys, offtsys;
1360 Vector<uChar> onflag, offflag;
1361 tsysCol.get(0, ontsys); tsysCol.get(1, offtsys);
1362 specCol.get(0, onspec); specCol.get(1, offspec);
1363 flagCol.get(0, onflag); flagCol.get(1, offflag);
1364 MaskedArray<Float> on = maskedArray(onspec, onflag);
1365 MaskedArray<Float> off = maskedArray(offspec, offflag);
1366 MaskedArray<Float> oncopy = on.copy();
1367
1368 on /= off; on -= 1.0f;
1369 on *= ontsys[0];
1370 off /= oncopy; off -= 1.0f;
1371 off *= offtsys[0];
1372 specCol.put(0, on.getArray());
1373 const Vector<Bool>& m0 = on.getMask();
1374 Vector<uChar> flags0(m0.shape());
1375 convertArray(flags0, !m0);
1376 flagCol.put(0, flags0);
1377
1378 specCol.put(1, off.getArray());
1379 const Vector<Bool>& m1 = off.getMask();
1380 Vector<uChar> flags1(m1.shape());
1381 convertArray(flags1, !m1);
1382 flagCol.put(1, flags1);
1383 ++iter;
1384 }
1385
1386 return out;
1387}
1388
1389std::vector< float > STMath::statistic( const CountedPtr< Scantable > & in,
1390 const std::vector< bool > & mask,
1391 const std::string& which )
1392{
1393
1394 Vector<Bool> m(mask);
1395 const Table& tab = in->table();
1396 ROArrayColumn<Float> specCol(tab, "SPECTRA");
1397 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1398 std::vector<float> out;
1399 for (uInt i=0; i < tab.nrow(); ++i ) {
1400 Vector<Float> spec; specCol.get(i, spec);
1401 Vector<uChar> flag; flagCol.get(i, flag);
1402 MaskedArray<Float> ma = maskedArray(spec, flag);
1403 float outstat = 0.0;
1404 if ( spec.nelements() == m.nelements() ) {
1405 outstat = mathutil::statistics(which, ma(m));
1406 } else {
1407 outstat = mathutil::statistics(which, ma);
1408 }
1409 out.push_back(outstat);
1410 }
1411 return out;
1412}
1413
1414std::vector< int > STMath::minMaxChan( const CountedPtr< Scantable > & in,
1415 const std::vector< bool > & mask,
1416 const std::string& which )
1417{
1418
1419 Vector<Bool> m(mask);
1420 const Table& tab = in->table();
1421 ROArrayColumn<Float> specCol(tab, "SPECTRA");
1422 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1423 std::vector<int> out;
1424 for (uInt i=0; i < tab.nrow(); ++i ) {
1425 Vector<Float> spec; specCol.get(i, spec);
1426 Vector<uChar> flag; flagCol.get(i, flag);
1427 MaskedArray<Float> ma = maskedArray(spec, flag);
1428 if (ma.ndim() != 1) {
1429 throw (ArrayError(
1430 "std::vector<int> STMath::minMaxChan("
1431 "ContedPtr<Scantable> &in, std::vector<bool> &mask, "
1432 " std::string &which)"
1433 " - MaskedArray is not 1D"));
1434 }
1435 IPosition outpos(1,0);
1436 if ( spec.nelements() == m.nelements() ) {
1437 outpos = mathutil::minMaxPos(which, ma(m));
1438 } else {
1439 outpos = mathutil::minMaxPos(which, ma);
1440 }
1441 out.push_back(outpos[0]);
1442 }
1443 return out;
1444}
1445
1446CountedPtr< Scantable > STMath::bin( const CountedPtr< Scantable > & in,
1447 int width )
1448{
1449 if ( !in->getSelection().empty() ) throw(AipsError("Can't bin subset of the data."));
1450 CountedPtr< Scantable > out = getScantable(in, false);
1451 Table& tout = out->table();
1452 out->frequencies().rescale(width, "BIN");
1453 ArrayColumn<Float> specCol(tout, "SPECTRA");
1454 ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
1455 for (uInt i=0; i < tout.nrow(); ++i ) {
1456 MaskedArray<Float> main = maskedArray(specCol(i), flagCol(i));
1457 MaskedArray<Float> maout;
1458 LatticeUtilities::bin(maout, main, 0, Int(width));
1459 /// @todo implement channel based tsys binning
1460 specCol.put(i, maout.getArray());
1461 flagCol.put(i, flagsFromMA(maout));
1462 // take only the first binned spectrum's length for the deprecated
1463 // global header item nChan
1464 if (i==0) tout.rwKeywordSet().define(String("nChan"),
1465 Int(maout.getArray().nelements()));
1466 }
1467 return out;
1468}
1469
1470CountedPtr< Scantable > STMath::resample( const CountedPtr< Scantable >& in,
1471 const std::string& method,
1472 float width )
1473//
1474// Should add the possibility of width being specified in km/s. This means
1475// that for each freqID (SpectralCoordinate) we will need to convert to an
1476// average channel width (say at the reference pixel). Then we would need
1477// to be careful to make sure each spectrum (of different freqID)
1478// is the same length.
1479//
1480{
1481 //InterpolateArray1D<Double,Float>::InterpolationMethod interp;
1482 Int interpMethod(stringToIMethod(method));
1483
1484 CountedPtr< Scantable > out = getScantable(in, false);
1485 Table& tout = out->table();
1486
1487// Resample SpectralCoordinates (one per freqID)
1488 out->frequencies().rescale(width, "RESAMPLE");
1489 TableIterator iter(tout, "IFNO");
1490 TableRow row(tout);
1491 while ( !iter.pastEnd() ) {
1492 Table tab = iter.table();
1493 ArrayColumn<Float> specCol(tab, "SPECTRA");
1494 //ArrayColumn<Float> tsysCol(tout, "TSYS");
1495 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1496 Vector<Float> spec;
1497 Vector<uChar> flag;
1498 specCol.get(0,spec); // the number of channels should be constant per IF
1499 uInt nChanIn = spec.nelements();
1500 Vector<Float> xIn(nChanIn); indgen(xIn);
1501 Int fac = Int(nChanIn/width);
1502 Vector<Float> xOut(fac+10); // 10 to be safe - resize later
1503 uInt k = 0;
1504 Float x = 0.0;
1505 while (x < Float(nChanIn) ) {
1506 xOut(k) = x;
1507 k++;
1508 x += width;
1509 }
1510 uInt nChanOut = k;
1511 xOut.resize(nChanOut, True);
1512 // process all rows for this IFNO
1513 Vector<Float> specOut;
1514 Vector<Bool> maskOut;
1515 Vector<uChar> flagOut;
1516 for (uInt i=0; i < tab.nrow(); ++i) {
1517 specCol.get(i, spec);
1518 flagCol.get(i, flag);
1519 Vector<Bool> mask(flag.nelements());
1520 convertArray(mask, flag);
1521
1522 IPosition shapeIn(spec.shape());
1523 //sh.nchan = nChanOut;
1524 InterpolateArray1D<Float,Float>::interpolate(specOut, maskOut, xOut,
1525 xIn, spec, mask,
1526 interpMethod, True, True);
1527 /// @todo do the same for channel based Tsys
1528 flagOut.resize(maskOut.nelements());
1529 convertArray(flagOut, maskOut);
1530 specCol.put(i, specOut);
1531 flagCol.put(i, flagOut);
1532 }
1533 ++iter;
1534 }
1535
1536 return out;
1537}
1538
1539STMath::imethod STMath::stringToIMethod(const std::string& in)
1540{
1541 static STMath::imap lookup;
1542
1543 // initialize the lookup table if necessary
1544 if ( lookup.empty() ) {
1545 lookup["nearest"] = InterpolateArray1D<Double,Float>::nearestNeighbour;
1546 lookup["linear"] = InterpolateArray1D<Double,Float>::linear;
1547 lookup["cubic"] = InterpolateArray1D<Double,Float>::cubic;
1548 lookup["spline"] = InterpolateArray1D<Double,Float>::spline;
1549 }
1550
1551 STMath::imap::const_iterator iter = lookup.find(in);
1552
1553 if ( lookup.end() == iter ) {
1554 std::string message = in;
1555 message += " is not a valid interpolation mode";
1556 throw(AipsError(message));
1557 }
1558 return iter->second;
1559}
1560
1561WeightType STMath::stringToWeight(const std::string& in)
1562{
1563 static std::map<std::string, WeightType> lookup;
1564
1565 // initialize the lookup table if necessary
1566 if ( lookup.empty() ) {
1567 lookup["NONE"] = asap::NONE;
1568 lookup["TINT"] = asap::TINT;
1569 lookup["TINTSYS"] = asap::TINTSYS;
1570 lookup["TSYS"] = asap::TSYS;
1571 lookup["VAR"] = asap::VAR;
1572 }
1573
1574 std::map<std::string, WeightType>::const_iterator iter = lookup.find(in);
1575
1576 if ( lookup.end() == iter ) {
1577 std::string message = in;
1578 message += " is not a valid weighting mode";
1579 throw(AipsError(message));
1580 }
1581 return iter->second;
1582}
1583
1584CountedPtr< Scantable > STMath::gainElevation( const CountedPtr< Scantable >& in,
1585 const vector< float > & coeff,
1586 const std::string & filename,
1587 const std::string& method)
1588{
1589 // Get elevation data from Scantable and convert to degrees
1590 CountedPtr< Scantable > out = getScantable(in, false);
1591 Table& tab = out->table();
1592 ROScalarColumn<Float> elev(tab, "ELEVATION");
1593 Vector<Float> x = elev.getColumn();
1594 x *= Float(180 / C::pi); // Degrees
1595
1596 Vector<Float> coeffs(coeff);
1597 const uInt nc = coeffs.nelements();
1598 if ( filename.length() > 0 && nc > 0 ) {
1599 throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
1600 }
1601
1602 // Correct
1603 if ( nc > 0 || filename.length() == 0 ) {
1604 // Find instrument
1605 Bool throwit = True;
1606 Instrument inst =
1607 STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
1608 throwit);
1609
1610 // Set polynomial
1611 Polynomial<Float>* ppoly = 0;
1612 Vector<Float> coeff;
1613 String msg;
1614 if ( nc > 0 ) {
1615 ppoly = new Polynomial<Float>(nc);
1616 coeff = coeffs;
1617 msg = String("user");
1618 } else {
1619 STAttr sdAttr;
1620 coeff = sdAttr.gainElevationPoly(inst);
1621 ppoly = new Polynomial<Float>(3);
1622 msg = String("built in");
1623 }
1624
1625 if ( coeff.nelements() > 0 ) {
1626 ppoly->setCoefficients(coeff);
1627 } else {
1628 delete ppoly;
1629 throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
1630 }
1631 ostringstream oss;
1632 oss << "Making polynomial correction with " << msg << " coefficients:" << endl;
1633 oss << " " << coeff;
1634 pushLog(String(oss));
1635 const uInt nrow = tab.nrow();
1636 Vector<Float> factor(nrow);
1637 for ( uInt i=0; i < nrow; ++i ) {
1638 factor[i] = 1.0 / (*ppoly)(x[i]);
1639 }
1640 delete ppoly;
1641 scaleByVector(tab, factor, true);
1642
1643 } else {
1644 // Read and correct
1645 pushLog("Making correction from ascii Table");
1646 scaleFromAsciiTable(tab, filename, method, x, true);
1647 }
1648 return out;
1649}
1650
1651void STMath::scaleFromAsciiTable(Table& in, const std::string& filename,
1652 const std::string& method,
1653 const Vector<Float>& xout, bool dotsys)
1654{
1655
1656// Read gain-elevation ascii file data into a Table.
1657
1658 String formatString;
1659 Table tbl = readAsciiTable(formatString, Table::Memory, filename, "", "", False);
1660 scaleFromTable(in, tbl, method, xout, dotsys);
1661}
1662
1663void STMath::scaleFromTable(Table& in,
1664 const Table& table,
1665 const std::string& method,
1666 const Vector<Float>& xout, bool dotsys)
1667{
1668
1669 ROScalarColumn<Float> geElCol(table, "ELEVATION");
1670 ROScalarColumn<Float> geFacCol(table, "FACTOR");
1671 Vector<Float> xin = geElCol.getColumn();
1672 Vector<Float> yin = geFacCol.getColumn();
1673 Vector<Bool> maskin(xin.nelements(),True);
1674
1675 // Interpolate (and extrapolate) with desired method
1676
1677 InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
1678
1679 Vector<Float> yout;
1680 Vector<Bool> maskout;
1681 InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout,
1682 xin, yin, maskin, interp,
1683 True, True);
1684
1685 scaleByVector(in, Float(1.0)/yout, dotsys);
1686}
1687
1688void STMath::scaleByVector( Table& in,
1689 const Vector< Float >& factor,
1690 bool dotsys )
1691{
1692 uInt nrow = in.nrow();
1693 if ( factor.nelements() != nrow ) {
1694 throw(AipsError("factors.nelements() != table.nelements()"));
1695 }
1696 ArrayColumn<Float> specCol(in, "SPECTRA");
1697 ArrayColumn<uChar> flagCol(in, "FLAGTRA");
1698 ArrayColumn<Float> tsysCol(in, "TSYS");
1699 for (uInt i=0; i < nrow; ++i) {
1700 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
1701 ma *= factor[i];
1702 specCol.put(i, ma.getArray());
1703 flagCol.put(i, flagsFromMA(ma));
1704 if ( dotsys ) {
1705 Vector<Float> tsys = tsysCol(i);
1706 tsys *= factor[i];
1707 tsysCol.put(i,tsys);
1708 }
1709 }
1710}
1711
1712CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in,
1713 float d, float etaap,
1714 float jyperk )
1715{
1716 CountedPtr< Scantable > out = getScantable(in, false);
1717 Table& tab = in->table();
1718 Unit fluxUnit(tab.keywordSet().asString("FluxUnit"));
1719 Unit K(String("K"));
1720 Unit JY(String("Jy"));
1721
1722 bool tokelvin = true;
1723 Double cfac = 1.0;
1724
1725 if ( fluxUnit == JY ) {
1726 pushLog("Converting to K");
1727 Quantum<Double> t(1.0,fluxUnit);
1728 Quantum<Double> t2 = t.get(JY);
1729 cfac = (t2 / t).getValue(); // value to Jy
1730
1731 tokelvin = true;
1732 out->setFluxUnit("K");
1733 } else if ( fluxUnit == K ) {
1734 pushLog("Converting to Jy");
1735 Quantum<Double> t(1.0,fluxUnit);
1736 Quantum<Double> t2 = t.get(K);
1737 cfac = (t2 / t).getValue(); // value to K
1738
1739 tokelvin = false;
1740 out->setFluxUnit("Jy");
1741 } else {
1742 throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
1743 }
1744 // Make sure input values are converted to either Jy or K first...
1745 Float factor = cfac;
1746
1747 // Select method
1748 if (jyperk > 0.0) {
1749 factor *= jyperk;
1750 if ( tokelvin ) factor = 1.0 / jyperk;
1751 ostringstream oss;
1752 oss << "Jy/K = " << jyperk;
1753 pushLog(String(oss));
1754 Vector<Float> factors(tab.nrow(), factor);
1755 scaleByVector(tab,factors, false);
1756 } else if ( etaap > 0.0) {
1757 if (d < 0) {
1758 Instrument inst =
1759 STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
1760 True);
1761 STAttr sda;
1762 d = sda.diameter(inst);
1763 }
1764 jyperk = STAttr::findJyPerK(etaap, d);
1765 ostringstream oss;
1766 oss << "Jy/K = " << jyperk;
1767 pushLog(String(oss));
1768 factor *= jyperk;
1769 if ( tokelvin ) {
1770 factor = 1.0 / factor;
1771 }
1772 Vector<Float> factors(tab.nrow(), factor);
1773 scaleByVector(tab, factors, False);
1774 } else {
1775
1776 // OK now we must deal with automatic look up of values.
1777 // We must also deal with the fact that the factors need
1778 // to be computed per IF and may be different and may
1779 // change per integration.
1780
1781 pushLog("Looking up conversion factors");
1782 convertBrightnessUnits(out, tokelvin, cfac);
1783 }
1784
1785 return out;
1786}
1787
1788void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in,
1789 bool tokelvin, float cfac )
1790{
1791 Table& table = in->table();
1792 Instrument inst =
1793 STAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True);
1794 TableIterator iter(table, "FREQ_ID");
1795 STFrequencies stfreqs = in->frequencies();
1796 STAttr sdAtt;
1797 while (!iter.pastEnd()) {
1798 Table tab = iter.table();
1799 ArrayColumn<Float> specCol(tab, "SPECTRA");
1800 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1801 ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID");
1802 MEpoch::ROScalarColumn timeCol(tab, "TIME");
1803
1804 uInt freqid; freqidCol.get(0, freqid);
1805 Vector<Float> tmpspec; specCol.get(0, tmpspec);
1806 // STAttr.JyPerK has a Vector interface... change sometime.
1807 Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements()));
1808 for ( uInt i=0; i<tab.nrow(); ++i) {
1809 Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0];
1810 Float factor = cfac * jyperk;
1811 if ( tokelvin ) factor = Float(1.0) / factor;
1812 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
1813 ma *= factor;
1814 specCol.put(i, ma.getArray());
1815 flagCol.put(i, flagsFromMA(ma));
1816 }
1817 ++iter;
1818 }
1819}
1820
1821CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
1822 float tau )
1823{
1824 CountedPtr< Scantable > out = getScantable(in, false);
1825
1826 Table tab = out->table();
1827 ROScalarColumn<Float> elev(tab, "ELEVATION");
1828 ArrayColumn<Float> specCol(tab, "SPECTRA");
1829 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1830 ArrayColumn<Float> tsysCol(tab, "TSYS");
1831 for ( uInt i=0; i<tab.nrow(); ++i) {
1832 Float zdist = Float(C::pi_2) - elev(i);
1833 Float factor = exp(tau/cos(zdist));
1834 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
1835 ma *= factor;
1836 specCol.put(i, ma.getArray());
1837 flagCol.put(i, flagsFromMA(ma));
1838 Vector<Float> tsys;
1839 tsysCol.get(i, tsys);
1840 tsys *= factor;
1841 tsysCol.put(i, tsys);
1842 }
1843 return out;
1844}
1845
1846CountedPtr< Scantable > STMath::smoothOther( const CountedPtr< Scantable >& in,
1847 const std::string& kernel,
1848 float width )
1849{
1850 CountedPtr< Scantable > out = getScantable(in, false);
1851 Table& table = out->table();
1852 ArrayColumn<Float> specCol(table, "SPECTRA");
1853 ArrayColumn<uChar> flagCol(table, "FLAGTRA");
1854 Vector<Float> spec;
1855 Vector<uChar> flag;
1856 for ( uInt i=0; i<table.nrow(); ++i) {
1857 specCol.get(i, spec);
1858 flagCol.get(i, flag);
1859 Vector<Bool> mask(flag.nelements());
1860 convertArray(mask, flag);
1861 Vector<Float> specout;
1862 Vector<Bool> maskout;
1863 if ( kernel == "hanning" ) {
1864 mathutil::hanning(specout, maskout, spec , !mask);
1865 convertArray(flag, !maskout);
1866 } else if ( kernel == "rmedian" ) {
1867 mathutil::runningMedian(specout, maskout, spec , mask, width);
1868 convertArray(flag, maskout);
1869 }
1870 flagCol.put(i, flag);
1871 specCol.put(i, specout);
1872 }
1873 return out;
1874}
1875
1876CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
1877 const std::string& kernel, float width )
1878{
1879 if (kernel == "rmedian" || kernel == "hanning") {
1880 return smoothOther(in, kernel, width);
1881 }
1882 CountedPtr< Scantable > out = getScantable(in, false);
1883 Table& table = out->table();
1884 VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
1885 // same IFNO should have same no of channels
1886 // this saves overhead
1887 TableIterator iter(table, "IFNO");
1888 while (!iter.pastEnd()) {
1889 Table tab = iter.table();
1890 ArrayColumn<Float> specCol(tab, "SPECTRA");
1891 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1892 Vector<Float> tmpspec; specCol.get(0, tmpspec);
1893 uInt nchan = tmpspec.nelements();
1894 Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
1895 Convolver<Float> conv(kvec, IPosition(1,nchan));
1896 Vector<Float> spec;
1897 Vector<uChar> flag;
1898 for ( uInt i=0; i<tab.nrow(); ++i) {
1899 specCol.get(i, spec);
1900 flagCol.get(i, flag);
1901 Vector<Bool> mask(flag.nelements());
1902 convertArray(mask, flag);
1903 Vector<Float> specout;
1904 mathutil::replaceMaskByZero(specout, mask);
1905 conv.linearConv(specout, spec);
1906 specCol.put(i, specout);
1907 }
1908 ++iter;
1909 }
1910 return out;
1911}
1912
1913CountedPtr< Scantable >
1914 STMath::merge( const std::vector< CountedPtr < Scantable > >& in )
1915{
1916 if ( in.size() < 2 ) {
1917 throw(AipsError("Need at least two scantables to perform a merge."));
1918 }
1919 std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
1920 bool insitu = insitu_;
1921 setInsitu(false);
1922 CountedPtr< Scantable > out = getScantable(*it, false);
1923 setInsitu(insitu);
1924 Table& tout = out->table();
1925 ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
1926 ScalarColumn<uInt> scannocol(tout,"SCANNO"), focusidcol(tout,"FOCUS_ID");
1927 // Renumber SCANNO to be 0-based
1928 Vector<uInt> scannos = scannocol.getColumn();
1929 uInt offset = min(scannos);
1930 scannos -= offset;
1931 scannocol.putColumn(scannos);
1932 uInt newscanno = max(scannos)+1;
1933 ++it;
1934 while ( it != in.end() ){
1935 if ( ! (*it)->conformant(*out) ) {
1936 // non conformant.
1937 pushLog(String("Warning: Can't merge scantables as header info differs."));
1938 }
1939 out->appendToHistoryTable((*it)->history());
1940 const Table& tab = (*it)->table();
1941 TableIterator scanit(tab, "SCANNO");
1942 while (!scanit.pastEnd()) {
1943 TableIterator freqit(scanit.table(), "FREQ_ID");
1944 while ( !freqit.pastEnd() ) {
1945 Table thetab = freqit.table();
1946 uInt nrow = tout.nrow();
1947 tout.addRow(thetab.nrow());
1948 TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
1949 ROTableRow row(thetab);
1950 for ( uInt i=0; i<thetab.nrow(); ++i) {
1951 uInt k = nrow+i;
1952 scannocol.put(k, newscanno);
1953 const TableRecord& rec = row.get(i);
1954 Double rv,rp,inc;
1955 (*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
1956 uInt id;
1957 id = out->frequencies().addEntry(rp, rv, inc);
1958 freqidcol.put(k,id);
1959 //String name,fname;Double rf;
1960 Vector<String> name,fname;Vector<Double> rf;
1961 (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
1962 id = out->molecules().addEntry(rf, name, fname);
1963 molidcol.put(k, id);
1964 Float frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
1965 (*it)->focus().getEntry(fax, ftan, frot, fhand,
1966 fmount,fuser, fxy, fxyp,
1967 rec.asuInt("FOCUS_ID"));
1968 id = out->focus().addEntry(fax, ftan, frot, fhand,
1969 fmount,fuser, fxy, fxyp);
1970 focusidcol.put(k, id);
1971 }
1972 ++freqit;
1973 }
1974 ++newscanno;
1975 ++scanit;
1976 }
1977 ++it;
1978 }
1979 return out;
1980}
1981
1982CountedPtr< Scantable >
1983 STMath::invertPhase( const CountedPtr < Scantable >& in )
1984{
1985 return applyToPol(in, &STPol::invertPhase, Float(0.0));
1986}
1987
1988CountedPtr< Scantable >
1989 STMath::rotateXYPhase( const CountedPtr < Scantable >& in, float phase )
1990{
1991 return applyToPol(in, &STPol::rotatePhase, Float(phase));
1992}
1993
1994CountedPtr< Scantable >
1995 STMath::rotateLinPolPhase( const CountedPtr < Scantable >& in, float phase )
1996{
1997 return applyToPol(in, &STPol::rotateLinPolPhase, Float(phase));
1998}
1999
2000CountedPtr< Scantable > STMath::applyToPol( const CountedPtr<Scantable>& in,
2001 STPol::polOperation fptr,
2002 Float phase )
2003{
2004 CountedPtr< Scantable > out = getScantable(in, false);
2005 Table& tout = out->table();
2006 Block<String> cols(4);
2007 cols[0] = String("SCANNO");
2008 cols[1] = String("BEAMNO");
2009 cols[2] = String("IFNO");
2010 cols[3] = String("CYCLENO");
2011 TableIterator iter(tout, cols);
2012 CountedPtr<STPol> stpol = STPol::getPolClass(out->factories_,
2013 out->getPolType() );
2014 while (!iter.pastEnd()) {
2015 Table t = iter.table();
2016 ArrayColumn<Float> speccol(t, "SPECTRA");
2017 ScalarColumn<uInt> focidcol(t, "FOCUS_ID");
2018 ScalarColumn<Float> parancol(t, "PARANGLE");
2019 Matrix<Float> pols(speccol.getColumn());
2020 try {
2021 stpol->setSpectra(pols);
2022 Float fang,fhand,parang;
2023 fang = in->focusTable_.getTotalFeedAngle(focidcol(0));
2024 fhand = in->focusTable_.getFeedHand(focidcol(0));
2025 parang = parancol(0);
2026 /// @todo re-enable this
2027 // disable total feed angle to support paralactifying Caswell style
2028 stpol->setPhaseCorrections(parang, -parang, fhand);
2029 // use a member function pointer in STPol. This only works on
2030 // the STPol pointer itself, not the Counted Pointer so
2031 // derefernce it.
2032 (&(*(stpol))->*fptr)(phase);
2033 speccol.putColumn(stpol->getSpectra());
2034 } catch (AipsError& e) {
2035 //delete stpol;stpol=0;
2036 throw(e);
2037 }
2038 ++iter;
2039 }
2040 //delete stpol;stpol=0;
2041 return out;
2042}
2043
2044CountedPtr< Scantable >
2045 STMath::swapPolarisations( const CountedPtr< Scantable > & in )
2046{
2047 CountedPtr< Scantable > out = getScantable(in, false);
2048 Table& tout = out->table();
2049 Table t0 = tout(tout.col("POLNO") == 0);
2050 Table t1 = tout(tout.col("POLNO") == 1);
2051 if ( t0.nrow() != t1.nrow() )
2052 throw(AipsError("Inconsistent number of polarisations"));
2053 ArrayColumn<Float> speccol0(t0, "SPECTRA");
2054 ArrayColumn<uChar> flagcol0(t0, "FLAGTRA");
2055 ArrayColumn<Float> speccol1(t1, "SPECTRA");
2056 ArrayColumn<uChar> flagcol1(t1, "FLAGTRA");
2057 Matrix<Float> s0 = speccol0.getColumn();
2058 Matrix<uChar> f0 = flagcol0.getColumn();
2059 speccol0.putColumn(speccol1.getColumn());
2060 flagcol0.putColumn(flagcol1.getColumn());
2061 speccol1.putColumn(s0);
2062 flagcol1.putColumn(f0);
2063 return out;
2064}
2065
2066CountedPtr< Scantable >
2067 STMath::averagePolarisations( const CountedPtr< Scantable > & in,
2068 const std::vector<bool>& mask,
2069 const std::string& weight )
2070{
2071 if (in->npol() < 2 )
2072 throw(AipsError("averagePolarisations can only be applied to two or more"
2073 "polarisations"));
2074 bool insitu = insitu_;
2075 setInsitu(false);
2076 CountedPtr< Scantable > pols = getScantable(in, true);
2077 setInsitu(insitu);
2078 Table& tout = pols->table();
2079 std::string taql = "SELECT FROM $1 WHERE POLNO IN [0,1]";
2080 Table tab = tableCommand(taql, in->table());
2081 if (tab.nrow() == 0 )
2082 throw(AipsError("Could not find any rows with POLNO==0 and POLNO==1"));
2083 TableCopy::copyRows(tout, tab);
2084 TableVector<uInt> vec(tout, "POLNO");
2085 vec = 0;
2086 pols->table_.rwKeywordSet().define("nPol", Int(1));
2087 //pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
2088 pols->table_.rwKeywordSet().define("POLTYPE", in->getPolType());
2089 std::vector<CountedPtr<Scantable> > vpols;
2090 vpols.push_back(pols);
2091 CountedPtr< Scantable > out = average(vpols, mask, weight, "SCAN");
2092 return out;
2093}
2094
2095CountedPtr< Scantable >
2096 STMath::averageBeams( const CountedPtr< Scantable > & in,
2097 const std::vector<bool>& mask,
2098 const std::string& weight )
2099{
2100 bool insitu = insitu_;
2101 setInsitu(false);
2102 CountedPtr< Scantable > beams = getScantable(in, false);
2103 setInsitu(insitu);
2104 Table& tout = beams->table();
2105 // give all rows the same BEAMNO
2106 TableVector<uInt> vec(tout, "BEAMNO");
2107 vec = 0;
2108 beams->table_.rwKeywordSet().define("nBeam", Int(1));
2109 std::vector<CountedPtr<Scantable> > vbeams;
2110 vbeams.push_back(beams);
2111 CountedPtr< Scantable > out = average(vbeams, mask, weight, "SCAN");
2112 return out;
2113}
2114
2115
2116CountedPtr< Scantable >
2117 asap::STMath::frequencyAlign( const CountedPtr< Scantable > & in,
2118 const std::string & refTime,
2119 const std::string & method)
2120{
2121 // clone as this is not working insitu
2122 bool insitu = insitu_;
2123 setInsitu(false);
2124 CountedPtr< Scantable > out = getScantable(in, false);
2125 setInsitu(insitu);
2126 Table& tout = out->table();
2127 // Get reference Epoch to time of first row or given String
2128 Unit DAY(String("d"));
2129 MEpoch::Ref epochRef(in->getTimeReference());
2130 MEpoch refEpoch;
2131 if (refTime.length()>0) {
2132 Quantum<Double> qt;
2133 if (MVTime::read(qt,refTime)) {
2134 MVEpoch mv(qt);
2135 refEpoch = MEpoch(mv, epochRef);
2136 } else {
2137 throw(AipsError("Invalid format for Epoch string"));
2138 }
2139 } else {
2140 refEpoch = in->timeCol_(0);
2141 }
2142 MPosition refPos = in->getAntennaPosition();
2143
2144 InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
2145 /*
2146 // Comment from MV.
2147 // the following code has been commented out because different FREQ_IDs have to be aligned together even
2148 // if the frame doesn't change. So far, lack of this check didn't cause any problems.
2149 // test if user frame is different to base frame
2150 if ( in->frequencies().getFrameString(true)
2151 == in->frequencies().getFrameString(false) ) {
2152 throw(AipsError("Can't convert as no output frame has been set"
2153 " (use set_freqframe) or it is aligned already."));
2154 }
2155 */
2156 MFrequency::Types system = in->frequencies().getFrame();
2157 MVTime mvt(refEpoch.getValue());
2158 String epochout = mvt.string(MVTime::YMD) + String(" (") + refEpoch.getRefString() + String(")");
2159 ostringstream oss;
2160 oss << "Aligned at reference Epoch " << epochout
2161 << " in frame " << MFrequency::showType(system);
2162 pushLog(String(oss));
2163 // set up the iterator
2164 Block<String> cols(4);
2165 // select by constant direction
2166 cols[0] = String("SRCNAME");
2167 cols[1] = String("BEAMNO");
2168 // select by IF ( no of channels varies over this )
2169 cols[2] = String("IFNO");
2170 // select by restfrequency
2171 cols[3] = String("MOLECULE_ID");
2172 TableIterator iter(tout, cols);
2173 while ( !iter.pastEnd() ) {
2174 Table t = iter.table();
2175 MDirection::ROScalarColumn dirCol(t, "DIRECTION");
2176 TableIterator fiter(t, "FREQ_ID");
2177 // determine nchan from the first row. This should work as
2178 // we are iterating over BEAMNO and IFNO // we should have constant direction
2179
2180 ROArrayColumn<Float> sCol(t, "SPECTRA");
2181 const MDirection direction = dirCol(0);
2182 const uInt nchan = sCol(0).nelements();
2183
2184 // skip operations if there is nothing to align
2185 if (fiter.pastEnd()) {
2186 continue;
2187 }
2188
2189 Table ftab = fiter.table();
2190 // align all frequency ids with respect to the first encountered id
2191 ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
2192 // get the SpectralCoordinate for the freqid, which we are iterating over
2193 SpectralCoordinate sC = in->frequencies().getSpectralCoordinate(freqidCol(0));
2194 FrequencyAligner<Float> fa( sC, nchan, refEpoch,
2195 direction, refPos, system );
2196 // realign the SpectralCoordinate and put into the output Scantable
2197 Vector<String> units(1);
2198 units = String("Hz");
2199 Bool linear=True;
2200 SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
2201 sc2.setWorldAxisUnits(units);
2202 const uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
2203 sc2.referenceValue()[0],
2204 sc2.increment()[0]);
2205 while ( !fiter.pastEnd() ) {
2206 ftab = fiter.table();
2207 // spectral coordinate for the current FREQ_ID
2208 ScalarColumn<uInt> freqidCol2(ftab, "FREQ_ID");
2209 sC = in->frequencies().getSpectralCoordinate(freqidCol2(0));
2210 // create the "global" abcissa for alignment with same FREQ_ID
2211 Vector<Double> abc(nchan);
2212 for (uInt i=0; i<nchan; i++) {
2213 Double w;
2214 sC.toWorld(w,Double(i));
2215 abc[i] = w;
2216 }
2217 TableVector<uInt> tvec(ftab, "FREQ_ID");
2218 // assign new frequency id to all rows
2219 tvec = id;
2220 // cache abcissa for same time stamps, so iterate over those
2221 TableIterator timeiter(ftab, "TIME");
2222 while ( !timeiter.pastEnd() ) {
2223 Table tab = timeiter.table();
2224 ArrayColumn<Float> specCol(tab, "SPECTRA");
2225 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2226 MEpoch::ROScalarColumn timeCol(tab, "TIME");
2227 // use align abcissa cache after the first row
2228 // these rows should be just be POLNO
2229 bool first = true;
2230 for (int i=0; i<int(tab.nrow()); ++i) {
2231 // input values
2232 Vector<uChar> flag = flagCol(i);
2233 Vector<Bool> mask(flag.shape());
2234 Vector<Float> specOut, spec;
2235 spec = specCol(i);
2236 Vector<Bool> maskOut;Vector<uChar> flagOut;
2237 convertArray(mask, flag);
2238 // alignment
2239 Bool ok = fa.align(specOut, maskOut, abc, spec,
2240 mask, timeCol(i), !first,
2241 interp, False);
2242 // back into scantable
2243 flagOut.resize(maskOut.nelements());
2244 convertArray(flagOut, maskOut);
2245 flagCol.put(i, flagOut);
2246 specCol.put(i, specOut);
2247 // start abcissa caching
2248 first = false;
2249 }
2250 // next timestamp
2251 ++timeiter;
2252 }
2253 // next FREQ_ID
2254 ++fiter;
2255 }
2256 // next aligner
2257 ++iter;
2258 }
2259 // set this afterwards to ensure we are doing insitu correctly.
2260 out->frequencies().setFrame(system, true);
2261 return out;
2262}
2263
2264CountedPtr<Scantable>
2265 asap::STMath::convertPolarisation( const CountedPtr<Scantable>& in,
2266 const std::string & newtype )
2267{
2268 if (in->npol() != 2 && in->npol() != 4)
2269 throw(AipsError("Can only convert two or four polarisations."));
2270 if ( in->getPolType() == newtype )
2271 throw(AipsError("No need to convert."));
2272 if ( ! in->selector_.empty() )
2273 throw(AipsError("Can only convert whole scantable. Unset the selection."));
2274 bool insitu = insitu_;
2275 setInsitu(false);
2276 CountedPtr< Scantable > out = getScantable(in, true);
2277 setInsitu(insitu);
2278 Table& tout = out->table();
2279 tout.rwKeywordSet().define("POLTYPE", String(newtype));
2280
2281 Block<String> cols(4);
2282 cols[0] = "SCANNO";
2283 cols[1] = "CYCLENO";
2284 cols[2] = "BEAMNO";
2285 cols[3] = "IFNO";
2286 TableIterator it(in->originalTable_, cols);
2287 String basetype = in->getPolType();
2288 STPol* stpol = STPol::getPolClass(in->factories_, basetype);
2289 try {
2290 while ( !it.pastEnd() ) {
2291 Table tab = it.table();
2292 uInt row = tab.rowNumbers()[0];
2293 stpol->setSpectra(in->getPolMatrix(row));
2294 Float fang,fhand,parang;
2295 fang = in->focusTable_.getTotalFeedAngle(in->mfocusidCol_(row));
2296 fhand = in->focusTable_.getFeedHand(in->mfocusidCol_(row));
2297 parang = in->paraCol_(row);
2298 /// @todo re-enable this
2299 // disable total feed angle to support paralactifying Caswell style
2300 stpol->setPhaseCorrections(parang, -parang, fhand);
2301 Int npolout = 0;
2302 for (uInt i=0; i<tab.nrow(); ++i) {
2303 Vector<Float> outvec = stpol->getSpectrum(i, newtype);
2304 if ( outvec.nelements() > 0 ) {
2305 tout.addRow();
2306 TableCopy::copyRows(tout, tab, tout.nrow()-1, 0, 1);
2307 ArrayColumn<Float> sCol(tout,"SPECTRA");
2308 ScalarColumn<uInt> pCol(tout,"POLNO");
2309 sCol.put(tout.nrow()-1 ,outvec);
2310 pCol.put(tout.nrow()-1 ,uInt(npolout));
2311 npolout++;
2312 }
2313 }
2314 tout.rwKeywordSet().define("nPol", npolout);
2315 ++it;
2316 }
2317 } catch (AipsError& e) {
2318 delete stpol;
2319 throw(e);
2320 }
2321 delete stpol;
2322 return out;
2323}
2324
2325CountedPtr< Scantable >
2326 asap::STMath::mxExtract( const CountedPtr< Scantable > & in,
2327 const std::string & scantype )
2328{
2329 bool insitu = insitu_;
2330 setInsitu(false);
2331 CountedPtr< Scantable > out = getScantable(in, true);
2332 setInsitu(insitu);
2333 Table& tout = out->table();
2334 std::string taql = "SELECT FROM $1 WHERE BEAMNO != REFBEAMNO";
2335 if (scantype == "on") {
2336 taql = "SELECT FROM $1 WHERE BEAMNO == REFBEAMNO";
2337 }
2338 Table tab = tableCommand(taql, in->table());
2339 TableCopy::copyRows(tout, tab);
2340 if (scantype == "on") {
2341 // re-index SCANNO to 0
2342 TableVector<uInt> vec(tout, "SCANNO");
2343 vec = 0;
2344 }
2345 return out;
2346}
2347
2348CountedPtr< Scantable >
2349 asap::STMath::lagFlag( const CountedPtr< Scantable > & in,
2350 double frequency, double width )
2351{
2352 CountedPtr< Scantable > out = getScantable(in, false);
2353 Table& tout = out->table();
2354 TableIterator iter(tout, "FREQ_ID");
2355 FFTServer<Float,Complex> ffts;
2356 while ( !iter.pastEnd() ) {
2357 Table tab = iter.table();
2358 Double rp,rv,inc;
2359 ROTableRow row(tab);
2360 const TableRecord& rec = row.get(0);
2361 uInt freqid = rec.asuInt("FREQ_ID");
2362 out->frequencies().getEntry(rp, rv, inc, freqid);
2363 ArrayColumn<Float> specCol(tab, "SPECTRA");
2364 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2365 for (int i=0; i<int(tab.nrow()); ++i) {
2366 Vector<Float> spec = specCol(i);
2367 Vector<uChar> flag = flagCol(i);
2368 Int lag0 = Int(spec.nelements()*abs(inc)/(frequency+width)+0.5);
2369 Int lag1 = Int(spec.nelements()*abs(inc)/(frequency-width)+0.5);
2370 for (int k=0; k < flag.nelements(); ++k ) {
2371 if (flag[k] > 0) {
2372 spec[k] = 0.0;
2373 }
2374 }
2375 Vector<Complex> lags;
2376 ffts.fft0(lags, spec);
2377 Int start = max(0, lag0);
2378 Int end = min(Int(lags.nelements()-1), lag1);
2379 if (start == end) {
2380 lags[start] = Complex(0.0);
2381 } else {
2382 for (int j=start; j <=end ;++j) {
2383 lags[j] = Complex(0.0);
2384 }
2385 }
2386 ffts.fft0(spec, lags);
2387 specCol.put(i, spec);
2388 }
2389 ++iter;
2390 }
2391 return out;
2392}
2393
2394// Averaging spectra with different channel/resolution
2395CountedPtr<Scantable>
2396STMath::new_average( const std::vector<CountedPtr<Scantable> >& in,
2397 const bool& compel,
2398 const std::vector<bool>& mask,
2399 const std::string& weight,
2400 const std::string& avmode )
2401 throw ( casa::AipsError )
2402{
2403 if ( avmode == "SCAN" && in.size() != 1 )
2404 throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
2405 "Use merge first."));
2406
2407 // check if OTF observation
2408 String obstype = in[0]->getHeader().obstype ;
2409 bool otfscan = false ;
2410 if ( obstype.find( "OTF" ) != String::npos ) {
2411 //cout << "OTF scan" << endl ;
2412 otfscan = true ;
2413 }
2414
2415 CountedPtr<Scantable> out ; // processed result
2416 if ( compel ) {
2417 std::vector< CountedPtr<Scantable> > newin ; // input for average process
2418 uInt insize = in.size() ; // number of input scantables
2419
2420 // TEST: do normal average in each table before IF grouping
2421 cout << "Do preliminary averaging" << endl ;
2422 vector< CountedPtr<Scantable> > tmpin( insize ) ;
2423 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2424 vector< CountedPtr<Scantable> > v( 1, in[itable] ) ;
2425 tmpin[itable] = average( v, mask, weight, avmode ) ;
2426 }
2427
2428 // warning
2429 cout << "Average spectra with different spectral resolution" << endl ;
2430 cout << endl ;
2431
2432 // temporarily set coordinfo
2433 vector<string> oldinfo( insize ) ;
2434 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2435 vector<string> coordinfo = in[itable]->getCoordInfo() ;
2436 oldinfo[itable] = coordinfo[0] ;
2437 coordinfo[0] = "Hz" ;
2438 tmpin[itable]->setCoordInfo( coordinfo ) ;
2439 }
2440
2441 // columns
2442 ScalarColumn<uInt> freqIDCol ;
2443 ScalarColumn<uInt> ifnoCol ;
2444 ScalarColumn<uInt> scannoCol ;
2445
2446
2447 // check IF frequency coverage
2448 // freqid: list of FREQ_ID, which is used, in each table
2449 // iffreq: list of minimum and maximum frequency for each FREQ_ID in
2450 // each table
2451 // freqid[insize][numIF]
2452 // freqid: [[id00, id01, ...],
2453 // [id10, id11, ...],
2454 // ...
2455 // [idn0, idn1, ...]]
2456 // iffreq[insize][numIF*2]
2457 // iffreq: [[min_id00, max_id00, min_id01, max_id01, ...],
2458 // [min_id10, max_id10, min_id11, max_id11, ...],
2459 // ...
2460 // [min_idn0, max_idn0, min_idn1, max_idn1, ...]]
2461 //cout << "Check IF settings in each table" << endl ;
2462 vector< vector<uInt> > freqid( insize );
2463 vector< vector<double> > iffreq( insize ) ;
2464 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2465 uInt rows = tmpin[itable]->nrow() ;
2466 uInt freqnrows = tmpin[itable]->frequencies().table().nrow() ;
2467 for ( uInt irow = 0 ; irow < rows ; irow++ ) {
2468 if ( freqid[itable].size() == freqnrows ) {
2469 break ;
2470 }
2471 else {
2472 freqIDCol.attach( tmpin[itable]->table(), "FREQ_ID" ) ;
2473 ifnoCol.attach( tmpin[itable]->table(), "IFNO" ) ;
2474 uInt id = freqIDCol( irow ) ;
2475 if ( freqid[itable].size() == 0 || count( freqid[itable].begin(), freqid[itable].end(), id ) == 0 ) {
2476 //cout << "itable = " << itable << ": IF " << id << " is included in the list" << endl ;
2477 vector<double> abcissa = tmpin[itable]->getAbcissa( irow ) ;
2478 freqid[itable].push_back( id ) ;
2479 iffreq[itable].push_back( abcissa[0] - 0.5 * ( abcissa[1] - abcissa[0] ) ) ;
2480 iffreq[itable].push_back( abcissa[abcissa.size()-1] + 0.5 * ( abcissa[1] - abcissa[0] ) ) ;
2481 }
2482 }
2483 }
2484 }
2485
2486 // debug
2487 //cout << "IF settings summary:" << endl ;
2488 //for ( uInt i = 0 ; i < freqid.size() ; i++ ) {
2489 //cout << " Table" << i << endl ;
2490 //for ( uInt j = 0 ; j < freqid[i].size() ; j++ ) {
2491 //cout << " id = " << freqid[i][j] << " (min,max) = (" << iffreq[i][2*j] << "," << iffreq[i][2*j+1] << ")" << endl ;
2492 //}
2493 //}
2494 //cout << endl ;
2495
2496 // IF grouping based on their frequency coverage
2497 // ifgrp: list of table index and FREQ_ID for all members in each IF group
2498 // ifgfreq: list of minimum and maximum frequency in each IF group
2499 // ifgrp[numgrp][nummember*2]
2500 // ifgrp: [[table00, freqrow00, table01, freqrow01, ...],
2501 // [table10, freqrow10, table11, freqrow11, ...],
2502 // ...
2503 // [tablen0, freqrown0, tablen1, freqrown1, ...]]
2504 // ifgfreq[numgrp*2]
2505 // ifgfreq: [min0_grp0, max0_grp0, min1_grp1, max1_grp1, ...]
2506 //cout << "IF grouping based on their frequency coverage" << endl ;
2507 vector< vector<uInt> > ifgrp ;
2508 vector<double> ifgfreq ;
2509
2510 // parameter for IF grouping
2511 // groupmode = OR retrieve all region
2512 // AND only retrieve overlaped region
2513 //string groupmode = "AND" ;
2514 string groupmode = "OR" ;
2515 uInt sizecr = 0 ;
2516 if ( groupmode == "AND" )
2517 sizecr = 2 ;
2518 else if ( groupmode == "OR" )
2519 sizecr = 0 ;
2520
2521 vector<double> sortedfreq ;
2522 for ( uInt i = 0 ; i < iffreq.size() ; i++ ) {
2523 for ( uInt j = 0 ; j < iffreq[i].size() ; j++ ) {
2524 if ( count( sortedfreq.begin(), sortedfreq.end(), iffreq[i][j] ) == 0 )
2525 sortedfreq.push_back( iffreq[i][j] ) ;
2526 }
2527 }
2528 sort( sortedfreq.begin(), sortedfreq.end() ) ;
2529 for ( vector<double>::iterator i = sortedfreq.begin() ; i != sortedfreq.end()-1 ; i++ ) {
2530 ifgfreq.push_back( *i ) ;
2531 ifgfreq.push_back( *(i+1) ) ;
2532 }
2533 ifgrp.resize( ifgfreq.size()/2 ) ;
2534 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2535 for ( uInt iif = 0 ; iif < freqid[itable].size() ; iif++ ) {
2536 double range0 = iffreq[itable][2*iif] ;
2537 double range1 = iffreq[itable][2*iif+1] ;
2538 for ( uInt j = 0 ; j < ifgrp.size() ; j++ ) {
2539 double fmin = max( range0, ifgfreq[2*j] ) ;
2540 double fmax = min( range1, ifgfreq[2*j+1] ) ;
2541 if ( fmin < fmax ) {
2542 ifgrp[j].push_back( itable ) ;
2543 ifgrp[j].push_back( freqid[itable][iif] ) ;
2544 }
2545 }
2546 }
2547 }
2548 vector< vector<uInt> >::iterator fiter = ifgrp.begin() ;
2549 vector<double>::iterator giter = ifgfreq.begin() ;
2550 while( fiter != ifgrp.end() ) {
2551 if ( fiter->size() <= sizecr ) {
2552 fiter = ifgrp.erase( fiter ) ;
2553 giter = ifgfreq.erase( giter ) ;
2554 giter = ifgfreq.erase( giter ) ;
2555 }
2556 else {
2557 fiter++ ;
2558 advance( giter, 2 ) ;
2559 }
2560 }
2561
2562 // Grouping continuous IF groups (without frequency gap)
2563 // freqgrp: list of IF group indexes in each frequency group
2564 // freqrange: list of minimum and maximum frequency in each frequency group
2565 // freqgrp[numgrp][nummember]
2566 // freqgrp: [[ifgrp00, ifgrp01, ifgrp02, ...],
2567 // [ifgrp10, ifgrp11, ifgrp12, ...],
2568 // ...
2569 // [ifgrpn0, ifgrpn1, ifgrpn2, ...]]
2570 // freqrange[numgrp*2]
2571 // freqrange: [min_grp0, max_grp0, min_grp1, max_grp1, ...]
2572 vector< vector<uInt> > freqgrp ;
2573 double freqrange = 0.0 ;
2574 uInt grpnum = 0 ;
2575 for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
2576 // Assumed that ifgfreq was sorted
2577 if ( grpnum != 0 && freqrange == ifgfreq[2*i] ) {
2578 freqgrp[grpnum-1].push_back( i ) ;
2579 }
2580 else {
2581 vector<uInt> grp0( 1, i ) ;
2582 freqgrp.push_back( grp0 ) ;
2583 grpnum++ ;
2584 }
2585 freqrange = ifgfreq[2*i+1] ;
2586 }
2587
2588
2589 // print IF groups
2590 cout << "IF Group summary: " << endl ;
2591 cout << " GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
2592 for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
2593 cout << " GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
2594 for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
2595 cout << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ;
2596 }
2597 cout << endl ;
2598 }
2599 cout << endl ;
2600
2601 // print frequency group
2602 cout << "Frequency Group summary: " << endl ;
2603 cout << " GROUP_ID [FREQ_MIN, FREQ_MAX]: IF_GROUP_ID" << endl ;
2604 for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
2605 cout << " GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*freqgrp[i][0]] << "," << ifgfreq[2*freqgrp[i][freqgrp[i].size()-1]+1] << "]: " ;
2606 for ( uInt j = 0 ; j < freqgrp[i].size() ; j++ ) {
2607 cout << freqgrp[i][j] << " " ;
2608 }
2609 cout << endl ;
2610 }
2611 cout << endl ;
2612
2613 // membership check
2614 // groups: list of IF group indexes whose frequency range overlaps with
2615 // that of each table and IF
2616 // groups[numtable][numIF][nummembership]
2617 // groups: [[[grp, grp,...], [grp, grp,...],...],
2618 // [[grp, grp,...], [grp, grp,...],...],
2619 // ...
2620 // [[grp, grp,...], [grp, grp,...],...]]
2621 vector< vector< vector<uInt> > > groups( insize ) ;
2622 for ( uInt i = 0 ; i < insize ; i++ ) {
2623 groups[i].resize( freqid[i].size() ) ;
2624 }
2625 for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
2626 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
2627 uInt tableid = ifgrp[igrp][2*imem] ;
2628 vector<uInt>::iterator iter = find( freqid[tableid].begin(), freqid[tableid].end(), ifgrp[igrp][2*imem+1] ) ;
2629 if ( iter != freqid[tableid].end() ) {
2630 uInt rowid = distance( freqid[tableid].begin(), iter ) ;
2631 groups[tableid][rowid].push_back( igrp ) ;
2632 }
2633 }
2634 }
2635
2636 // print membership
2637 //for ( uInt i = 0 ; i < insize ; i++ ) {
2638 //cout << "Table " << i << endl ;
2639 //for ( uInt j = 0 ; j < groups[i].size() ; j++ ) {
2640 //cout << " FREQ_ID " << setw( 2 ) << freqid[i][j] << ": " ;
2641 //for ( uInt k = 0 ; k < groups[i][j].size() ; k++ ) {
2642 //cout << setw( 2 ) << groups[i][j][k] << " " ;
2643 //}
2644 //cout << endl ;
2645 //}
2646 //}
2647
2648 // set back coordinfo
2649 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2650 vector<string> coordinfo = tmpin[itable]->getCoordInfo() ;
2651 coordinfo[0] = oldinfo[itable] ;
2652 tmpin[itable]->setCoordInfo( coordinfo ) ;
2653 }
2654
2655 // Create additional table if needed
2656 bool oldInsitu = insitu_ ;
2657 setInsitu( false ) ;
2658 vector< vector<uInt> > addrow( insize ) ;
2659 vector<uInt> addtable( insize, 0 ) ;
2660 vector<uInt> newtableids( insize ) ;
2661 vector<uInt> newifids( insize, 0 ) ;
2662 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2663 //cout << "Table " << setw(2) << itable << ": " ;
2664 for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
2665 addrow[itable].push_back( groups[itable][ifrow].size()-1 ) ;
2666 //cout << addrow[itable][ifrow] << " " ;
2667 }
2668 addtable[itable] = *max_element( addrow[itable].begin(), addrow[itable].end() ) ;
2669 //cout << "(" << addtable[itable] << ")" << endl ;
2670 }
2671 newin.resize( insize ) ;
2672 copy( tmpin.begin(), tmpin.end(), newin.begin() ) ;
2673 for ( uInt i = 0 ; i < insize ; i++ ) {
2674 newtableids[i] = i ;
2675 }
2676 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2677 for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
2678 CountedPtr<Scantable> add = getScantable( newin[itable], false ) ;
2679 vector<int> freqidlist ;
2680 for ( uInt i = 0 ; i < groups[itable].size() ; i++ ) {
2681 if ( groups[itable][i].size() > iadd + 1 ) {
2682 freqidlist.push_back( freqid[itable][i] ) ;
2683 }
2684 }
2685 stringstream taqlstream ;
2686 taqlstream << "SELECT FROM $1 WHERE FREQ_ID IN [" ;
2687 for ( uInt i = 0 ; i < freqidlist.size() ; i++ ) {
2688 taqlstream << i ;
2689 if ( i < freqidlist.size() - 1 )
2690 taqlstream << "," ;
2691 else
2692 taqlstream << "]" ;
2693 }
2694 string taql = taqlstream.str() ;
2695 //cout << "taql = " << taql << endl ;
2696 STSelector selector = STSelector() ;
2697 selector.setTaQL( taql ) ;
2698 add->setSelection( selector ) ;
2699 newin.push_back( add ) ;
2700 newtableids.push_back( itable ) ;
2701 newifids.push_back( iadd + 1 ) ;
2702 }
2703 }
2704
2705 // udpate ifgrp
2706 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2707 for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
2708 for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
2709 if ( groups[itable][ifrow].size() > iadd + 1 ) {
2710 uInt igrp = groups[itable][ifrow][iadd+1] ;
2711 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
2712 if ( ifgrp[igrp][2*imem] == newtableids[iadd+insize] && ifgrp[igrp][2*imem+1] == freqid[newtableids[iadd+insize]][ifrow] ) {
2713 ifgrp[igrp][2*imem] = insize + iadd ;
2714 }
2715 }
2716 }
2717 }
2718 }
2719 }
2720
2721 // print IF groups again for debug
2722 //cout << "IF Group summary: " << endl ;
2723 //cout << " GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
2724 //for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
2725 //cout << " GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
2726 //for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
2727 //cout << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ;
2728 //}
2729 //cout << endl ;
2730 //}
2731 //cout << endl ;
2732
2733 // reset SCANNO and IFNO/FREQ_ID: IF is reset by the result of sortation
2734 cout << "All scan number is set to 0" << endl ;
2735 //cout << "All IF number is set to IF group index" << endl ;
2736 cout << endl ;
2737 insize = newin.size() ;
2738 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2739 uInt rows = newin[itable]->nrow() ;
2740 Table &tmpt = newin[itable]->table() ;
2741 freqIDCol.attach( tmpt, "FREQ_ID" ) ;
2742 scannoCol.attach( tmpt, "SCANNO" ) ;
2743 ifnoCol.attach( tmpt, "IFNO" ) ;
2744 for ( uInt irow=0 ; irow < rows ; irow++ ) {
2745 scannoCol.put( irow, 0 ) ;
2746 uInt freqID = freqIDCol( irow ) ;
2747 vector<uInt>::iterator iter = find( freqid[newtableids[itable]].begin(), freqid[newtableids[itable]].end(), freqID ) ;
2748 if ( iter != freqid[newtableids[itable]].end() ) {
2749 uInt index = distance( freqid[newtableids[itable]].begin(), iter ) ;
2750 ifnoCol.put( irow, groups[newtableids[itable]][index][newifids[itable]] ) ;
2751 }
2752 else {
2753 throw(AipsError("IF grouping was wrong in additional tables.")) ;
2754 }
2755 }
2756 }
2757 oldinfo.resize( insize ) ;
2758 setInsitu( oldInsitu ) ;
2759
2760 // temporarily set coordinfo
2761 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2762 vector<string> coordinfo = newin[itable]->getCoordInfo() ;
2763 oldinfo[itable] = coordinfo[0] ;
2764 coordinfo[0] = "Hz" ;
2765 newin[itable]->setCoordInfo( coordinfo ) ;
2766 }
2767
2768 // save column values in the vector
2769 vector< vector<uInt> > freqTableIdVec( insize ) ;
2770 vector< vector<uInt> > freqIdVec( insize ) ;
2771 vector< vector<uInt> > ifNoVec( insize ) ;
2772 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2773 ScalarColumn<uInt> freqIDs ;
2774 freqIDs.attach( newin[itable]->frequencies().table(), "ID" ) ;
2775 ifnoCol.attach( newin[itable]->table(), "IFNO" ) ;
2776 freqIDCol.attach( newin[itable]->table(), "FREQ_ID" ) ;
2777 for ( uInt irow = 0 ; irow < newin[itable]->frequencies().table().nrow() ; irow++ ) {
2778 freqTableIdVec[itable].push_back( freqIDs( irow ) ) ;
2779 }
2780 for ( uInt irow = 0 ; irow < newin[itable]->table().nrow() ; irow++ ) {
2781 freqIdVec[itable].push_back( freqIDCol( irow ) ) ;
2782 ifNoVec[itable].push_back( ifnoCol( irow ) ) ;
2783 }
2784 }
2785
2786 // reset spectra and flagtra: pick up common part of frequency coverage
2787 //cout << "Pick common frequency range and align resolution" << endl ;
2788 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2789 uInt rows = newin[itable]->nrow() ;
2790 int nminchan = -1 ;
2791 int nmaxchan = -1 ;
2792 vector<uInt> freqIdUpdate ;
2793 for ( uInt irow = 0 ; irow < rows ; irow++ ) {
2794 uInt ifno = ifNoVec[itable][irow] ; // IFNO is reset by group index
2795 double minfreq = ifgfreq[2*ifno] ;
2796 double maxfreq = ifgfreq[2*ifno+1] ;
2797 //cout << "frequency range: [" << minfreq << "," << maxfreq << "]" << endl ;
2798 vector<double> abcissa = newin[itable]->getAbcissa( irow ) ;
2799 int nchan = abcissa.size() ;
2800 double resol = abcissa[1] - abcissa[0] ;
2801 //cout << "abcissa range : [" << abcissa[0] << "," << abcissa[nchan-1] << "]" << endl ;
2802 if ( minfreq <= abcissa[0] )
2803 nminchan = 0 ;
2804 else {
2805 //double cfreq = ( minfreq - abcissa[0] ) / resol ;
2806 double cfreq = ( minfreq - abcissa[0] + 0.5 * resol ) / resol ;
2807 nminchan = int(cfreq) + ( ( cfreq - int(cfreq) <= 0.5 ) ? 0 : 1 ) ;
2808 }
2809 if ( maxfreq >= abcissa[abcissa.size()-1] )
2810 nmaxchan = abcissa.size() - 1 ;
2811 else {
2812 //double cfreq = ( abcissa[abcissa.size()-1] - maxfreq ) / resol ;
2813 double cfreq = ( abcissa[abcissa.size()-1] - maxfreq + 0.5 * resol ) / resol ;
2814 nmaxchan = abcissa.size() - 1 - int(cfreq) - ( ( cfreq - int(cfreq) >= 0.5 ) ? 1 : 0 ) ;
2815 }
2816 //cout << "channel range (" << irow << "): [" << nminchan << "," << nmaxchan << "]" << endl ;
2817 if ( nmaxchan > nminchan ) {
2818 newin[itable]->reshapeSpectrum( nminchan, nmaxchan, irow ) ;
2819 int newchan = nmaxchan - nminchan + 1 ;
2820 if ( count( freqIdUpdate.begin(), freqIdUpdate.end(), freqIdVec[itable][irow] ) == 0 && newchan < nchan )
2821 freqIdUpdate.push_back( freqIdVec[itable][irow] ) ;
2822 }
2823 else {
2824 throw(AipsError("Failed to pick up common part of frequency range.")) ;
2825 }
2826 }
2827 for ( uInt i = 0 ; i < freqIdUpdate.size() ; i++ ) {
2828 uInt freqId = freqIdUpdate[i] ;
2829 Double refpix ;
2830 Double refval ;
2831 Double increment ;
2832
2833 // update row
2834 newin[itable]->frequencies().getEntry( refpix, refval, increment, freqId ) ;
2835 refval = refval - ( refpix - nminchan ) * increment ;
2836 refpix = 0 ;
2837 newin[itable]->frequencies().setEntry( refpix, refval, increment, freqId ) ;
2838 }
2839 }
2840
2841
2842 // reset spectra and flagtra: align spectral resolution
2843 //cout << "Align spectral resolution" << endl ;
2844 // gmaxdnu: the coarsest frequency resolution in the frequency group
2845 // gmemid: member index that have a resolution equal to gmaxdnu
2846 // gmaxdnu[numfreqgrp]
2847 // gmaxdnu: [dnu0, dnu1, ...]
2848 // gmemid[numfreqgrp]
2849 // gmemid: [id0, id1, ...]
2850 vector<double> gmaxdnu( freqgrp.size(), 0.0 ) ;
2851 vector<uInt> gmemid( freqgrp.size(), 0 ) ;
2852 for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
2853 double maxdnu = 0.0 ; // maximum (coarsest) frequency resolution
2854 int minchan = INT_MAX ; // minimum channel number
2855 Double refpixref = -1 ; // reference of 'reference pixel'
2856 Double refvalref = -1 ; // reference of 'reference frequency'
2857 Double refinc = -1 ; // reference frequency resolution
2858 uInt refreqid ;
2859 uInt reftable = INT_MAX;
2860 // process only if group member > 1
2861 if ( ifgrp[igrp].size() > 2 ) {
2862 // find minchan and maxdnu in each group
2863 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
2864 uInt tableid = ifgrp[igrp][2*imem] ;
2865 uInt rowid = ifgrp[igrp][2*imem+1] ;
2866 vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
2867 if ( iter != freqIdVec[tableid].end() ) {
2868 uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
2869 vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
2870 int nchan = abcissa.size() ;
2871 double dnu = abcissa[1] - abcissa[0] ;
2872 //cout << "GROUP " << igrp << " (" << tableid << "," << rowid << "): nchan = " << nchan << " (minchan = " << minchan << ")" << endl ;
2873 if ( nchan < minchan ) {
2874 minchan = nchan ;
2875 maxdnu = dnu ;
2876 newin[tableid]->frequencies().getEntry( refpixref, refvalref, refinc, rowid ) ;
2877 refreqid = rowid ;
2878 reftable = tableid ;
2879 }
2880 }
2881 }
2882 // regrid spectra in each group
2883 cout << "GROUP " << igrp << endl ;
2884 cout << " Channel number is adjusted to " << minchan << endl ;
2885 cout << " Corresponding frequency resolution is " << maxdnu << "Hz" << endl ;
2886 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
2887 uInt tableid = ifgrp[igrp][2*imem] ;
2888 uInt rowid = ifgrp[igrp][2*imem+1] ;
2889 freqIDCol.attach( newin[tableid]->table(), "FREQ_ID" ) ;
2890 //cout << "tableid = " << tableid << " rowid = " << rowid << ": " << endl ;
2891 //cout << " regridChannel applied to " ;
2892 if ( tableid != reftable )
2893 refreqid = newin[tableid]->frequencies().addEntry( refpixref, refvalref, refinc ) ;
2894 for ( uInt irow = 0 ; irow < newin[tableid]->table().nrow() ; irow++ ) {
2895 uInt tfreqid = freqIdVec[tableid][irow] ;
2896 if ( tfreqid == rowid ) {
2897 //cout << irow << " " ;
2898 newin[tableid]->regridChannel( minchan, maxdnu, irow ) ;
2899 freqIDCol.put( irow, refreqid ) ;
2900 freqIdVec[tableid][irow] = refreqid ;
2901 }
2902 }
2903 //cout << endl ;
2904 }
2905 }
2906 else {
2907 uInt tableid = ifgrp[igrp][0] ;
2908 uInt rowid = ifgrp[igrp][1] ;
2909 vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
2910 if ( iter != freqIdVec[tableid].end() ) {
2911 uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
2912 vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
2913 minchan = abcissa.size() ;
2914 maxdnu = abcissa[1] - abcissa[0] ;
2915 }
2916 }
2917 for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
2918 if ( count( freqgrp[i].begin(), freqgrp[i].end(), igrp ) > 0 ) {
2919 if ( maxdnu > gmaxdnu[i] ) {
2920 gmaxdnu[i] = maxdnu ;
2921 gmemid[i] = igrp ;
2922 }
2923 break ;
2924 }
2925 }
2926 }
2927 cout << endl ;
2928
2929 // set back coordinfo
2930 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2931 vector<string> coordinfo = newin[itable]->getCoordInfo() ;
2932 coordinfo[0] = oldinfo[itable] ;
2933 newin[itable]->setCoordInfo( coordinfo ) ;
2934 }
2935
2936 // accumulate all rows into the first table
2937 // NOTE: assumed in.size() = 1
2938 vector< CountedPtr<Scantable> > tmp( 1 ) ;
2939 if ( newin.size() == 1 )
2940 tmp[0] = newin[0] ;
2941 else
2942 tmp[0] = merge( newin ) ;
2943
2944 //return tmp[0] ;
2945
2946 // average
2947 CountedPtr<Scantable> tmpout = average( tmp, mask, weight, avmode ) ;
2948
2949 //return tmpout ;
2950
2951 // combine frequency group
2952 cout << "Combine spectra based on frequency grouping" << endl ;
2953 cout << "IFNO is renumbered as frequency group ID (see above)" << endl ;
2954 vector<string> coordinfo = tmpout->getCoordInfo() ;
2955 oldinfo[0] = coordinfo[0] ;
2956 coordinfo[0] = "Hz" ;
2957 tmpout->setCoordInfo( coordinfo ) ;
2958 // create proformas of output table
2959 stringstream taqlstream ;
2960 taqlstream << "SELECT FROM $1 WHERE IFNO IN [" ;
2961 for ( uInt i = 0 ; i < gmemid.size() ; i++ ) {
2962 taqlstream << gmemid[i] ;
2963 if ( i < gmemid.size() - 1 )
2964 taqlstream << "," ;
2965 else
2966 taqlstream << "]" ;
2967 }
2968 string taql = taqlstream.str() ;
2969 //cout << "taql = " << taql << endl ;
2970 STSelector selector = STSelector() ;
2971 selector.setTaQL( taql ) ;
2972 oldInsitu = insitu_ ;
2973 setInsitu( false ) ;
2974 out = getScantable( tmpout, false ) ;
2975 setInsitu( oldInsitu ) ;
2976 out->setSelection( selector ) ;
2977 // regrid rows
2978 ifnoCol.attach( tmpout->table(), "IFNO" ) ;
2979 for ( uInt irow = 0 ; irow < tmpout->table().nrow() ; irow++ ) {
2980 uInt ifno = ifnoCol( irow ) ;
2981 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
2982 if ( count( freqgrp[igrp].begin(), freqgrp[igrp].end(), ifno ) > 0 ) {
2983 vector<double> abcissa = tmpout->getAbcissa( irow ) ;
2984 double bw = ( abcissa[1] - abcissa[0] ) * abcissa.size() ;
2985 int nchan = (int)( bw / gmaxdnu[igrp] ) ;
2986 tmpout->regridChannel( nchan, gmaxdnu[igrp], irow ) ;
2987 break ;
2988 }
2989 }
2990 }
2991 // combine spectra
2992 ArrayColumn<Float> specColOut ;
2993 specColOut.attach( out->table(), "SPECTRA" ) ;
2994 ArrayColumn<uChar> flagColOut ;
2995 flagColOut.attach( out->table(), "FLAGTRA" ) ;
2996 ScalarColumn<uInt> ifnoColOut ;
2997 ifnoColOut.attach( out->table(), "IFNO" ) ;
2998 ScalarColumn<uInt> polnoColOut ;
2999 polnoColOut.attach( out->table(), "POLNO" ) ;
3000 ScalarColumn<uInt> freqidColOut ;
3001 freqidColOut.attach( out->table(), "FREQ_ID" ) ;
3002 MDirection::ScalarColumn dirColOut ;
3003 dirColOut.attach( out->table(), "DIRECTION" ) ;
3004 Table &tab = tmpout->table() ;
3005 Block<String> cols(1);
3006 cols[0] = String("POLNO") ;
3007 TableIterator iter( tab, cols ) ;
3008 bool done = false ;
3009 vector< vector<uInt> > sizes( freqgrp.size() ) ;
3010 while( !iter.pastEnd() ) {
3011 vector< vector<Float> > specout( freqgrp.size() ) ;
3012 vector< vector<uChar> > flagout( freqgrp.size() ) ;
3013 ArrayColumn<Float> specCols ;
3014 specCols.attach( iter.table(), "SPECTRA" ) ;
3015 ArrayColumn<uChar> flagCols ;
3016 flagCols.attach( iter.table(), "FLAGTRA" ) ;
3017 ifnoCol.attach( iter.table(), "IFNO" ) ;
3018 ScalarColumn<uInt> polnos ;
3019 polnos.attach( iter.table(), "POLNO" ) ;
3020 MDirection::ScalarColumn dircol ;
3021 dircol.attach( iter.table(), "DIRECTION" ) ;
3022 uInt polno = polnos( 0 ) ;
3023 //cout << "POLNO iteration: " << polno << endl ;
3024// for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3025// sizes[igrp].resize( freqgrp[igrp].size() ) ;
3026// for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
3027// for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
3028// uInt ifno = ifnoCol( irow ) ;
3029// if ( ifno == freqgrp[igrp][imem] ) {
3030// Vector<Float> spec = specCols( irow ) ;
3031// Vector<uChar> flag = flagCols( irow ) ;
3032// vector<Float> svec ;
3033// spec.tovector( svec ) ;
3034// vector<uChar> fvec ;
3035// flag.tovector( fvec ) ;
3036// //cout << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << endl ;
3037// specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ;
3038// flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ;
3039// //cout << "specout[" << igrp << "].size() = " << specout[igrp].size() << endl ;
3040// sizes[igrp][imem] = spec.nelements() ;
3041// }
3042// }
3043// }
3044// for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
3045// uInt ifout = ifnoColOut( irow ) ;
3046// uInt polout = polnoColOut( irow ) ;
3047// if ( ifout == gmemid[igrp] && polout == polno ) {
3048// // set SPECTRA and FRAGTRA
3049// Vector<Float> newspec( specout[igrp] ) ;
3050// Vector<uChar> newflag( flagout[igrp] ) ;
3051// specColOut.put( irow, newspec ) ;
3052// flagColOut.put( irow, newflag ) ;
3053// // IFNO renumbering
3054// ifnoColOut.put( irow, igrp ) ;
3055// }
3056// }
3057// }
3058 // get a list of number of channels for each frequency group member
3059 if ( !done ) {
3060 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3061 sizes[igrp].resize( freqgrp[igrp].size() ) ;
3062 for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
3063 for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
3064 uInt ifno = ifnoCol( irow ) ;
3065 if ( ifno == freqgrp[igrp][imem] ) {
3066 Vector<Float> spec = specCols( irow ) ;
3067 sizes[igrp][imem] = spec.nelements() ;
3068 break ;
3069 }
3070 }
3071 }
3072 }
3073 done = true ;
3074 }
3075 // combine spectra
3076 for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
3077 uInt polout = polnoColOut( irow ) ;
3078 if ( polout == polno ) {
3079 uInt ifout = ifnoColOut( irow ) ;
3080 Vector<Double> direction = dirColOut(irow).getAngle(Unit(String("rad"))).getValue() ;
3081 uInt igrp ;
3082 for ( uInt jgrp = 0 ; jgrp < freqgrp.size() ; jgrp++ ) {
3083 if ( ifout == gmemid[jgrp] ) {
3084 igrp = jgrp ;
3085 break ;
3086 }
3087 }
3088 for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
3089 for ( uInt jrow = 0 ; jrow < iter.table().nrow() ; jrow++ ) {
3090 uInt ifno = ifnoCol( jrow ) ;
3091 Vector<Double> tdir = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
3092 //if ( ifno == freqgrp[igrp][imem] && allTrue( tdir == direction ) ) {
3093 if ( ifno == freqgrp[igrp][imem] && allNearAbs( tdir, direction, tol ) ) {
3094 Vector<Float> spec = specCols( jrow ) ;
3095 Vector<uChar> flag = flagCols( jrow ) ;
3096 vector<Float> svec ;
3097 spec.tovector( svec ) ;
3098 vector<uChar> fvec ;
3099 flag.tovector( fvec ) ;
3100 //cout << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << endl ;
3101 specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ;
3102 flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ;
3103 //cout << "specout[" << igrp << "].size() = " << specout[igrp].size() << endl ;
3104 }
3105 }
3106 }
3107 // set SPECTRA and FRAGTRA
3108 Vector<Float> newspec( specout[igrp] ) ;
3109 Vector<uChar> newflag( flagout[igrp] ) ;
3110 specColOut.put( irow, newspec ) ;
3111 flagColOut.put( irow, newflag ) ;
3112 // IFNO renumbering
3113 ifnoColOut.put( irow, igrp ) ;
3114 }
3115 }
3116 iter++ ;
3117 }
3118 // update FREQUENCIES subtable
3119 vector<bool> updated( freqgrp.size(), false ) ;
3120 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3121 uInt index = 0 ;
3122 uInt pixShift = 0 ;
3123 while ( freqgrp[igrp][index] != gmemid[igrp] ) {
3124 pixShift += sizes[igrp][index++] ;
3125 }
3126 for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
3127 if ( ifnoColOut( irow ) == gmemid[igrp] && !updated[igrp] ) {
3128 uInt freqidOut = freqidColOut( irow ) ;
3129 //cout << "freqgrp " << igrp << " freqidOut = " << freqidOut << endl ;
3130 double refpix ;
3131 double refval ;
3132 double increm ;
3133 out->frequencies().getEntry( refpix, refval, increm, freqidOut ) ;
3134 refpix += pixShift ;
3135 out->frequencies().setEntry( refpix, refval, increm, freqidOut ) ;
3136 updated[igrp] = true ;
3137 }
3138 }
3139 }
3140
3141 //out = tmpout ;
3142
3143 coordinfo = tmpout->getCoordInfo() ;
3144 coordinfo[0] = oldinfo[0] ;
3145 tmpout->setCoordInfo( coordinfo ) ;
3146 }
3147 else {
3148 // simple average
3149 out = average( in, mask, weight, avmode ) ;
3150 }
3151
3152 return out ;
3153}
Note: See TracBrowser for help on using the repository browser.