source: trunk/src/STMath.cpp@ 2918

Last change on this file since 2918 was 2918, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: test_tsdaverage etc.

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Replaced STIdxIterExAcc with STIdxIter2.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 155.9 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 <sstream>
14
15#include <casa/iomanip.h>
16#include <casa/Arrays/MaskArrLogi.h>
17#include <casa/Arrays/MaskArrMath.h>
18#include <casa/Arrays/ArrayLogical.h>
19#include <casa/Arrays/ArrayMath.h>
20#include <casa/Arrays/Slice.h>
21#include <casa/Arrays/Slicer.h>
22#include <casa/BasicSL/String.h>
23#include <casa/Containers/Block.h>
24#include <casa/Containers/RecordField.h>
25#include <casa/Exceptions/Error.h>
26#include <casa/Logging/LogIO.h>
27#include <casa/Quanta/Quantum.h>
28
29#include <coordinates/Coordinates/CoordinateSystem.h>
30#include <coordinates/Coordinates/CoordinateUtil.h>
31#include <coordinates/Coordinates/FrequencyAligner.h>
32#include <coordinates/Coordinates/SpectralCoordinate.h>
33
34#include <lattices/Lattices/LatticeUtilities.h>
35
36#include <scimath/Functionals/Polynomial.h>
37#include <scimath/Mathematics/Convolver.h>
38#include <scimath/Mathematics/VectorKernel.h>
39
40#include <tables/Tables/ExprNode.h>
41#include <tables/Tables/ReadAsciiTable.h>
42#include <tables/Tables/TableCopy.h>
43#include <tables/Tables/TableIter.h>
44#include <tables/Tables/TableParse.h>
45#include <tables/Tables/TableRecord.h>
46#include <tables/Tables/TableRow.h>
47#include <tables/Tables/TableVector.h>
48#include <tables/Tables/TabVecMath.h>
49
50#include <atnf/PKSIO/SrcType.h>
51
52#include "RowAccumulator.h"
53#include "STAttr.h"
54#include "STMath.h"
55#include "STSelector.h"
56#include "Accelerator.h"
57#include "STIdxIter.h"
58
59#include "CalibrationHelper.h"
60
61using namespace casa;
62using namespace asap;
63
64// namespace {
65// class CalibrationIterator
66// {
67// public:
68// CalibrationIterator()
69// : table_(table)
70// {}
71// ~CalibrationIterator() {}
72// void Iterate() {
73// STSelector sel;
74// vector<string> cols( 3 ) ;
75// cols[0] = "BEAMNO" ;
76// cols[1] = "POLNO" ;
77// cols[2] = "IFNO" ;
78// STIdxIter2 iter(out, cols);
79// while ( !iter.pastEnd() ) {
80// Record ids = iter.currentValue();
81// stringstream ss ;
82// ss << "SELECT FROM $1 WHERE "
83// << "BEAMNO==" << ids.asuInt("BEAMNO") << "&&"
84// << "POLNO==" << ids.asuInt("POLNO") << "&&"
85// << "IFNO==" << ids.asuInt("IFNO") ;
86// //cout << "TaQL string: " << ss.str() << endl ;
87// sel.setTaQL( ss.str() ) ;
88// ////
89// // begin individual processing
90// aoff->setSelection( sel ) ;
91// ahot->setSelection( sel ) ;
92// asky->setSelection( sel ) ;
93// Vector<uInt> rows = iter.getRows( SHARE ) ;
94// // out should be an exact copy of s except that SPECTRA column is empty
95// calibrateCW( out, s, aoff, asky, ahot, acold, rows, antname ) ;
96// aoff->unsetSelection() ;
97// ahot->unsetSelection() ;
98// asky->unsetSelection() ;
99// // end individual processing
100// ////
101// sel.reset() ;
102// iter.next() ;
103// }
104// };
105// protected:
106// CountedPtr<Scantable> table_;
107// };
108
109
110// } // anonymous namespace
111
112// 2012/02/17 TN
113// Since STGrid is implemented, average doesn't consider direction
114// when accumulating
115// tolerance for direction comparison (rad)
116// #define TOL_OTF 1.0e-15
117// #define TOL_POINT 2.9088821e-4 // 1 arcmin
118
119STMath::STMath(bool insitu) :
120 insitu_(insitu)
121{
122}
123
124
125STMath::~STMath()
126{
127}
128
129CountedPtr<Scantable>
130STMath::average( const std::vector<CountedPtr<Scantable> >& in,
131 const std::vector<bool>& mask,
132 const std::string& weight,
133 const std::string& avmode)
134{
135// double t0, t1 ;
136// t0 = mathutil::gettimeofday_sec() ;
137
138 LogIO os( LogOrigin( "STMath", "average()", WHERE ) ) ;
139 if ( avmode == "SCAN" && in.size() != 1 )
140 throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
141 "Use merge first."));
142 WeightType wtype = stringToWeight(weight);
143
144 // 2012/02/17 TN
145 // Since STGrid is implemented, average doesn't consider direction
146 // when accumulating
147 // check if OTF observation
148// String obstype = in[0]->getHeader().obstype ;
149// Double tol = 0.0 ;
150// if ( (obstype.find( "OTF" ) != String::npos) || (obstype.find( "OBSERVE_TARGET" ) != String::npos) ) {
151// tol = TOL_OTF ;
152// }
153// else {
154// tol = TOL_POINT ;
155// }
156
157 // output
158 // clone as this is non insitu
159 bool insitu = insitu_;
160 setInsitu(false);
161 CountedPtr< Scantable > out = getScantable(in[0], true);
162 setInsitu(insitu);
163 std::vector<CountedPtr<Scantable> >::const_iterator stit = in.begin();
164 ++stit;
165 while ( stit != in.end() ) {
166 out->appendToHistoryTable((*stit)->history());
167 ++stit;
168 }
169
170 Table& tout = out->table();
171
172 /// @todo check if all scantables are conformant
173
174 ArrayColumn<Float> specColOut(tout,"SPECTRA");
175 ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
176 ArrayColumn<Float> tsysColOut(tout,"TSYS");
177 ScalarColumn<Double> mjdColOut(tout,"TIME");
178 ScalarColumn<Double> intColOut(tout,"INTERVAL");
179 ScalarColumn<uInt> cycColOut(tout,"CYCLENO");
180 ScalarColumn<uInt> scanColOut(tout,"SCANNO");
181 ScalarColumn<uInt> flagRowColOut(tout,"FLAGROW");
182
183 // set up the output table rows. These are based on the structure of the
184 // FIRST scantable in the vector
185 const Table& baset = in[0]->table();
186
187 RowAccumulator acc(wtype);
188 Vector<Bool> cmask(mask);
189 acc.setUserMask(cmask);
190// ROTableRow row(tout);
191 ROArrayColumn<Float> specCol, tsysCol;
192 ROArrayColumn<uChar> flagCol;
193 ROScalarColumn<Double> mjdCol, intCol;
194 ROScalarColumn<Int> scanIDCol;
195 ROScalarColumn<uInt> flagRowCol;
196
197 //Vector<uInt> rowstodelete;
198 Block<uInt> rowstodelB( in[0]->nrow() ) ;
199 uInt nrowdel = 0 ;
200
201// Block<String> cols(3);
202 vector<string> cols(3) ;
203 cols[0] = String("BEAMNO");
204 cols[1] = String("IFNO");
205 cols[2] = String("POLNO");
206 if ( avmode == "SOURCE" ) {
207 cols.resize(4);
208 cols[3] = String("SRCNAME");
209 }
210 if ( avmode == "SCAN" && in.size() == 1) {
211 //cols.resize(4);
212 //cols[3] = String("SCANNO");
213 cols.resize(5);
214 cols[3] = String("SRCNAME");
215 cols[4] = String("SCANNO");
216 }
217 uInt outrowCount = 0;
218 // use STIdxIter2 instead of TableIterator
219 STIdxIter2 iter( in[0], cols ) ;
220// double t2 = 0 ;
221// double t3 = 0 ;
222// double t4 = 0 ;
223// double t5 = 0 ;
224// TableIterator iter(baset, cols);
225// int count = 0 ;
226 while (!iter.pastEnd()) {
227 Vector<uInt> rows = iter.getRows( SHARE ) ;
228 if ( rows.nelements() == 0 ) {
229 iter.next() ;
230 continue ;
231 }
232 Record current = iter.currentValue() ;
233 //Table subt = iter.table();
234 // copy the first row of this selection into the new table
235 tout.addRow();
236// t4 = mathutil::gettimeofday_sec() ;
237 // skip to copy SPECTRA, FLAGTRA, and TSYS since those heavy columns are
238 // overwritten in the following process
239 copyRows( tout, baset, outrowCount, rows[0], 1, False, False, False ) ;
240// t5 += mathutil::gettimeofday_sec() - t4 ;
241 // re-index to 0
242 if ( avmode != "SCAN" && avmode != "SOURCE" ) {
243 scanColOut.put(outrowCount, uInt(0));
244 }
245
246 // 2012/02/17 TN
247 // Since STGrid is implemented, average doesn't consider direction
248 // when accumulating
249// MDirection::ScalarColumn dircol ;
250// dircol.attach( subt, "DIRECTION" ) ;
251// Int length = subt.nrow() ;
252// vector< Vector<Double> > dirs ;
253// vector<int> indexes ;
254// for ( Int i = 0 ; i < length ; i++ ) {
255// Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
256// //os << << count++ << ": " ;
257// //os << "[" << t[0] << "," << t[1] << "]" << LogIO::POST ;
258// bool adddir = true ;
259// for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
260// //if ( allTrue( t == dirs[j] ) ) {
261// Double dx = t[0] - dirs[j][0] ;
262// Double dy = t[1] - dirs[j][1] ;
263// Double dd = sqrt( dx * dx + dy * dy ) ;
264// //if ( allNearAbs( t, dirs[j], tol ) ) {
265// if ( dd <= tol ) {
266// adddir = false ;
267// break ;
268// }
269// }
270// if ( adddir ) {
271// dirs.push_back( t ) ;
272// indexes.push_back( i ) ;
273// }
274// }
275// uInt rowNum = dirs.size() ;
276// tout.addRow( rowNum ) ;
277// for ( uInt i = 0 ; i < rowNum ; i++ ) {
278// TableCopy::copyRows( tout, subt, outrowCount+i, indexes[i], 1 ) ;
279// // re-index to 0
280// if ( avmode != "SCAN" && avmode != "SOURCE" ) {
281// scanColOut.put(outrowCount+i, uInt(0));
282// }
283// }
284// outrowCount += rowNum ;
285
286 // merge loop
287 uInt i = outrowCount ;
288 // in[0] is already selected by iterator
289 specCol.attach(baset,"SPECTRA");
290 flagCol.attach(baset,"FLAGTRA");
291 tsysCol.attach(baset,"TSYS");
292 intCol.attach(baset,"INTERVAL");
293 mjdCol.attach(baset,"TIME");
294 flagRowCol.attach(baset,"FLAGROW");
295 Vector<Float> spec,tsys;
296 Vector<uChar> flag;
297 Double inter,time;
298 uInt flagRow;
299
300 for (uInt l = 0; l < rows.nelements(); ++l ) {
301 uInt k = rows[l] ;
302 flagCol.get(k, flag);
303 Vector<Bool> bflag(flag.shape());
304 flagRowCol.get(k, flagRow);
305 if (flagRow > 0)
306 bflag = true;
307 else
308 convertArray(bflag, flag);
309 /*
310 if ( allEQ(bflag, True) ) {
311 continue;//don't accumulate
312 }
313 */
314 specCol.get(k, spec);
315 tsysCol.get(k, tsys);
316 intCol.get(k, inter);
317 mjdCol.get(k, time);
318 // spectrum has to be added last to enable weighting by the other values
319// t2 = mathutil::gettimeofday_sec() ;
320 acc.add(spec, !bflag, tsys, inter, time);
321// t3 += mathutil::gettimeofday_sec() - t2 ;
322
323 }
324
325
326 // in[0] is already selected by TableIterator so that index is
327 // started from 1
328 for ( int j=1; j < int(in.size()); ++j ) {
329 const Table& tin = in[j]->table();
330 //const TableRecord& rec = row.get(i);
331 ROScalarColumn<Double> tmp(tin, "TIME");
332 Double td;tmp.get(0,td);
333
334#if 1
335 static char const*const colNames1[] = { "IFNO", "BEAMNO", "POLNO" };
336 //uInt const values1[] = { rec.asuInt("IFNO"), rec.asuInt("BEAMNO"), rec.asuInt("POLNO") };
337 uInt const values1[] = { current.asuInt("IFNO"), current.asuInt("BEAMNO"), current.asuInt("POLNO") };
338 SingleTypeEqPredicate<uInt, 3> myPred(tin, colNames1, values1);
339 CustomTableExprNodeRep myNodeRep(tin, myPred);
340 myNodeRep.link(); // to avoid automatic delete when myExpr is destructed.
341 CustomTableExprNode myExpr(myNodeRep);
342 Table basesubt = tin(myExpr);
343#else
344// Table basesubt = tin( tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
345// && tin.col("IFNO") == Int(rec.asuInt("IFNO"))
346// && tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
347 Table basesubt = tin( tin.col("BEAMNO") == current.asuInt("BEAMNO")
348 && tin.col("IFNO") == current.asuInt("IFNO")
349 && tin.col("POLNO") == current.asuInt("POLNO") );
350#endif
351 Table subt;
352 if ( avmode == "SOURCE") {
353// subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME"));
354 subt = basesubt( basesubt.col("SRCNAME") == current.asString("SRCNAME") );
355
356 } else if (avmode == "SCAN") {
357// subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME")
358// && basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
359 subt = basesubt( basesubt.col("SRCNAME") == current.asString("SRCNAME")
360 && basesubt.col("SCANNO") == current.asuInt("SCANNO") );
361 } else {
362 subt = basesubt;
363 }
364
365 // 2012/02/17 TN
366 // Since STGrid is implemented, average doesn't consider direction
367 // when accumulating
368// vector<uInt> removeRows ;
369// uInt nrsubt = subt.nrow() ;
370// for ( uInt irow = 0 ; irow < nrsubt ; irow++ ) {
371// //if ( !allTrue((subt.col("DIRECTION").getArrayDouble(TableExprId(irow)))==rec.asArrayDouble("DIRECTION")) ) {
372// Vector<Double> x0 = (subt.col("DIRECTION").getArrayDouble(TableExprId(irow))) ;
373// Vector<Double> x1 = rec.asArrayDouble("DIRECTION") ;
374// double dx = x0[0] - x1[0];
375// double dy = x0[1] - x1[1];
376// Double dd = sqrt( dx * dx + dy * dy ) ;
377// //if ( !allNearAbs((subt.col("DIRECTION").getArrayDouble(TableExprId(irow))), rec.asArrayDouble("DIRECTION"), tol ) ) {
378// if ( dd > tol ) {
379// removeRows.push_back( irow ) ;
380// }
381// }
382// if ( removeRows.size() != 0 ) {
383// subt.removeRow( removeRows ) ;
384// }
385
386// if ( nrsubt == removeRows.size() )
387// throw(AipsError("Averaging data is empty.")) ;
388
389 specCol.attach(subt,"SPECTRA");
390 flagCol.attach(subt,"FLAGTRA");
391 tsysCol.attach(subt,"TSYS");
392 intCol.attach(subt,"INTERVAL");
393 mjdCol.attach(subt,"TIME");
394 flagRowCol.attach(subt,"FLAGROW");
395 for (uInt k = 0; k < subt.nrow(); ++k ) {
396 flagCol.get(k, flag);
397 Vector<Bool> bflag(flag.shape());
398 flagRowCol.get(k, flagRow);
399 if (flagRow > 0)
400 bflag = true;
401 else
402 convertArray(bflag, flag);
403 /*
404 if ( allEQ(bflag, True) ) {
405 continue;//don't accumulate
406 }
407 */
408 specCol.get(k, spec);
409 //tsysCol.get(k, tsys);
410 tsys.assign( tsysCol(k) );
411 intCol.get(k, inter);
412 mjdCol.get(k, time);
413 // spectrum has to be added last to enable weighting by the other values
414// t2 = mathutil::gettimeofday_sec() ;
415 acc.add(spec, !bflag, tsys, inter, time);
416// t3 += mathutil::gettimeofday_sec() - t2 ;
417 }
418
419 }
420 const Vector<Bool>& msk = acc.getMask();
421 if ( allEQ(msk, False) ) {
422 rowstodelB[nrowdel] = i ;
423 nrowdel++ ;
424 continue;
425 }
426 //write out
427 if (acc.state()) {
428 // If there exists a channel at which all the input spectra are masked,
429 // spec has 'nan' values for that channel and it may affect the following
430 // processes. To avoid this, replacing 'nan' values in spec with
431 // weighted-mean of all spectra in the following line.
432 // (done for CAS-2776, 2011/04/07 by Wataru Kawasaki)
433 acc.replaceNaN();
434
435 Vector<uChar> flg(msk.shape());
436 convertArray(flg, !msk);
437 for (uInt k = 0; k < flg.nelements(); ++k) {
438 uChar userFlag = 1 << 7;
439 if (msk[k]==True) userFlag = 0 << 7;
440 flg(k) = userFlag;
441 }
442
443 flagColOut.put(i, flg);
444 specColOut.put(i, acc.getSpectrum());
445 tsysColOut.put(i, acc.getTsys());
446 intColOut.put(i, acc.getInterval());
447 mjdColOut.put(i, acc.getTime());
448 // we should only have one cycle now -> reset it to be 0
449 // frequency switched data has different CYCLENO for different IFNO
450 // which requires resetting this value
451 cycColOut.put(i, uInt(0));
452 // completely flagged rows are removed anyway
453 flagRowColOut.put(i, uInt(0));
454 } else {
455 os << "For output row="<<i<<", all input rows of data are flagged. no averaging" << LogIO::POST;
456 }
457 acc.reset();
458
459 // merge with while loop for preparing out table
460 ++outrowCount;
461// ++iter ;
462 iter.next() ;
463 }
464
465 if ( nrowdel > 0 ) {
466 Vector<uInt> rowstodelete( IPosition(1,nrowdel), rowstodelB.storage(), SHARE ) ;
467 os << rowstodelete << LogIO::POST ;
468 tout.removeRow(rowstodelete);
469 if (tout.nrow() == 0) {
470 throw(AipsError("Can't average fully flagged data."));
471 }
472 }
473
474// t1 = mathutil::gettimeofday_sec() ;
475// cout << "elapsed time for average(): " << t1-t0 << " sec" << endl ;
476// cout << " elapsed time for acc.add(): " << t3 << " sec" << endl ;
477// cout << " elapsed time for copyRows(): " << t5 << " sec" << endl ;
478
479 return out;
480}
481
482CountedPtr< Scantable >
483STMath::averageChannel( const CountedPtr < Scantable > & in,
484 const std::string & mode,
485 const std::string& avmode )
486{
487 (void) mode; // currently unused
488 // 2012/02/17 TN
489 // Since STGrid is implemented, average doesn't consider direction
490 // when accumulating
491 // check if OTF observation
492// String obstype = in->getHeader().obstype ;
493// Double tol = 0.0 ;
494// if ( obstype.find( "OTF" ) != String::npos ) {
495// tol = TOL_OTF ;
496// }
497// else {
498// tol = TOL_POINT ;
499// }
500
501 // clone as this is non insitu
502 bool insitu = insitu_;
503 setInsitu(false);
504 CountedPtr< Scantable > out = getScantable(in, true);
505 setInsitu(insitu);
506 Table& tout = out->table();
507 ArrayColumn<Float> specColOut(tout,"SPECTRA");
508 ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
509 ArrayColumn<Float> tsysColOut(tout,"TSYS");
510 ScalarColumn<uInt> scanColOut(tout,"SCANNO");
511 ScalarColumn<Double> intColOut(tout, "INTERVAL");
512 Table tmp = in->table().sort("BEAMNO");
513 Block<String> cols(3);
514 cols[0] = String("BEAMNO");
515 cols[1] = String("IFNO");
516 cols[2] = String("POLNO");
517 if ( avmode == "SCAN") {
518 cols.resize(4);
519 cols[3] = String("SCANNO");
520 }
521 uInt outrowCount = 0;
522 uChar userflag = 1 << 7;
523 TableIterator iter(tmp, cols);
524 while (!iter.pastEnd()) {
525 Table subt = iter.table();
526 ROArrayColumn<Float> specCol, tsysCol;
527 ROArrayColumn<uChar> flagCol;
528 ROScalarColumn<Double> intCol(subt, "INTERVAL");
529 specCol.attach(subt,"SPECTRA");
530 flagCol.attach(subt,"FLAGTRA");
531 tsysCol.attach(subt,"TSYS");
532
533 tout.addRow();
534 TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
535 if ( avmode != "SCAN") {
536 scanColOut.put(outrowCount, uInt(0));
537 }
538 Vector<Float> tmp;
539 specCol.get(0, tmp);
540 uInt nchan = tmp.nelements();
541 // have to do channel by channel here as MaskedArrMath
542 // doesn't have partialMedians
543 Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
544 Vector<Float> outspec(nchan);
545 Vector<uChar> outflag(nchan,0);
546 Vector<Float> outtsys(1);/// @fixme when tsys is channel based
547 for (uInt i=0; i<nchan; ++i) {
548 Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
549 MaskedArray<Float> ma = maskedArray(specs,flags);
550 outspec[i] = median(ma);
551 if ( allEQ(ma.getMask(), False) )
552 outflag[i] = userflag;// flag data
553 }
554 outtsys[0] = median(tsysCol.getColumn());
555 specColOut.put(outrowCount, outspec);
556 flagColOut.put(outrowCount, outflag);
557 tsysColOut.put(outrowCount, outtsys);
558 Double intsum = sum(intCol.getColumn());
559 intColOut.put(outrowCount, intsum);
560 ++outrowCount;
561 ++iter;
562
563 // 2012/02/17 TN
564 // Since STGrid is implemented, average doesn't consider direction
565 // when accumulating
566// MDirection::ScalarColumn dircol ;
567// dircol.attach( subt, "DIRECTION" ) ;
568// Int length = subt.nrow() ;
569// vector< Vector<Double> > dirs ;
570// vector<int> indexes ;
571// // Handle MX mode averaging
572// if (in->nbeam() > 1 ) {
573// length = 1;
574// }
575// for ( Int i = 0 ; i < length ; i++ ) {
576// Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
577// bool adddir = true ;
578// for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
579// //if ( allTrue( t == dirs[j] ) ) {
580// Double dx = t[0] - dirs[j][0] ;
581// Double dy = t[1] - dirs[j][1] ;
582// Double dd = sqrt( dx * dx + dy * dy ) ;
583// //if ( allNearAbs( t, dirs[j], tol ) ) {
584// if ( dd <= tol ) {
585// adddir = false ;
586// break ;
587// }
588// }
589// if ( adddir ) {
590// dirs.push_back( t ) ;
591// indexes.push_back( i ) ;
592// }
593// }
594// uInt rowNum = dirs.size() ;
595// tout.addRow( rowNum );
596// for ( uInt i = 0 ; i < rowNum ; i++ ) {
597// TableCopy::copyRows(tout, subt, outrowCount+i, indexes[i], 1) ;
598// // Handle MX mode averaging
599// if ( avmode != "SCAN") {
600// scanColOut.put(outrowCount+i, uInt(0));
601// }
602// }
603// MDirection::ScalarColumn dircolOut ;
604// dircolOut.attach( tout, "DIRECTION" ) ;
605// for ( uInt irow = 0 ; irow < rowNum ; irow++ ) {
606// Vector<Double> t = \
607// dircolOut(outrowCount+irow).getAngle(Unit(String("rad"))).getValue() ;
608// Vector<Float> tmp;
609// specCol.get(0, tmp);
610// uInt nchan = tmp.nelements();
611// // have to do channel by channel here as MaskedArrMath
612// // doesn't have partialMedians
613// Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
614// // mask spectra for different DIRECTION
615// for ( uInt jrow = 0 ; jrow < subt.nrow() ; jrow++ ) {
616// Vector<Double> direction = \
617// dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
618// //if ( t[0] != direction[0] || t[1] != direction[1] ) {
619// Double dx = t[0] - direction[0];
620// Double dy = t[1] - direction[1];
621// Double dd = sqrt(dx*dx + dy*dy);
622// //if ( !allNearAbs( t, direction, tol ) ) {
623// if ( dd > tol && in->nbeam() < 2 ) {
624// flags[jrow] = userflag ;
625// }
626// }
627// Vector<Float> outspec(nchan);
628// Vector<uChar> outflag(nchan,0);
629// Vector<Float> outtsys(1);/// @fixme when tsys is channel based
630// for (uInt i=0; i<nchan; ++i) {
631// Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
632// MaskedArray<Float> ma = maskedArray(specs,flags);
633// outspec[i] = median(ma);
634// if ( allEQ(ma.getMask(), False) )
635// outflag[i] = userflag;// flag data
636// }
637// outtsys[0] = median(tsysCol.getColumn());
638// specColOut.put(outrowCount+irow, outspec);
639// flagColOut.put(outrowCount+irow, outflag);
640// tsysColOut.put(outrowCount+irow, outtsys);
641// Vector<Double> integ = intCol.getColumn() ;
642// MaskedArray<Double> mi = maskedArray( integ, flags ) ;
643// Double intsum = sum(mi);
644// intColOut.put(outrowCount+irow, intsum);
645// }
646// outrowCount += rowNum ;
647// ++iter;
648 }
649 return out;
650}
651
652CountedPtr< Scantable > STMath::getScantable(const CountedPtr< Scantable >& in,
653 bool droprows)
654{
655 if (insitu_) {
656 return in;
657 }
658 else {
659 // clone
660 return CountedPtr<Scantable>(new Scantable(*in, Bool(droprows)));
661 }
662}
663
664CountedPtr< Scantable > STMath::unaryOperate( const CountedPtr< Scantable >& in,
665 float val,
666 const std::string& mode,
667 bool tsys )
668{
669 CountedPtr< Scantable > out = getScantable(in, false);
670 Table& tab = out->table();
671 ArrayColumn<Float> specCol(tab,"SPECTRA");
672 ArrayColumn<Float> tsysCol(tab,"TSYS");
673 if (mode=="DIV") val = 1.0/val ;
674 else if (mode=="SUB") val *= -1.0 ;
675 for (uInt i=0; i<tab.nrow(); ++i) {
676 Vector<Float> spec;
677 Vector<Float> ts;
678 specCol.get(i, spec);
679 tsysCol.get(i, ts);
680 if (mode == "MUL" || mode == "DIV") {
681 //if (mode == "DIV") val = 1.0/val;
682 spec *= val;
683 specCol.put(i, spec);
684 if ( tsys ) {
685 ts *= val;
686 tsysCol.put(i, ts);
687 }
688 } else if ( mode == "ADD" || mode == "SUB") {
689 //if (mode == "SUB") val *= -1.0;
690 spec += val;
691 specCol.put(i, spec);
692 if ( tsys ) {
693 ts += val;
694 tsysCol.put(i, ts);
695 }
696 }
697 }
698 return out;
699}
700
701CountedPtr< Scantable > STMath::arrayOperate( const CountedPtr< Scantable >& in,
702 const std::vector<float> val,
703 const std::string& mode,
704 const std::string& opmode,
705 bool tsys )
706{
707 CountedPtr< Scantable > out ;
708 if ( opmode == "channel" ) {
709 out = arrayOperateChannel( in, val, mode, tsys ) ;
710 }
711 else if ( opmode == "row" ) {
712 out = arrayOperateRow( in, val, mode, tsys ) ;
713 }
714 else {
715 throw( AipsError( "Unknown array operation mode." ) ) ;
716 }
717 return out ;
718}
719
720CountedPtr< Scantable > STMath::arrayOperateChannel( const CountedPtr< Scantable >& in,
721 const std::vector<float> val,
722 const std::string& mode,
723 bool tsys )
724{
725 if ( val.size() == 1 ){
726 return unaryOperate( in, val[0], mode, tsys ) ;
727 }
728
729 // conformity of SPECTRA and TSYS
730 if ( tsys ) {
731 TableIterator titer(in->table(), "IFNO");
732 while ( !titer.pastEnd() ) {
733 ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ;
734 ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ;
735 Array<Float> spec = specCol.getColumn() ;
736 Array<Float> ts = tsysCol.getColumn() ;
737 if ( !spec.conform( ts ) ) {
738 throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ;
739 }
740 titer.next() ;
741 }
742 }
743
744 // check if all spectra in the scantable have the same number of channel
745 vector<uInt> nchans;
746 vector<uInt> ifnos = in->getIFNos() ;
747 for ( uInt i = 0 ; i < ifnos.size() ; i++ ) {
748 nchans.push_back( in->nchan( ifnos[i] ) ) ;
749 }
750 Vector<uInt> mchans( nchans ) ;
751 if ( anyNE( mchans, mchans[0] ) ) {
752 throw( AipsError("All spectra in the input scantable must have the same number of channel for vector operation." ) ) ;
753 }
754
755 // check if vector size is equal to nchan
756 Vector<Float> fact( val ) ;
757 if ( fact.nelements() != mchans[0] ) {
758 throw( AipsError("Vector size must be 1 or be same as number of channel.") ) ;
759 }
760
761 // check divided by zero
762 if ( ( mode == "DIV" ) && anyEQ( fact, (float)0.0 ) ) {
763 throw( AipsError("Divided by zero is not recommended." ) ) ;
764 }
765
766 CountedPtr< Scantable > out = getScantable(in, false);
767 Table& tab = out->table();
768 ArrayColumn<Float> specCol(tab,"SPECTRA");
769 ArrayColumn<Float> tsysCol(tab,"TSYS");
770 if (mode == "DIV") fact = (float)1.0 / fact;
771 else if (mode == "SUB") fact *= (float)-1.0 ;
772 for (uInt i=0; i<tab.nrow(); ++i) {
773 Vector<Float> spec;
774 Vector<Float> ts;
775 specCol.get(i, spec);
776 tsysCol.get(i, ts);
777 if (mode == "MUL" || mode == "DIV") {
778 //if (mode == "DIV") fact = (float)1.0 / fact;
779 spec *= fact;
780 specCol.put(i, spec);
781 if ( tsys ) {
782 ts *= fact;
783 tsysCol.put(i, ts);
784 }
785 } else if ( mode == "ADD" || mode == "SUB") {
786 //if (mode == "SUB") fact *= (float)-1.0 ;
787 spec += fact;
788 specCol.put(i, spec);
789 if ( tsys ) {
790 ts += fact;
791 tsysCol.put(i, ts);
792 }
793 }
794 }
795 return out;
796}
797
798CountedPtr< Scantable > STMath::arrayOperateRow( const CountedPtr< Scantable >& in,
799 const std::vector<float> val,
800 const std::string& mode,
801 bool tsys )
802{
803 if ( val.size() == 1 ) {
804 return unaryOperate( in, val[0], mode, tsys ) ;
805 }
806
807 // conformity of SPECTRA and TSYS
808 if ( tsys ) {
809 TableIterator titer(in->table(), "IFNO");
810 while ( !titer.pastEnd() ) {
811 ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ;
812 ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ;
813 Array<Float> spec = specCol.getColumn() ;
814 Array<Float> ts = tsysCol.getColumn() ;
815 if ( !spec.conform( ts ) ) {
816 throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ;
817 }
818 titer.next() ;
819 }
820 }
821
822 // check if vector size is equal to nrow
823 Vector<Float> fact( val ) ;
824 if (fact.nelements() != uInt(in->nrow())) {
825 throw( AipsError("Vector size must be 1 or be same as number of row.") ) ;
826 }
827
828 // check divided by zero
829 if ( ( mode == "DIV" ) && anyEQ( fact, (float)0.0 ) ) {
830 throw( AipsError("Divided by zero is not recommended." ) ) ;
831 }
832
833 CountedPtr< Scantable > out = getScantable(in, false);
834 Table& tab = out->table();
835 ArrayColumn<Float> specCol(tab,"SPECTRA");
836 ArrayColumn<Float> tsysCol(tab,"TSYS");
837 if (mode == "DIV") fact = (float)1.0 / fact;
838 if (mode == "SUB") fact *= (float)-1.0 ;
839 for (uInt i=0; i<tab.nrow(); ++i) {
840 Vector<Float> spec;
841 Vector<Float> ts;
842 specCol.get(i, spec);
843 tsysCol.get(i, ts);
844 if (mode == "MUL" || mode == "DIV") {
845 spec *= fact[i];
846 specCol.put(i, spec);
847 if ( tsys ) {
848 ts *= fact[i];
849 tsysCol.put(i, ts);
850 }
851 } else if ( mode == "ADD" || mode == "SUB") {
852 spec += fact[i];
853 specCol.put(i, spec);
854 if ( tsys ) {
855 ts += fact[i];
856 tsysCol.put(i, ts);
857 }
858 }
859 }
860 return out;
861}
862
863CountedPtr< Scantable > STMath::array2dOperate( const CountedPtr< Scantable >& in,
864 const std::vector< std::vector<float> > val,
865 const std::string& mode,
866 bool tsys )
867{
868 // conformity of SPECTRA and TSYS
869 if ( tsys ) {
870 TableIterator titer(in->table(), "IFNO");
871 while ( !titer.pastEnd() ) {
872 ArrayColumn<Float> specCol( in->table(), "SPECTRA" ) ;
873 ArrayColumn<Float> tsysCol( in->table(), "TSYS" ) ;
874 Array<Float> spec = specCol.getColumn() ;
875 Array<Float> ts = tsysCol.getColumn() ;
876 if ( !spec.conform( ts ) ) {
877 throw( AipsError( "SPECTRA and TSYS must conform in shape if you want to apply operation on Tsys." ) ) ;
878 }
879 titer.next() ;
880 }
881 }
882
883 // some checks
884 vector<uInt> nchans;
885 for (Int i = 0 ; i < in->nrow() ; i++) {
886 nchans.push_back((in->getSpectrum(i)).size());
887 }
888 //Vector<uInt> mchans( nchans ) ;
889 vector< Vector<Float> > facts ;
890 for ( uInt i = 0 ; i < nchans.size() ; i++ ) {
891 Vector<Float> tmp( val[i] ) ;
892 // check divided by zero
893 if ( ( mode == "DIV" ) && anyEQ( tmp, (float)0.0 ) ) {
894 throw( AipsError("Divided by zero is not recommended." ) ) ;
895 }
896 // conformity check
897 if ( tmp.nelements() != nchans[i] ) {
898 stringstream ss ;
899 ss << "Row " << i << ": Vector size must be same as number of channel." ;
900 throw( AipsError( ss.str() ) ) ;
901 }
902 facts.push_back( tmp ) ;
903 }
904
905
906 CountedPtr< Scantable > out = getScantable(in, false);
907 Table& tab = out->table();
908 ArrayColumn<Float> specCol(tab,"SPECTRA");
909 ArrayColumn<Float> tsysCol(tab,"TSYS");
910 for (uInt i=0; i<tab.nrow(); ++i) {
911 Vector<Float> fact = facts[i] ;
912 Vector<Float> spec;
913 Vector<Float> ts;
914 specCol.get(i, spec);
915 tsysCol.get(i, ts);
916 if (mode == "MUL" || mode == "DIV") {
917 if (mode == "DIV") fact = (float)1.0 / fact;
918 spec *= fact;
919 specCol.put(i, spec);
920 if ( tsys ) {
921 ts *= fact;
922 tsysCol.put(i, ts);
923 }
924 } else if ( mode == "ADD" || mode == "SUB") {
925 if (mode == "SUB") fact *= (float)-1.0 ;
926 spec += fact;
927 specCol.put(i, spec);
928 if ( tsys ) {
929 ts += fact;
930 tsysCol.put(i, ts);
931 }
932 }
933 }
934 return out;
935}
936
937CountedPtr<Scantable> STMath::binaryOperate(const CountedPtr<Scantable>& left,
938 const CountedPtr<Scantable>& right,
939 const std::string& mode)
940{
941 bool insitu = insitu_;
942 if ( ! left->conformant(*right) ) {
943 throw(AipsError("'left' and 'right' scantables are not conformant."));
944 }
945 setInsitu(false);
946 CountedPtr< Scantable > out = getScantable(left, false);
947 setInsitu(insitu);
948 Table& tout = out->table();
949 Block<String> coln(5);
950 coln[0] = "SCANNO"; coln[1] = "CYCLENO"; coln[2] = "BEAMNO";
951 coln[3] = "IFNO"; coln[4] = "POLNO";
952 Table tmpl = tout.sort(coln);
953 Table tmpr = right->table().sort(coln);
954 ArrayColumn<Float> lspecCol(tmpl,"SPECTRA");
955 ROArrayColumn<Float> rspecCol(tmpr,"SPECTRA");
956 ArrayColumn<uChar> lflagCol(tmpl,"FLAGTRA");
957 ROArrayColumn<uChar> rflagCol(tmpr,"FLAGTRA");
958
959 for (uInt i=0; i<tout.nrow(); ++i) {
960 Vector<Float> lspecvec, rspecvec;
961 Vector<uChar> lflagvec, rflagvec;
962 lspecvec = lspecCol(i); rspecvec = rspecCol(i);
963 lflagvec = lflagCol(i); rflagvec = rflagCol(i);
964 MaskedArray<Float> mleft = maskedArray(lspecvec, lflagvec);
965 MaskedArray<Float> mright = maskedArray(rspecvec, rflagvec);
966 if (mode == "ADD") {
967 mleft += mright;
968 } else if ( mode == "SUB") {
969 mleft -= mright;
970 } else if ( mode == "MUL") {
971 mleft *= mright;
972 } else if ( mode == "DIV") {
973 mleft /= mright;
974 } else {
975 throw(AipsError("Illegal binary operator"));
976 }
977 lspecCol.put(i, mleft.getArray());
978 }
979 return out;
980}
981
982
983
984MaskedArray<Float> STMath::maskedArray( const Vector<Float>& s,
985 const Vector<uChar>& f)
986{
987 Vector<Bool> mask;
988 mask.resize(f.shape());
989 convertArray(mask, f);
990 return MaskedArray<Float>(s,!mask);
991}
992
993MaskedArray<Double> STMath::maskedArray( const Vector<Double>& s,
994 const Vector<uChar>& f)
995{
996 Vector<Bool> mask;
997 mask.resize(f.shape());
998 convertArray(mask, f);
999 return MaskedArray<Double>(s,!mask);
1000}
1001
1002Vector<uChar> STMath::flagsFromMA(const MaskedArray<Float>& ma)
1003{
1004 const Vector<Bool>& m = ma.getMask();
1005 Vector<uChar> flags(m.shape());
1006 convertArray(flags, !m);
1007 return flags;
1008}
1009
1010CountedPtr< Scantable > STMath::autoQuotient( const CountedPtr< Scantable >& in,
1011 const std::string & mode,
1012 bool preserve )
1013{
1014 /// @todo make other modes available
1015 /// modes should be "nearest", "pair"
1016 // make this operation non insitu
1017 (void) mode; //currently unused
1018 const Table& tin = in->table();
1019 Table ons = tin(tin.col("SRCTYPE") == Int(SrcType::PSON));
1020 Table offs = tin(tin.col("SRCTYPE") == Int(SrcType::PSOFF));
1021 if ( offs.nrow() == 0 )
1022 throw(AipsError("No 'off' scans present."));
1023 // put all "on" scans into output table
1024
1025 bool insitu = insitu_;
1026 setInsitu(false);
1027 CountedPtr< Scantable > out = getScantable(in, true);
1028 setInsitu(insitu);
1029 Table& tout = out->table();
1030
1031 TableCopy::copyRows(tout, ons);
1032 TableRow row(tout);
1033 ROScalarColumn<Double> offtimeCol(offs, "TIME");
1034 ArrayColumn<Float> outspecCol(tout, "SPECTRA");
1035 ROArrayColumn<Float> outtsysCol(tout, "TSYS");
1036 ArrayColumn<uChar> outflagCol(tout, "FLAGTRA");
1037 for (uInt i=0; i < tout.nrow(); ++i) {
1038 const TableRecord& rec = row.get(i);
1039 Double ontime = rec.asDouble("TIME");
1040 Table presel = offs(offs.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
1041 && offs.col("IFNO") == Int(rec.asuInt("IFNO"))
1042 && offs.col("POLNO") == Int(rec.asuInt("POLNO")) );
1043 ROScalarColumn<Double> offtimeCol(presel, "TIME");
1044
1045 Double mindeltat = min(abs(offtimeCol.getColumn() - ontime));
1046 // Timestamp may vary within a cycle ???!!!
1047 // increase this by 0.01 sec in case of rounding errors...
1048 // There might be a better way to do this.
1049 // fix to this fix. TIME is MJD, so 1.0d not 1.0s
1050 mindeltat += 0.01/24./60./60.;
1051 Table sel = presel( abs(presel.col("TIME")-ontime) <= mindeltat);
1052
1053 if ( sel.nrow() < 1 ) {
1054 throw(AipsError("No closest in time found... This could be a rounding "
1055 "issue. Try quotient instead."));
1056 }
1057 TableRow offrow(sel);
1058 const TableRecord& offrec = offrow.get(0);//should only be one row
1059 RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
1060 RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
1061 RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
1062 /// @fixme this assumes tsys is a scalar not vector
1063 Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
1064 Vector<Float> specon, tsyson;
1065 outtsysCol.get(i, tsyson);
1066 outspecCol.get(i, specon);
1067 Vector<uChar> flagon;
1068 outflagCol.get(i, flagon);
1069 MaskedArray<Float> mon = maskedArray(specon, flagon);
1070 MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
1071 MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
1072 if (preserve) {
1073 quot -= tsysoffscalar;
1074 } else {
1075 quot -= tsyson[0];
1076 }
1077 outspecCol.put(i, quot.getArray());
1078 outflagCol.put(i, flagsFromMA(quot));
1079 }
1080 // renumber scanno
1081 TableIterator it(tout, "SCANNO");
1082 uInt i = 0;
1083 while ( !it.pastEnd() ) {
1084 Table t = it.table();
1085 TableVector<uInt> vec(t, "SCANNO");
1086 vec = i;
1087 ++i;
1088 ++it;
1089 }
1090 return out;
1091}
1092
1093
1094CountedPtr< Scantable > STMath::quotient( const CountedPtr< Scantable > & on,
1095 const CountedPtr< Scantable > & off,
1096 bool preserve )
1097{
1098 bool insitu = insitu_;
1099 if ( ! on->conformant(*off) ) {
1100 throw(AipsError("'on' and 'off' scantables are not conformant."));
1101 }
1102 setInsitu(false);
1103 CountedPtr< Scantable > out = getScantable(on, false);
1104 setInsitu(insitu);
1105 Table& tout = out->table();
1106 const Table& toff = off->table();
1107 TableIterator sit(tout, "SCANNO");
1108 TableIterator s2it(toff, "SCANNO");
1109 while ( !sit.pastEnd() ) {
1110 Table ton = sit.table();
1111 TableRow row(ton);
1112 Table t = s2it.table();
1113 ArrayColumn<Float> outspecCol(ton, "SPECTRA");
1114 ROArrayColumn<Float> outtsysCol(ton, "TSYS");
1115 ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
1116 for (uInt i=0; i < ton.nrow(); ++i) {
1117 const TableRecord& rec = row.get(i);
1118 Table offsel = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
1119 && t.col("IFNO") == Int(rec.asuInt("IFNO"))
1120 && t.col("POLNO") == Int(rec.asuInt("POLNO")) );
1121 if ( offsel.nrow() == 0 )
1122 throw AipsError("STMath::quotient: no matching off");
1123 TableRow offrow(offsel);
1124 const TableRecord& offrec = offrow.get(0);//should be ncycles - take first
1125 RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
1126 RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
1127 RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
1128 Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
1129 Vector<Float> specon, tsyson;
1130 outtsysCol.get(i, tsyson);
1131 outspecCol.get(i, specon);
1132 Vector<uChar> flagon;
1133 outflagCol.get(i, flagon);
1134 MaskedArray<Float> mon = maskedArray(specon, flagon);
1135 MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
1136 MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
1137 if (preserve) {
1138 quot -= tsysoffscalar;
1139 } else {
1140 quot -= tsyson[0];
1141 }
1142 outspecCol.put(i, quot.getArray());
1143 outflagCol.put(i, flagsFromMA(quot));
1144 }
1145 ++sit;
1146 ++s2it;
1147 // take the first off for each on scan which doesn't have a
1148 // matching off scan
1149 // non <= noff: matching pairs, non > noff matching pairs then first off
1150 if ( s2it.pastEnd() ) s2it.reset();
1151 }
1152 return out;
1153}
1154
1155// dototalpower (migration of GBTIDL procedure dototalpower.pro)
1156// calibrate the CAL on-off pair. It calculate Tsys and average CAL on-off subintegrations
1157// do it for each cycles in a specific scan.
1158CountedPtr< Scantable > STMath::dototalpower( const CountedPtr< Scantable >& calon,
1159 const CountedPtr< Scantable >& caloff, Float tcal )
1160{
1161 if ( ! calon->conformant(*caloff) ) {
1162 throw(AipsError("'CAL on' and 'CAL off' scantables are not conformant."));
1163 }
1164 setInsitu(false);
1165 CountedPtr< Scantable > out = getScantable(caloff, false);
1166 Table& tout = out->table();
1167 const Table& tcon = calon->table();
1168 Vector<Float> tcalout;
1169
1170 std::map<uInt,uInt> tcalIdToRecNoMap;
1171 const Table& calOffTcalTable = caloff->tcal().table();
1172 {
1173 ROScalarColumn<uInt> calOffTcalTable_IDcol(calOffTcalTable, "ID");
1174 const Vector<uInt> tcalIds(calOffTcalTable_IDcol.getColumn());
1175 size_t tcalIdsEnd = tcalIds.nelements();
1176 for (uInt i = 0; i < tcalIdsEnd; i++) {
1177 tcalIdToRecNoMap[tcalIds[i]] = i;
1178 }
1179 }
1180 ROArrayColumn<Float> calOffTcalTable_TCALcol(calOffTcalTable, "TCAL");
1181
1182 if ( tout.nrow() != tcon.nrow() ) {
1183 throw(AipsError("Mismatch in number of rows to form cal on - off pair."));
1184 }
1185 // iteration by scanno or cycle no.
1186 TableIterator sit(tout, "SCANNO");
1187 TableIterator s2it(tcon, "SCANNO");
1188 while ( !sit.pastEnd() ) {
1189 Table toff = sit.table();
1190 TableRow row(toff);
1191 Table t = s2it.table();
1192 ScalarColumn<Double> outintCol(toff, "INTERVAL");
1193 ArrayColumn<Float> outspecCol(toff, "SPECTRA");
1194 ArrayColumn<Float> outtsysCol(toff, "TSYS");
1195 ArrayColumn<uChar> outflagCol(toff, "FLAGTRA");
1196 ROScalarColumn<uInt> outtcalIdCol(toff, "TCAL_ID");
1197 ROScalarColumn<uInt> outpolCol(toff, "POLNO");
1198 ROScalarColumn<Double> onintCol(t, "INTERVAL");
1199 ROArrayColumn<Float> onspecCol(t, "SPECTRA");
1200 ROArrayColumn<Float> ontsysCol(t, "TSYS");
1201 ROArrayColumn<uChar> onflagCol(t, "FLAGTRA");
1202 //ROScalarColumn<uInt> ontcalIdCol(t, "TCAL_ID");
1203
1204 for (uInt i=0; i < toff.nrow(); ++i) {
1205 //skip these checks -> assumes the data order are the same between the cal on off pairs
1206 //
1207 Vector<Float> specCalon, specCaloff;
1208 // to store scalar (mean) tsys
1209 Vector<Float> tsysout(1);
1210 uInt tcalId, polno;
1211 Double offint, onint;
1212 outpolCol.get(i, polno);
1213 outspecCol.get(i, specCaloff);
1214 onspecCol.get(i, specCalon);
1215 Vector<uChar> flagCaloff, flagCalon;
1216 outflagCol.get(i, flagCaloff);
1217 onflagCol.get(i, flagCalon);
1218 outtcalIdCol.get(i, tcalId);
1219 outintCol.get(i, offint);
1220 onintCol.get(i, onint);
1221 // caluculate mean Tsys
1222 uInt nchan = specCaloff.nelements();
1223 // percentage of edge cut off
1224 uInt pc = 10;
1225 uInt bchan = nchan/pc;
1226 uInt echan = nchan-bchan;
1227
1228 Slicer chansl(IPosition(1,bchan-1), IPosition(1,echan-1), IPosition(1,1),Slicer::endIsLast);
1229 Vector<Float> testsubsp = specCaloff(chansl);
1230 MaskedArray<Float> spoff = maskedArray( specCaloff(chansl),flagCaloff(chansl) );
1231 MaskedArray<Float> spon = maskedArray( specCalon(chansl),flagCalon(chansl) );
1232 MaskedArray<Float> spdiff = spon-spoff;
1233 uInt noff = spoff.nelementsValid();
1234 //uInt non = spon.nelementsValid();
1235 uInt ndiff = spdiff.nelementsValid();
1236 Float meantsys;
1237
1238/**
1239 Double subspec, subdiff;
1240 uInt usednchan;
1241 subspec = 0;
1242 subdiff = 0;
1243 usednchan = 0;
1244 for(uInt k=(bchan-1); k<echan; k++) {
1245 subspec += specCaloff[k];
1246 subdiff += static_cast<Double>(specCalon[k]-specCaloff[k]);
1247 ++usednchan;
1248 }
1249**/
1250 // get tcal if input tcal <= 0
1251 Float tcalUsed;
1252 tcalUsed = tcal;
1253 if ( tcal <= 0.0 ) {
1254 uInt tcalRecNo = tcalIdToRecNoMap[tcalId];
1255 calOffTcalTable_TCALcol.get(tcalRecNo, tcalout);
1256// if (polno<=3) {
1257// tcalUsed = tcalout[polno];
1258// }
1259// else {
1260// tcalUsed = tcalout[0];
1261// }
1262 if ( tcalout.size() == 1 )
1263 tcalUsed = tcalout[0] ;
1264 else if ( tcalout.size() == nchan )
1265 tcalUsed = mean(tcalout) ;
1266 else {
1267 uInt ipol = polno ;
1268 if ( ipol > 3 ) ipol = 0 ;
1269 tcalUsed = tcalout[ipol] ;
1270 }
1271 }
1272
1273 Float meanoff;
1274 Float meandiff;
1275 if (noff && ndiff) {
1276 //Debug
1277 //if(noff!=ndiff) cerr<<"noff and ndiff is not equal"<<endl;
1278 //LogIO os( LogOrigin( "STMath", "dototalpower()", WHERE ) ) ;
1279 //if(noff!=ndiff) os<<"noff and ndiff is not equal"<<LogIO::POST;
1280 meanoff = sum(spoff)/noff;
1281 meandiff = sum(spdiff)/ndiff;
1282 meantsys= (meanoff/meandiff )*tcalUsed + tcalUsed/2;
1283 }
1284 else {
1285 meantsys=1;
1286 }
1287
1288 tsysout[0] = Float(meantsys);
1289 MaskedArray<Float> mcaloff = maskedArray(specCaloff, flagCaloff);
1290 MaskedArray<Float> mcalon = maskedArray(specCalon, flagCalon);
1291 MaskedArray<Float> sig = Float(0.5) * (mcaloff + mcalon);
1292 //uInt ncaloff = mcaloff.nelementsValid();
1293 //uInt ncalon = mcalon.nelementsValid();
1294
1295 outintCol.put(i, offint+onint);
1296 outspecCol.put(i, sig.getArray());
1297 outflagCol.put(i, flagsFromMA(sig));
1298 outtsysCol.put(i, tsysout);
1299 }
1300 ++sit;
1301 ++s2it;
1302 }
1303 return out;
1304}
1305
1306//dosigref - migrated from GBT IDL's dosigref.pro, do calibration of position switch
1307// observatiions.
1308// input: sig and ref scantables, and an optional boxcar smoothing width(default width=0,
1309// no smoothing).
1310// output: resultant scantable [= (sig-ref/ref)*tsys]
1311CountedPtr< Scantable > STMath::dosigref( const CountedPtr < Scantable >& sig,
1312 const CountedPtr < Scantable >& ref,
1313 int smoothref,
1314 casa::Float tsysv,
1315 casa::Float tau )
1316{
1317 LogIO os( casa::LogOrigin( "STMath", "dosigref()"));
1318if ( ! ref->conformant(*sig) ) {
1319 throw(AipsError("'sig' and 'ref' scantables are not conformant."));
1320 }
1321 setInsitu(false);
1322 CountedPtr< Scantable > out = getScantable(sig, false);
1323 CountedPtr< Scantable > smref;
1324 if ( smoothref > 1 ) {
1325 float fsmoothref = static_cast<float>(smoothref);
1326 std::string inkernel = "boxcar";
1327 smref = smooth(ref, inkernel, fsmoothref );
1328 ostringstream oss;
1329 os <<"Applied smoothing of "<<fsmoothref<<" on the reference."
1330 << LogIO::POST;
1331 }
1332 else {
1333 smref = ref;
1334 }
1335 Table& tout = out->table();
1336 const Table& tref = smref->table();
1337 if ( tout.nrow() != tref.nrow() ) {
1338 throw(AipsError("Mismatch in number of rows to form on-source and reference pair."));
1339 }
1340 // iteration by scanno? or cycle no.
1341 TableIterator sit(tout, "SCANNO");
1342 TableIterator s2it(tref, "SCANNO");
1343 while ( !sit.pastEnd() ) {
1344 Table ton = sit.table();
1345 Table t = s2it.table();
1346 ScalarColumn<Double> outintCol(ton, "INTERVAL");
1347 ArrayColumn<Float> outspecCol(ton, "SPECTRA");
1348 ArrayColumn<Float> outtsysCol(ton, "TSYS");
1349 ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
1350 ArrayColumn<Float> refspecCol(t, "SPECTRA");
1351 ROScalarColumn<Double> refintCol(t, "INTERVAL");
1352 ROArrayColumn<Float> reftsysCol(t, "TSYS");
1353 ArrayColumn<uChar> refflagCol(t, "FLAGTRA");
1354 ROScalarColumn<Float> refelevCol(t, "ELEVATION");
1355 for (uInt i=0; i < ton.nrow(); ++i) {
1356
1357 Double onint, refint;
1358 Vector<Float> specon, specref;
1359 // to store scalar (mean) tsys
1360 Vector<Float> tsysref;
1361 outintCol.get(i, onint);
1362 refintCol.get(i, refint);
1363 outspecCol.get(i, specon);
1364 refspecCol.get(i, specref);
1365 Vector<uChar> flagref, flagon;
1366 outflagCol.get(i, flagon);
1367 refflagCol.get(i, flagref);
1368 reftsysCol.get(i, tsysref);
1369
1370 Float tsysrefscalar;
1371 if ( tsysv > 0.0 ) {
1372 ostringstream oss;
1373 Float elev;
1374 refelevCol.get(i, elev);
1375 os << "user specified Tsys = " << tsysv;
1376 // do recalc elevation if EL = 0
1377 if ( elev == 0 ) {
1378 throw(AipsError("EL=0, elevation data is missing."));
1379 } else {
1380 if ( tau <= 0.0 ) {
1381 throw(AipsError("Valid tau is not supplied."));
1382 } else {
1383 tsysrefscalar = tsysv * exp(tau/elev);
1384 }
1385 }
1386 os << ", corrected (for El) tsys= "<<tsysrefscalar;
1387 }
1388 else {
1389 tsysrefscalar = tsysref[0];
1390 }
1391 //get quotient spectrum
1392 MaskedArray<Float> mref = maskedArray(specref, flagref);
1393 MaskedArray<Float> mon = maskedArray(specon, flagon);
1394 MaskedArray<Float> specres = tsysrefscalar*((mon - mref)/mref);
1395 Double resint = onint*refint*smoothref/(onint+refint*smoothref);
1396
1397 //Debug
1398 //cerr<<"Tsys used="<<tsysrefscalar<<endl;
1399 //LogIO os( LogOrigin( "STMath", "dosigref", WHERE ) ) ;
1400 //os<<"Tsys used="<<tsysrefscalar<<LogIO::POST;
1401 // fill the result, replay signal tsys by reference tsys
1402 outintCol.put(i, resint);
1403 outspecCol.put(i, specres.getArray());
1404 outflagCol.put(i, flagsFromMA(specres));
1405 outtsysCol.put(i, tsysref);
1406 }
1407 ++sit;
1408 ++s2it;
1409 }
1410
1411 out->setFluxUnit("K");
1412
1413 return out;
1414}
1415
1416CountedPtr< Scantable > STMath::donod(const casa::CountedPtr<Scantable>& s,
1417 const std::vector<int>& scans,
1418 int smoothref,
1419 casa::Float tsysv,
1420 casa::Float tau,
1421 casa::Float tcal )
1422
1423{
1424 setInsitu(false);
1425 LogIO os( casa::LogOrigin( "STMath", "donod()"));
1426 STSelector sel;
1427 std::vector<int> scan1, scan2, beams, types;
1428 std::vector< vector<int> > scanpair;
1429 //std::vector<string> calstate;
1430 std::vector<int> calstate;
1431 String msg;
1432
1433 CountedPtr< Scantable > s1b1on, s1b1off, s1b2on, s1b2off;
1434 CountedPtr< Scantable > s2b1on, s2b1off, s2b2on, s2b2off;
1435
1436 std::vector< CountedPtr< Scantable > > sctables;
1437 sctables.push_back(s1b1on);
1438 sctables.push_back(s1b1off);
1439 sctables.push_back(s1b2on);
1440 sctables.push_back(s1b2off);
1441 sctables.push_back(s2b1on);
1442 sctables.push_back(s2b1off);
1443 sctables.push_back(s2b2on);
1444 sctables.push_back(s2b2off);
1445
1446 //check scanlist
1447 int n=s->checkScanInfo(scans);
1448 if (n==1) {
1449 throw(AipsError("Incorrect scan pairs. "));
1450 }
1451
1452 // Assume scans contain only a pair of consecutive scan numbers.
1453 // It is assumed that first beam, b1, is on target.
1454 // There is no check if the first beam is on or not.
1455 if ( scans.size()==1 ) {
1456 scan1.push_back(scans[0]);
1457 scan2.push_back(scans[0]+1);
1458 } else if ( scans.size()==2 ) {
1459 scan1.push_back(scans[0]);
1460 scan2.push_back(scans[1]);
1461 } else {
1462 if ( scans.size()%2 == 0 ) {
1463 for (uInt i=0; i<scans.size(); i++) {
1464 if (i%2 == 0) {
1465 scan1.push_back(scans[i]);
1466 }
1467 else {
1468 scan2.push_back(scans[i]);
1469 }
1470 }
1471 } else {
1472 throw(AipsError("Odd numbers of scans, cannot form pairs."));
1473 }
1474 }
1475 scanpair.push_back(scan1);
1476 scanpair.push_back(scan2);
1477 //calstate.push_back("*calon");
1478 //calstate.push_back("*[^calon]");
1479 calstate.push_back(SrcType::NODCAL);
1480 calstate.push_back(SrcType::NOD);
1481 CountedPtr< Scantable > ws = getScantable(s, false);
1482 uInt l=0;
1483 while ( l < sctables.size() ) {
1484 for (uInt i=0; i < 2; i++) {
1485 for (uInt j=0; j < 2; j++) {
1486 for (uInt k=0; k < 2; k++) {
1487 sel.reset();
1488 sel.setScans(scanpair[i]);
1489 //sel.setName(calstate[k]);
1490 types.clear();
1491 types.push_back(calstate[k]);
1492 sel.setTypes(types);
1493 beams.clear();
1494 beams.push_back(j);
1495 sel.setBeams(beams);
1496 ws->setSelection(sel);
1497 sctables[l]= getScantable(ws, false);
1498 l++;
1499 }
1500 }
1501 }
1502 }
1503
1504 // replace here by splitData or getData functionality
1505 CountedPtr< Scantable > sig1;
1506 CountedPtr< Scantable > ref1;
1507 CountedPtr< Scantable > sig2;
1508 CountedPtr< Scantable > ref2;
1509 CountedPtr< Scantable > calb1;
1510 CountedPtr< Scantable > calb2;
1511
1512 msg=String("Processing dototalpower for subset of the data");
1513 os << msg << LogIO::POST;
1514 // Debug for IRC CS data
1515 //float tcal1=7.0;
1516 //float tcal2=4.0;
1517 sig1 = dototalpower(sctables[0], sctables[1], tcal=tcal);
1518 ref1 = dototalpower(sctables[2], sctables[3], tcal=tcal);
1519 ref2 = dototalpower(sctables[4], sctables[5], tcal=tcal);
1520 sig2 = dototalpower(sctables[6], sctables[7], tcal=tcal);
1521
1522 // correction of user-specified tsys for elevation here
1523
1524 // dosigref calibration
1525 msg=String("Processing dosigref for subset of the data");
1526 os << msg << endl;
1527 calb1=dosigref(sig1,ref2,smoothref,tsysv,tau);
1528 calb2=dosigref(sig2,ref1,smoothref,tsysv,tau);
1529
1530 // iteration by scanno or cycle no.
1531 Table& tcalb1 = calb1->table();
1532 Table& tcalb2 = calb2->table();
1533 TableIterator sit(tcalb1, "SCANNO");
1534 TableIterator s2it(tcalb2, "SCANNO");
1535 while ( !sit.pastEnd() ) {
1536 Table t1 = sit.table();
1537 Table t2= s2it.table();
1538 ArrayColumn<Float> outspecCol(t1, "SPECTRA");
1539 ArrayColumn<Float> outtsysCol(t1, "TSYS");
1540 ArrayColumn<uChar> outflagCol(t1, "FLAGTRA");
1541 ScalarColumn<Double> outintCol(t1, "INTERVAL");
1542 ArrayColumn<Float> t2specCol(t2, "SPECTRA");
1543 ROArrayColumn<Float> t2tsysCol(t2, "TSYS");
1544 ArrayColumn<uChar> t2flagCol(t2, "FLAGTRA");
1545 ROScalarColumn<Double> t2intCol(t2, "INTERVAL");
1546 for (uInt i=0; i < t1.nrow(); ++i) {
1547 Vector<Float> spec1, spec2;
1548 // to store scalar (mean) tsys
1549 Vector<Float> tsys1, tsys2;
1550 Vector<uChar> flag1, flag2;
1551 Double tint1, tint2;
1552 outspecCol.get(i, spec1);
1553 t2specCol.get(i, spec2);
1554 outflagCol.get(i, flag1);
1555 t2flagCol.get(i, flag2);
1556 outtsysCol.get(i, tsys1);
1557 t2tsysCol.get(i, tsys2);
1558 outintCol.get(i, tint1);
1559 t2intCol.get(i, tint2);
1560 // average
1561 // assume scalar tsys for weights
1562 Float wt1, wt2, tsyssq1, tsyssq2;
1563 tsyssq1 = tsys1[0]*tsys1[0];
1564 tsyssq2 = tsys2[0]*tsys2[0];
1565 wt1 = Float(tint1)/tsyssq1;
1566 wt2 = Float(tint2)/tsyssq2;
1567 Float invsumwt=1/(wt1+wt2);
1568 MaskedArray<Float> mspec1 = maskedArray(spec1, flag1);
1569 MaskedArray<Float> mspec2 = maskedArray(spec2, flag2);
1570 MaskedArray<Float> avspec = invsumwt * (wt1*mspec1 + wt2*mspec2);
1571 //Array<Float> avtsys = Float(0.5) * (tsys1 + tsys2);
1572 // cerr<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<endl;
1573 // LogIO os( LogOrigin( "STMath", "donod", WHERE ) ) ;
1574 // os<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<LogIO::POST;
1575 tsys1[0] = sqrt(tsyssq1 + tsyssq2);
1576 Array<Float> avtsys = tsys1;
1577
1578 outspecCol.put(i, avspec.getArray());
1579 outflagCol.put(i, flagsFromMA(avspec));
1580 outtsysCol.put(i, avtsys);
1581 }
1582 ++sit;
1583 ++s2it;
1584 }
1585
1586 calb1->setFluxUnit("K");
1587
1588 return calb1;
1589}
1590
1591//GBTIDL version of frequency switched data calibration
1592CountedPtr< Scantable > STMath::dofs( const CountedPtr< Scantable >& s,
1593 const std::vector<int>& scans,
1594 int smoothref,
1595 casa::Float tsysv,
1596 casa::Float tau,
1597 casa::Float tcal )
1598{
1599
1600
1601 (void) scans; //currently unused
1602 STSelector sel;
1603 CountedPtr< Scantable > ws = getScantable(s, false);
1604 CountedPtr< Scantable > sig, sigwcal, ref, refwcal;
1605 CountedPtr< Scantable > calsig, calref, out, out1, out2;
1606 Bool nofold=False;
1607 vector<int> types ;
1608
1609 //split the data
1610 //sel.setName("*_fs");
1611 types.push_back( SrcType::FSON ) ;
1612 sel.setTypes( types ) ;
1613 ws->setSelection(sel);
1614 sig = getScantable(ws,false);
1615 sel.reset();
1616 types.clear() ;
1617 //sel.setName("*_fs_calon");
1618 types.push_back( SrcType::FONCAL ) ;
1619 sel.setTypes( types ) ;
1620 ws->setSelection(sel);
1621 sigwcal = getScantable(ws,false);
1622 sel.reset();
1623 types.clear() ;
1624 //sel.setName("*_fsr");
1625 types.push_back( SrcType::FSOFF ) ;
1626 sel.setTypes( types ) ;
1627 ws->setSelection(sel);
1628 ref = getScantable(ws,false);
1629 sel.reset();
1630 types.clear() ;
1631 //sel.setName("*_fsr_calon");
1632 types.push_back( SrcType::FOFFCAL ) ;
1633 sel.setTypes( types ) ;
1634 ws->setSelection(sel);
1635 refwcal = getScantable(ws,false);
1636 sel.reset() ;
1637 types.clear() ;
1638
1639 calsig = dototalpower(sigwcal, sig, tcal=tcal);
1640 calref = dototalpower(refwcal, ref, tcal=tcal);
1641
1642 out1=dosigref(calsig,calref,smoothref,tsysv,tau);
1643 out2=dosigref(calref,calsig,smoothref,tsysv,tau);
1644
1645 Table& tabout1=out1->table();
1646 Table& tabout2=out2->table();
1647 ROScalarColumn<uInt> freqidCol1(tabout1, "FREQ_ID");
1648 ScalarColumn<uInt> freqidCol2(tabout2, "FREQ_ID");
1649 ROArrayColumn<Float> specCol(tabout2, "SPECTRA");
1650 Vector<Float> spec; specCol.get(0, spec);
1651 uInt nchan = spec.nelements();
1652 uInt freqid1; freqidCol1.get(0,freqid1);
1653 uInt freqid2; freqidCol2.get(0,freqid2);
1654 Double rp1, rp2, rv1, rv2, inc1, inc2;
1655 out1->frequencies().getEntry(rp1, rv1, inc1, freqid1);
1656 out2->frequencies().getEntry(rp2, rv2, inc2, freqid2);
1657 //cerr << out1->frequencies().table().nrow() << " " << out2->frequencies().table().nrow() << endl ;
1658 //LogIO os( LogOrigin( "STMath", "dofs()", WHERE ) ) ;
1659 //os << out1->frequencies().table().nrow() << " " << out2->frequencies().table().nrow() << LogIO::POST ;
1660 if (rp1==rp2) {
1661 Double foffset = rv1 - rv2;
1662 uInt choffset = static_cast<uInt>(foffset/abs(inc2));
1663 if (choffset >= nchan) {
1664 //cerr<<"out-band frequency switching, no folding"<<endl;
1665 LogIO os( LogOrigin( "STMath", "dofs()", WHERE ) ) ;
1666 os<<"out-band frequency switching, no folding"<<LogIO::POST;
1667 nofold = True;
1668 }
1669 }
1670
1671 if (nofold) {
1672 std::vector< CountedPtr< Scantable > > tabs;
1673 tabs.push_back(out1);
1674 tabs.push_back(out2);
1675 out = merge(tabs);
1676 }
1677 else {
1678 //out = out1;
1679 Double choffset = ( rv1 - rv2 ) / inc2 ;
1680 out = dofold( out1, out2, choffset ) ;
1681 }
1682
1683 out->setFluxUnit("K");
1684
1685 return out;
1686}
1687
1688CountedPtr<Scantable> STMath::dofold( const CountedPtr<Scantable> &sig,
1689 const CountedPtr<Scantable> &ref,
1690 Double choffset,
1691 Double choffset2 )
1692{
1693 LogIO os( LogOrigin( "STMath", "dofold", WHERE ) ) ;
1694 os << "choffset=" << choffset << " choffset2=" << choffset2 << LogIO::POST ;
1695
1696 // output scantable
1697 CountedPtr<Scantable> out = getScantable( sig, false ) ;
1698
1699 // separate choffset to integer part and decimal part
1700 Int ioffset = (Int)choffset ;
1701 Double doffset = choffset - ioffset ;
1702 Int ioffset2 = (Int)choffset2 ;
1703 Double doffset2 = choffset2 - ioffset2 ;
1704 os << "ioffset=" << ioffset << " doffset=" << doffset << LogIO::POST ;
1705 os << "ioffset2=" << ioffset2 << " doffset2=" << doffset2 << LogIO::POST ;
1706
1707 // get column
1708 ROArrayColumn<Float> specCol1( sig->table(), "SPECTRA" ) ;
1709 ROArrayColumn<Float> specCol2( ref->table(), "SPECTRA" ) ;
1710 ROArrayColumn<Float> tsysCol1( sig->table(), "TSYS" ) ;
1711 ROArrayColumn<Float> tsysCol2( ref->table(), "TSYS" ) ;
1712 ROArrayColumn<uChar> flagCol1( sig->table(), "FLAGTRA" ) ;
1713 ROArrayColumn<uChar> flagCol2( ref->table(), "FLAGTRA" ) ;
1714 ROScalarColumn<Double> mjdCol1( sig->table(), "TIME" ) ;
1715 ROScalarColumn<Double> mjdCol2( ref->table(), "TIME" ) ;
1716 ROScalarColumn<Double> intervalCol1( sig->table(), "INTERVAL" ) ;
1717 ROScalarColumn<Double> intervalCol2( ref->table(), "INTERVAL" ) ;
1718
1719 // check
1720 if ( ioffset == 0 ) {
1721 LogIO os( LogOrigin( "STMath", "dofold()", WHERE ) ) ;
1722 os << "channel offset is zero, no folding" << LogIO::POST ;
1723 return out ;
1724 }
1725 int nchan = ref->nchan() ;
1726 if ( abs(ioffset) >= nchan ) {
1727 LogIO os( LogOrigin( "STMath", "dofold()", WHERE ) ) ;
1728 os << "out-band frequency switching, no folding" << LogIO::POST ;
1729 return out ;
1730 }
1731
1732 // attach column for output scantable
1733 ArrayColumn<Float> specColOut( out->table(), "SPECTRA" ) ;
1734 ArrayColumn<uChar> flagColOut( out->table(), "FLAGTRA" ) ;
1735 ArrayColumn<Float> tsysColOut( out->table(), "TSYS" ) ;
1736 ScalarColumn<Double> mjdColOut( out->table(), "TIME" ) ;
1737 ScalarColumn<Double> intervalColOut( out->table(), "INTERVAL" ) ;
1738 ScalarColumn<uInt> fidColOut( out->table(), "FREQ_ID" ) ;
1739
1740 // for each row
1741 // assume that the data order are same between sig and ref
1742 RowAccumulator acc( asap::W_TINTSYS ) ;
1743 for ( int i = 0 ; i < sig->nrow() ; i++ ) {
1744 // get values
1745 Vector<Float> spsig ;
1746 specCol1.get( i, spsig ) ;
1747 Vector<Float> spref ;
1748 specCol2.get( i, spref ) ;
1749 Vector<Float> tsyssig ;
1750 tsysCol1.get( i, tsyssig ) ;
1751 Vector<Float> tsysref ;
1752 tsysCol2.get( i, tsysref ) ;
1753 Vector<uChar> flagsig ;
1754 flagCol1.get( i, flagsig ) ;
1755 Vector<uChar> flagref ;
1756 flagCol2.get( i, flagref ) ;
1757 Double timesig ;
1758 mjdCol1.get( i, timesig ) ;
1759 Double timeref ;
1760 mjdCol2.get( i, timeref ) ;
1761 Double intsig ;
1762 intervalCol1.get( i, intsig ) ;
1763 Double intref ;
1764 intervalCol2.get( i, intref ) ;
1765
1766 // shift reference spectra
1767 int refchan = spref.nelements() ;
1768 Vector<Float> sspref( spref.nelements() ) ;
1769 Vector<Float> stsysref( tsysref.nelements() ) ;
1770 Vector<uChar> sflagref( flagref.nelements() ) ;
1771 if ( ioffset > 0 ) {
1772 // SPECTRA and FLAGTRA
1773 for ( int j = 0 ; j < refchan-ioffset ; j++ ) {
1774 sspref[j] = spref[j+ioffset] ;
1775 sflagref[j] = flagref[j+ioffset] ;
1776 }
1777 for ( int j = refchan-ioffset ; j < refchan ; j++ ) {
1778 sspref[j] = spref[j-refchan+ioffset] ;
1779 sflagref[j] = flagref[j-refchan+ioffset] ;
1780 }
1781 spref = sspref.copy() ;
1782 flagref = sflagref.copy() ;
1783 for ( int j = 0 ; j < refchan - 1 ; j++ ) {
1784 sspref[j] = doffset * spref[j+1] + ( 1.0 - doffset ) * spref[j] ;
1785 sflagref[j] = flagref[j+1] + flagref[j] ;
1786 }
1787 sspref[refchan-1] = doffset * spref[0] + ( 1.0 - doffset ) * spref[refchan-1] ;
1788 sflagref[refchan-1] = flagref[0] + flagref[refchan-1] ;
1789
1790 // TSYS
1791 if ( spref.nelements() == tsysref.nelements() ) {
1792 for ( int j = 0 ; j < refchan-ioffset ; j++ ) {
1793 stsysref[j] = tsysref[j+ioffset] ;
1794 }
1795 for ( int j = refchan-ioffset ; j < refchan ; j++ ) {
1796 stsysref[j] = tsysref[j-refchan+ioffset] ;
1797 }
1798 tsysref = stsysref.copy() ;
1799 for ( int j = 0 ; j < refchan - 1 ; j++ ) {
1800 stsysref[j] = doffset * tsysref[j+1] + ( 1.0 - doffset ) * tsysref[j] ;
1801 }
1802 stsysref[refchan-1] = doffset * tsysref[0] + ( 1.0 - doffset ) * tsysref[refchan-1] ;
1803 }
1804 }
1805 else {
1806 // SPECTRA and FLAGTRA
1807 for ( int j = 0 ; j < abs(ioffset) ; j++ ) {
1808 sspref[j] = spref[refchan+ioffset+j] ;
1809 sflagref[j] = flagref[refchan+ioffset+j] ;
1810 }
1811 for ( int j = abs(ioffset) ; j < refchan ; j++ ) {
1812 sspref[j] = spref[j+ioffset] ;
1813 sflagref[j] = flagref[j+ioffset] ;
1814 }
1815 spref = sspref.copy() ;
1816 flagref = sflagref.copy() ;
1817 sspref[0] = doffset * spref[refchan-1] + ( 1.0 - doffset ) * spref[0] ;
1818 sflagref[0] = flagref[0] + flagref[refchan-1] ;
1819 for ( int j = 1 ; j < refchan ; j++ ) {
1820 sspref[j] = doffset * spref[j-1] + ( 1.0 - doffset ) * spref[j] ;
1821 sflagref[j] = flagref[j-1] + flagref[j] ;
1822 }
1823 // TSYS
1824 if ( spref.nelements() == tsysref.nelements() ) {
1825 for ( int j = 0 ; j < abs(ioffset) ; j++ ) {
1826 stsysref[j] = tsysref[refchan+ioffset+j] ;
1827 }
1828 for ( int j = abs(ioffset) ; j < refchan ; j++ ) {
1829 stsysref[j] = tsysref[j+ioffset] ;
1830 }
1831 tsysref = stsysref.copy() ;
1832 stsysref[0] = doffset * tsysref[refchan-1] + ( 1.0 - doffset ) * tsysref[0] ;
1833 for ( int j = 1 ; j < refchan ; j++ ) {
1834 stsysref[j] = doffset * tsysref[j-1] + ( 1.0 - doffset ) * tsysref[j] ;
1835 }
1836 }
1837 }
1838
1839 // shift signal spectra if necessary (only for APEX?)
1840 if ( choffset2 != 0.0 ) {
1841 int sigchan = spsig.nelements() ;
1842 Vector<Float> sspsig( spsig.nelements() ) ;
1843 Vector<Float> stsyssig( tsyssig.nelements() ) ;
1844 Vector<uChar> sflagsig( flagsig.nelements() ) ;
1845 if ( ioffset2 > 0 ) {
1846 // SPECTRA and FLAGTRA
1847 for ( int j = 0 ; j < sigchan-ioffset2 ; j++ ) {
1848 sspsig[j] = spsig[j+ioffset2] ;
1849 sflagsig[j] = flagsig[j+ioffset2] ;
1850 }
1851 for ( int j = sigchan-ioffset2 ; j < sigchan ; j++ ) {
1852 sspsig[j] = spsig[j-sigchan+ioffset2] ;
1853 sflagsig[j] = flagsig[j-sigchan+ioffset2] ;
1854 }
1855 spsig = sspsig.copy() ;
1856 flagsig = sflagsig.copy() ;
1857 for ( int j = 0 ; j < sigchan - 1 ; j++ ) {
1858 sspsig[j] = doffset2 * spsig[j+1] + ( 1.0 - doffset2 ) * spsig[j] ;
1859 sflagsig[j] = flagsig[j+1] || flagsig[j] ;
1860 }
1861 sspsig[sigchan-1] = doffset2 * spsig[0] + ( 1.0 - doffset2 ) * spsig[sigchan-1] ;
1862 sflagsig[sigchan-1] = flagsig[0] || flagsig[sigchan-1] ;
1863 // TSTS
1864 if ( spsig.nelements() == tsyssig.nelements() ) {
1865 for ( int j = 0 ; j < sigchan-ioffset2 ; j++ ) {
1866 stsyssig[j] = tsyssig[j+ioffset2] ;
1867 }
1868 for ( int j = sigchan-ioffset2 ; j < sigchan ; j++ ) {
1869 stsyssig[j] = tsyssig[j-sigchan+ioffset2] ;
1870 }
1871 tsyssig = stsyssig.copy() ;
1872 for ( int j = 0 ; j < sigchan - 1 ; j++ ) {
1873 stsyssig[j] = doffset2 * tsyssig[j+1] + ( 1.0 - doffset2 ) * tsyssig[j] ;
1874 }
1875 stsyssig[sigchan-1] = doffset2 * tsyssig[0] + ( 1.0 - doffset2 ) * tsyssig[sigchan-1] ;
1876 }
1877 }
1878 else {
1879 // SPECTRA and FLAGTRA
1880 for ( int j = 0 ; j < abs(ioffset2) ; j++ ) {
1881 sspsig[j] = spsig[sigchan+ioffset2+j] ;
1882 sflagsig[j] = flagsig[sigchan+ioffset2+j] ;
1883 }
1884 for ( int j = abs(ioffset2) ; j < sigchan ; j++ ) {
1885 sspsig[j] = spsig[j+ioffset2] ;
1886 sflagsig[j] = flagsig[j+ioffset2] ;
1887 }
1888 spsig = sspsig.copy() ;
1889 flagsig = sflagsig.copy() ;
1890 sspsig[0] = doffset2 * spsig[sigchan-1] + ( 1.0 - doffset2 ) * spsig[0] ;
1891 sflagsig[0] = flagsig[0] + flagsig[sigchan-1] ;
1892 for ( int j = 1 ; j < sigchan ; j++ ) {
1893 sspsig[j] = doffset2 * spsig[j-1] + ( 1.0 - doffset2 ) * spsig[j] ;
1894 sflagsig[j] = flagsig[j-1] + flagsig[j] ;
1895 }
1896 // TSYS
1897 if ( spsig.nelements() == tsyssig.nelements() ) {
1898 for ( int j = 0 ; j < abs(ioffset2) ; j++ ) {
1899 stsyssig[j] = tsyssig[sigchan+ioffset2+j] ;
1900 }
1901 for ( int j = abs(ioffset2) ; j < sigchan ; j++ ) {
1902 stsyssig[j] = tsyssig[j+ioffset2] ;
1903 }
1904 tsyssig = stsyssig.copy() ;
1905 stsyssig[0] = doffset2 * tsyssig[sigchan-1] + ( 1.0 - doffset2 ) * tsyssig[0] ;
1906 for ( int j = 1 ; j < sigchan ; j++ ) {
1907 stsyssig[j] = doffset2 * tsyssig[j-1] + ( 1.0 - doffset2 ) * tsyssig[j] ;
1908 }
1909 }
1910 }
1911 }
1912
1913 // folding
1914 acc.add( spsig, !flagsig, tsyssig, intsig, timesig ) ;
1915 acc.add( sspref, !sflagref, stsysref, intref, timeref ) ;
1916
1917 // put result
1918 specColOut.put( i, acc.getSpectrum() ) ;
1919 const Vector<Bool> &msk = acc.getMask() ;
1920 Vector<uChar> flg( msk.shape() ) ;
1921 convertArray( flg, !msk ) ;
1922 flagColOut.put( i, flg ) ;
1923 tsysColOut.put( i, acc.getTsys() ) ;
1924 intervalColOut.put( i, acc.getInterval() ) ;
1925 mjdColOut.put( i, acc.getTime() ) ;
1926 // change FREQ_ID to unshifted IF setting (only for APEX?)
1927 if ( choffset2 != 0.0 ) {
1928 uInt freqid = fidColOut( 0 ) ; // assume single-IF data
1929 double refpix, refval, increment ;
1930 out->frequencies().getEntry( refpix, refval, increment, freqid ) ;
1931 refval -= choffset * increment ;
1932 uInt newfreqid = out->frequencies().addEntry( refpix, refval, increment ) ;
1933 Vector<uInt> freqids = fidColOut.getColumn() ;
1934 for ( uInt j = 0 ; j < freqids.nelements() ; j++ ) {
1935 if ( freqids[j] == freqid )
1936 freqids[j] = newfreqid ;
1937 }
1938 fidColOut.putColumn( freqids ) ;
1939 }
1940
1941 acc.reset() ;
1942 }
1943
1944 return out ;
1945}
1946
1947
1948CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
1949{
1950 // make copy or reference
1951 CountedPtr< Scantable > out = getScantable(in, false);
1952 Table& tout = out->table();
1953 Block<String> cols(4);
1954 cols[0] = String("SCANNO");
1955 cols[1] = String("CYCLENO");
1956 cols[2] = String("BEAMNO");
1957 cols[3] = String("POLNO");
1958 TableIterator iter(tout, cols);
1959 while (!iter.pastEnd()) {
1960 Table subt = iter.table();
1961 // this should leave us with two rows for the two IFs....if not ignore
1962 if (subt.nrow() != 2 ) {
1963 continue;
1964 }
1965 ArrayColumn<Float> specCol(subt, "SPECTRA");
1966 ArrayColumn<Float> tsysCol(subt, "TSYS");
1967 ArrayColumn<uChar> flagCol(subt, "FLAGTRA");
1968 Vector<Float> onspec,offspec, ontsys, offtsys;
1969 Vector<uChar> onflag, offflag;
1970 tsysCol.get(0, ontsys); tsysCol.get(1, offtsys);
1971 specCol.get(0, onspec); specCol.get(1, offspec);
1972 flagCol.get(0, onflag); flagCol.get(1, offflag);
1973 MaskedArray<Float> on = maskedArray(onspec, onflag);
1974 MaskedArray<Float> off = maskedArray(offspec, offflag);
1975 MaskedArray<Float> oncopy = on.copy();
1976
1977 on /= off; on -= 1.0f;
1978 on *= ontsys[0];
1979 off /= oncopy; off -= 1.0f;
1980 off *= offtsys[0];
1981 specCol.put(0, on.getArray());
1982 const Vector<Bool>& m0 = on.getMask();
1983 Vector<uChar> flags0(m0.shape());
1984 convertArray(flags0, !m0);
1985 flagCol.put(0, flags0);
1986
1987 specCol.put(1, off.getArray());
1988 const Vector<Bool>& m1 = off.getMask();
1989 Vector<uChar> flags1(m1.shape());
1990 convertArray(flags1, !m1);
1991 flagCol.put(1, flags1);
1992 ++iter;
1993 }
1994
1995 return out;
1996}
1997
1998std::vector< float > STMath::statistic( const CountedPtr< Scantable > & in,
1999 const std::vector< bool > & mask,
2000 const std::string& which )
2001{
2002
2003 Vector<Bool> m(mask);
2004 const Table& tab = in->table();
2005 ROArrayColumn<Float> specCol(tab, "SPECTRA");
2006 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2007 std::vector<float> out;
2008 for (uInt i=0; i < tab.nrow(); ++i ) {
2009 Vector<Float> spec; specCol.get(i, spec);
2010 Vector<uChar> flag; flagCol.get(i, flag);
2011 MaskedArray<Float> ma = maskedArray(spec, flag);
2012 float outstat = 0.0;
2013 if ( spec.nelements() == m.nelements() ) {
2014 outstat = mathutil::statistics(which, ma(m));
2015 } else {
2016 outstat = mathutil::statistics(which, ma);
2017 }
2018 out.push_back(outstat);
2019 }
2020 return out;
2021}
2022
2023std::vector< float > STMath::statisticRow( const CountedPtr< Scantable > & in,
2024 const std::vector< bool > & mask,
2025 const std::string& which,
2026 int row )
2027{
2028
2029 Vector<Bool> m(mask);
2030 const Table& tab = in->table();
2031 ROArrayColumn<Float> specCol(tab, "SPECTRA");
2032 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2033 std::vector<float> out;
2034
2035 Vector<Float> spec; specCol.get(row, spec);
2036 Vector<uChar> flag; flagCol.get(row, flag);
2037 MaskedArray<Float> ma = maskedArray(spec, flag);
2038 float outstat = 0.0;
2039 if ( spec.nelements() == m.nelements() ) {
2040 outstat = mathutil::statistics(which, ma(m));
2041 } else {
2042 outstat = mathutil::statistics(which, ma);
2043 }
2044 out.push_back(outstat);
2045
2046 return out;
2047}
2048
2049std::vector< int > STMath::minMaxChan( const CountedPtr< Scantable > & in,
2050 const std::vector< bool > & mask,
2051 const std::string& which )
2052{
2053
2054 Vector<Bool> m(mask);
2055 const Table& tab = in->table();
2056 ROArrayColumn<Float> specCol(tab, "SPECTRA");
2057 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2058 std::vector<int> out;
2059 for (uInt i=0; i < tab.nrow(); ++i ) {
2060 Vector<Float> spec; specCol.get(i, spec);
2061 Vector<uChar> flag; flagCol.get(i, flag);
2062 MaskedArray<Float> ma = maskedArray(spec, flag);
2063 if (ma.ndim() != 1) {
2064 throw (ArrayError(
2065 "std::vector<int> STMath::minMaxChan("
2066 "ContedPtr<Scantable> &in, std::vector<bool> &mask, "
2067 " std::string &which)"
2068 " - MaskedArray is not 1D"));
2069 }
2070 IPosition outpos(1,0);
2071 if ( spec.nelements() == m.nelements() ) {
2072 outpos = mathutil::minMaxPos(which, ma(m));
2073 } else {
2074 outpos = mathutil::minMaxPos(which, ma);
2075 }
2076 out.push_back(outpos[0]);
2077 }
2078 return out;
2079}
2080
2081CountedPtr< Scantable > STMath::bin( const CountedPtr< Scantable > & in,
2082 int width )
2083{
2084 if ( !in->getSelection().empty() ) throw(AipsError("Can't bin subset of the data."));
2085 CountedPtr< Scantable > out = getScantable(in, false);
2086 Table& tout = out->table();
2087 out->frequencies().rescale(width, "BIN");
2088 ArrayColumn<Float> specCol(tout, "SPECTRA");
2089 ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
2090 ArrayColumn<Float> tsysCol(tout, "TSYS");
2091
2092 for (uInt i=0; i < tout.nrow(); ++i ) {
2093 MaskedArray<Float> main = maskedArray(specCol(i), flagCol(i));
2094 MaskedArray<Float> maout;
2095 LatticeUtilities::bin(maout, main, 0, Int(width));
2096 specCol.put(i, maout.getArray());
2097 flagCol.put(i, flagsFromMA(maout));
2098 if (tsysCol(i).nelements() == specCol(i).nelements()) {
2099 MaskedArray<Float> matsysin = maskedArray(tsysCol(i), flagCol(i));
2100 MaskedArray<Float> matsysout;
2101 LatticeUtilities::bin(matsysout, matsysin, 0, Int(width));
2102 tsysCol.put(i, matsysout.getArray());
2103 }
2104 // take only the first binned spectrum's length for the deprecated
2105 // global header item nChan
2106 if (i==0) tout.rwKeywordSet().define(String("nChan"),
2107 Int(maout.getArray().nelements()));
2108 }
2109 return out;
2110}
2111
2112CountedPtr< Scantable > STMath::resample( const CountedPtr< Scantable >& in,
2113 const std::string& method,
2114 float width )
2115//
2116// Should add the possibility of width being specified in km/s. This means
2117// that for each freqID (SpectralCoordinate) we will need to convert to an
2118// average channel width (say at the reference pixel). Then we would need
2119// to be careful to make sure each spectrum (of different freqID)
2120// is the same length.
2121//
2122{
2123 //InterpolateArray1D<Double,Float>::InterpolationMethod interp;
2124 Int interpMethod(stringToIMethod(method));
2125
2126 CountedPtr< Scantable > out = getScantable(in, false);
2127 Table& tout = out->table();
2128
2129// Resample SpectralCoordinates (one per freqID)
2130 out->frequencies().rescale(width, "RESAMPLE");
2131 TableIterator iter(tout, "IFNO");
2132 TableRow row(tout);
2133 while ( !iter.pastEnd() ) {
2134 Table tab = iter.table();
2135 ArrayColumn<Float> specCol(tab, "SPECTRA");
2136 //ArrayColumn<Float> tsysCol(tout, "TSYS");
2137 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2138 Vector<Float> spec;
2139 Vector<uChar> flag;
2140 specCol.get(0,spec); // the number of channels should be constant per IF
2141 uInt nChanIn = spec.nelements();
2142 Vector<Float> xIn(nChanIn); indgen(xIn);
2143 Int fac = Int(nChanIn/width);
2144 Vector<Float> xOut(fac+10); // 10 to be safe - resize later
2145 uInt k = 0;
2146 Float x = 0.0;
2147 while (x < Float(nChanIn) ) {
2148 xOut(k) = x;
2149 k++;
2150 x += width;
2151 }
2152 uInt nChanOut = k;
2153 xOut.resize(nChanOut, True);
2154 // process all rows for this IFNO
2155 Vector<Float> specOut;
2156 Vector<Bool> maskOut;
2157 Vector<uChar> flagOut;
2158 for (uInt i=0; i < tab.nrow(); ++i) {
2159 specCol.get(i, spec);
2160 flagCol.get(i, flag);
2161 Vector<Bool> mask(flag.nelements());
2162 convertArray(mask, flag);
2163
2164 IPosition shapeIn(spec.shape());
2165 //sh.nchan = nChanOut;
2166 InterpolateArray1D<Float,Float>::interpolate(specOut, maskOut, xOut,
2167 xIn, spec, mask,
2168 interpMethod, True, True);
2169 /// @todo do the same for channel based Tsys
2170 flagOut.resize(maskOut.nelements());
2171 convertArray(flagOut, maskOut);
2172 specCol.put(i, specOut);
2173 flagCol.put(i, flagOut);
2174 }
2175 ++iter;
2176 }
2177
2178 return out;
2179}
2180
2181STMath::imethod STMath::stringToIMethod(const std::string& in)
2182{
2183 static STMath::imap lookup;
2184
2185 // initialize the lookup table if necessary
2186 if ( lookup.empty() ) {
2187 lookup["nearest"] = InterpolateArray1D<Double,Float>::nearestNeighbour;
2188 lookup["linear"] = InterpolateArray1D<Double,Float>::linear;
2189 lookup["cubic"] = InterpolateArray1D<Double,Float>::cubic;
2190 lookup["spline"] = InterpolateArray1D<Double,Float>::spline;
2191 }
2192
2193 STMath::imap::const_iterator iter = lookup.find(in);
2194
2195 if ( lookup.end() == iter ) {
2196 std::string message = in;
2197 message += " is not a valid interpolation mode";
2198 throw(AipsError(message));
2199 }
2200 return iter->second;
2201}
2202
2203WeightType STMath::stringToWeight(const std::string& in)
2204{
2205 static std::map<std::string, WeightType> lookup;
2206
2207 // initialize the lookup table if necessary
2208 if ( lookup.empty() ) {
2209 lookup["NONE"] = asap::W_NONE;
2210 lookup["TINT"] = asap::W_TINT;
2211 lookup["TINTSYS"] = asap::W_TINTSYS;
2212 lookup["TSYS"] = asap::W_TSYS;
2213 lookup["VAR"] = asap::W_VAR;
2214 }
2215
2216 std::map<std::string, WeightType>::const_iterator iter = lookup.find(in);
2217
2218 if ( lookup.end() == iter ) {
2219 std::string message = in;
2220 message += " is not a valid weighting mode";
2221 throw(AipsError(message));
2222 }
2223 return iter->second;
2224}
2225
2226CountedPtr< Scantable > STMath::gainElevation( const CountedPtr< Scantable >& in,
2227 const vector< float > & coeff,
2228 const std::string & filename,
2229 const std::string& method)
2230{
2231 LogIO os( LogOrigin( "STMath", "gainElevation", WHERE ) ) ;
2232 // Get elevation data from Scantable and convert to degrees
2233 CountedPtr< Scantable > out = getScantable(in, false);
2234 Table& tab = out->table();
2235 ROScalarColumn<Float> elev(tab, "ELEVATION");
2236 Vector<Float> x = elev.getColumn();
2237 x *= Float(180 / C::pi); // Degrees
2238
2239 Vector<Float> coeffs(coeff);
2240 const uInt nc = coeffs.nelements();
2241 if ( filename.length() > 0 && nc > 0 ) {
2242 throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
2243 }
2244
2245 // Correct
2246 if ( nc > 0 || filename.length() == 0 ) {
2247 // Find instrument
2248 Bool throwit = True;
2249 Instrument inst =
2250 STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
2251 throwit);
2252
2253 // Set polynomial
2254 Polynomial<Float>* ppoly = 0;
2255 Vector<Float> coeff;
2256 String msg;
2257 if ( nc > 0 ) {
2258 ppoly = new Polynomial<Float>(nc-1);
2259 coeff = coeffs;
2260 msg = String("user");
2261 } else {
2262 STAttr sdAttr;
2263 coeff = sdAttr.gainElevationPoly(inst);
2264 ppoly = new Polynomial<Float>(coeff.nelements()-1);
2265 msg = String("built in");
2266 }
2267
2268 if ( coeff.nelements() > 0 ) {
2269 ppoly->setCoefficients(coeff);
2270 } else {
2271 delete ppoly;
2272 throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
2273 }
2274 os << "Making polynomial correction with " << msg << " coefficients:" << endl;
2275 os << " " << coeff << LogIO::POST;
2276 const uInt nrow = tab.nrow();
2277 Vector<Float> factor(nrow);
2278 for ( uInt i=0; i < nrow; ++i ) {
2279 factor[i] = 1.0 / (*ppoly)(x[i]);
2280 }
2281 delete ppoly;
2282 scaleByVector(tab, factor, true);
2283
2284 } else {
2285 // Read and correct
2286 os << "Making correction from ascii Table" << LogIO::POST;
2287 scaleFromAsciiTable(tab, filename, method, x, true);
2288 }
2289 return out;
2290}
2291
2292void STMath::scaleFromAsciiTable(Table& in, const std::string& filename,
2293 const std::string& method,
2294 const Vector<Float>& xout, bool dotsys)
2295{
2296
2297// Read gain-elevation ascii file data into a Table.
2298
2299 String formatString;
2300 Table tbl = readAsciiTable(formatString, Table::Memory, filename, "", "", False);
2301 scaleFromTable(in, tbl, method, xout, dotsys);
2302}
2303
2304void STMath::scaleFromTable(Table& in,
2305 const Table& table,
2306 const std::string& method,
2307 const Vector<Float>& xout, bool dotsys)
2308{
2309
2310 ROScalarColumn<Float> geElCol(table, "ELEVATION");
2311 ROScalarColumn<Float> geFacCol(table, "FACTOR");
2312 Vector<Float> xin = geElCol.getColumn();
2313 Vector<Float> yin = geFacCol.getColumn();
2314 Vector<Bool> maskin(xin.nelements(),True);
2315
2316 // Interpolate (and extrapolate) with desired method
2317
2318 InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
2319
2320 Vector<Float> yout;
2321 Vector<Bool> maskout;
2322 InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout,
2323 xin, yin, maskin, interp,
2324 True, True);
2325
2326 scaleByVector(in, Float(1.0)/yout, dotsys);
2327}
2328
2329void STMath::scaleByVector( Table& in,
2330 const Vector< Float >& factor,
2331 bool dotsys )
2332{
2333 uInt nrow = in.nrow();
2334 if ( factor.nelements() != nrow ) {
2335 throw(AipsError("factors.nelements() != table.nelements()"));
2336 }
2337 ArrayColumn<Float> specCol(in, "SPECTRA");
2338 ArrayColumn<uChar> flagCol(in, "FLAGTRA");
2339 ArrayColumn<Float> tsysCol(in, "TSYS");
2340 for (uInt i=0; i < nrow; ++i) {
2341 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
2342 ma *= factor[i];
2343 specCol.put(i, ma.getArray());
2344 flagCol.put(i, flagsFromMA(ma));
2345 if ( dotsys ) {
2346 Vector<Float> tsys = tsysCol(i);
2347 tsys *= factor[i];
2348 tsysCol.put(i,tsys);
2349 }
2350 }
2351}
2352
2353CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in,
2354 float d, float etaap,
2355 float jyperk )
2356{
2357 LogIO os( LogOrigin( "STMath", "convertFlux", WHERE ) ) ;
2358
2359 CountedPtr< Scantable > out = getScantable(in, false);
2360 Table& tab = in->table();
2361 Table& outtab = out->table();
2362 Unit fluxUnit(tab.keywordSet().asString("FluxUnit"));
2363 Unit K(String("K"));
2364 Unit JY(String("Jy"));
2365
2366 bool tokelvin = true;
2367 Double cfac = 1.0;
2368
2369 if ( fluxUnit == JY ) {
2370 os << "Converting to K" << LogIO::POST;
2371 Quantum<Double> t(1.0,fluxUnit);
2372 Quantum<Double> t2 = t.get(JY);
2373 cfac = (t2 / t).getValue(); // value to Jy
2374
2375 tokelvin = true;
2376 out->setFluxUnit("K");
2377 } else if ( fluxUnit == K ) {
2378 os << "Converting to Jy" << LogIO::POST;
2379 Quantum<Double> t(1.0,fluxUnit);
2380 Quantum<Double> t2 = t.get(K);
2381 cfac = (t2 / t).getValue(); // value to K
2382
2383 tokelvin = false;
2384 out->setFluxUnit("Jy");
2385 } else {
2386 throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
2387 }
2388 // Make sure input values are converted to either Jy or K first...
2389 Float factor = cfac;
2390
2391 // Select method
2392 if (jyperk > 0.0) {
2393 factor *= jyperk;
2394 if ( tokelvin ) factor = 1.0 / jyperk;
2395 os << "Jy/K = " << jyperk << LogIO::POST;
2396 Vector<Float> factors(outtab.nrow(), factor);
2397 scaleByVector(outtab,factors, false);
2398 } else if ( etaap > 0.0) {
2399 if (d < 0) {
2400 Instrument inst =
2401 STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
2402 True);
2403 STAttr sda;
2404 d = sda.diameter(inst);
2405 }
2406 jyperk = STAttr::findJyPerK(etaap, d);
2407 os << "Jy/K = " << jyperk << LogIO::POST;
2408 factor *= jyperk;
2409 if ( tokelvin ) {
2410 factor = 1.0 / factor;
2411 }
2412 Vector<Float> factors(outtab.nrow(), factor);
2413 scaleByVector(outtab, factors, False);
2414 } else {
2415
2416 // OK now we must deal with automatic look up of values.
2417 // We must also deal with the fact that the factors need
2418 // to be computed per IF and may be different and may
2419 // change per integration.
2420
2421 os <<"Looking up conversion factors" << LogIO::POST;
2422 convertBrightnessUnits(out, tokelvin, cfac);
2423 }
2424
2425 return out;
2426}
2427
2428void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in,
2429 bool tokelvin, float cfac )
2430{
2431 Table& table = in->table();
2432 Instrument inst =
2433 STAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True);
2434 TableIterator iter(table, "FREQ_ID");
2435 STFrequencies stfreqs = in->frequencies();
2436 STAttr sdAtt;
2437 while (!iter.pastEnd()) {
2438 Table tab = iter.table();
2439 ArrayColumn<Float> specCol(tab, "SPECTRA");
2440 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2441 ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID");
2442 MEpoch::ROScalarColumn timeCol(tab, "TIME");
2443
2444 uInt freqid; freqidCol.get(0, freqid);
2445 Vector<Float> tmpspec; specCol.get(0, tmpspec);
2446 // STAttr.JyPerK has a Vector interface... change sometime.
2447 Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements()));
2448 for ( uInt i=0; i<tab.nrow(); ++i) {
2449 Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0];
2450 Float factor = cfac * jyperk;
2451 if ( tokelvin ) factor = Float(1.0) / factor;
2452 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
2453 ma *= factor;
2454 specCol.put(i, ma.getArray());
2455 flagCol.put(i, flagsFromMA(ma));
2456 }
2457 ++iter;
2458 }
2459}
2460
2461CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
2462 const std::vector<float>& tau )
2463{
2464 CountedPtr< Scantable > out = getScantable(in, false);
2465
2466 Table outtab = out->table();
2467
2468 const Int ntau = uInt(tau.size());
2469 std::vector<float>::const_iterator tauit = tau.begin();
2470 AlwaysAssert((ntau == 1 || ntau == in->nif() || ntau == in->nif() * in->npol()),
2471 AipsError);
2472 TableIterator iiter(outtab, "IFNO");
2473 while ( !iiter.pastEnd() ) {
2474 Table itab = iiter.table();
2475 TableIterator piter(itab, "POLNO");
2476 while ( !piter.pastEnd() ) {
2477 Table tab = piter.table();
2478 ROScalarColumn<Float> elev(tab, "ELEVATION");
2479 ArrayColumn<Float> specCol(tab, "SPECTRA");
2480 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2481 ArrayColumn<Float> tsysCol(tab, "TSYS");
2482 for ( uInt i=0; i<tab.nrow(); ++i) {
2483 Float zdist = Float(C::pi_2) - elev(i);
2484 Float factor = exp(*tauit/cos(zdist));
2485 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
2486 ma *= factor;
2487 specCol.put(i, ma.getArray());
2488 flagCol.put(i, flagsFromMA(ma));
2489 Vector<Float> tsys;
2490 tsysCol.get(i, tsys);
2491 tsys *= factor;
2492 tsysCol.put(i, tsys);
2493 }
2494 if (ntau == in->nif()*in->npol() ) {
2495 tauit++;
2496 }
2497 piter++;
2498 }
2499 if (ntau >= in->nif() ) {
2500 tauit++;
2501 }
2502 iiter++;
2503 }
2504 return out;
2505}
2506
2507CountedPtr< Scantable > STMath::smoothOther( const CountedPtr< Scantable >& in,
2508 const std::string& kernel,
2509 float width, int order)
2510{
2511 CountedPtr< Scantable > out = getScantable(in, false);
2512 Table table = out->table();
2513
2514 TableIterator iter(table, "IFNO");
2515 while (!iter.pastEnd()) {
2516 Table tab = iter.table();
2517 ArrayColumn<Float> specCol(tab, "SPECTRA");
2518 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2519 Vector<Float> spec;
2520 Vector<uChar> flag;
2521 for (uInt i = 0; i < tab.nrow(); ++i) {
2522 specCol.get(i, spec);
2523 flagCol.get(i, flag);
2524 Vector<Bool> mask(flag.nelements());
2525 convertArray(mask, flag);
2526 Vector<Float> specout;
2527 Vector<Bool> maskout;
2528 if (kernel == "hanning") {
2529 mathutil::hanning(specout, maskout, spec, !mask);
2530 } else if (kernel == "rmedian") {
2531 mathutil::runningMedian(specout, maskout, spec , mask, width);
2532 } else if (kernel == "poly") {
2533 mathutil::polyfit(specout, maskout, spec, !mask, width, order);
2534 }
2535
2536 for (uInt j = 0; j < flag.nelements(); ++j) {
2537 uChar userFlag = 1 << 7;
2538 if (maskout[j]==True) userFlag = 0 << 7;
2539 flag(j) = userFlag;
2540 }
2541
2542 flagCol.put(i, flag);
2543 specCol.put(i, specout);
2544 }
2545 ++iter;
2546 }
2547 return out;
2548}
2549
2550CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
2551 const std::string& kernel, float width,
2552 int order)
2553{
2554 if (kernel == "rmedian" || kernel == "hanning" || kernel == "poly") {
2555 return smoothOther(in, kernel, width, order);
2556 }
2557 CountedPtr< Scantable > out = getScantable(in, false);
2558 Table& table = out->table();
2559 VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
2560 // same IFNO should have same no of channels
2561 // this saves overhead
2562 TableIterator iter(table, "IFNO");
2563 while (!iter.pastEnd()) {
2564 Table tab = iter.table();
2565 ArrayColumn<Float> specCol(tab, "SPECTRA");
2566 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2567 Vector<Float> spec = specCol( 0 );
2568 uInt nchan = spec.nelements();
2569 Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
2570 Convolver<Float> conv(kvec, IPosition(1,nchan));
2571 Vector<uChar> flag;
2572 Vector<Bool> mask(nchan);
2573 for ( uInt i=0; i<tab.nrow(); ++i) {
2574 specCol.get(i, spec);
2575 flagCol.get(i, flag);
2576 convertArray(mask, flag);
2577 Vector<Float> specout;
2578 //mathutil::replaceMaskByZero(specout, mask);
2579 mathutil::replaceMaskByZero(spec, !mask);
2580 //std::vector<bool> vmask;
2581 //(!mask).tovector(vmask);
2582 //mathutil::doZeroOrderInterpolation(spec, vmask);
2583 conv.linearConv(specout, spec);
2584 specCol.put(i, specout);
2585 }
2586 ++iter;
2587 }
2588 return out;
2589}
2590
2591CountedPtr< Scantable >
2592STMath::merge( const std::vector< CountedPtr < Scantable > >& in,
2593 const std::string &freqTol )
2594{
2595 Double freqTolInHz = 0.0; // default is 0.0Hz (merge only when exact match)
2596 if (freqTol.size() > 0) {
2597 Quantum<Double> freqTolInQuantity;
2598 if (!Quantum<Double>::read(freqTolInQuantity, freqTol)) {
2599 throw(AipsError("Failed to convert freqTol string to quantity"));
2600 }
2601 if (!freqTolInQuantity.isConform("Hz")) {
2602 throw(AipsError("Invalid freqTol string"));
2603 }
2604 freqTolInHz = freqTolInQuantity.getValue("Hz");
2605 LogIO os(LogOrigin("STMath", "merge", WHERE));
2606 os << "frequency tolerance = " << freqTolInHz << "Hz" << LogIO::POST;
2607 }
2608
2609 if ( in.size() < 2 ) {
2610 throw(AipsError("Need at least two scantables to perform a merge."));
2611 }
2612 std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
2613 bool insitu = insitu_;
2614 setInsitu(false);
2615 CountedPtr< Scantable > out = getScantable(*it, false);
2616 setInsitu(insitu);
2617 Table& tout = out->table();
2618 ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
2619 ScalarColumn<uInt> scannocol(tout,"SCANNO"), focusidcol(tout,"FOCUS_ID");
2620 ScalarColumn<uInt> ifnocol(tout, "IFNO");
2621 // Renumber SCANNO to be 0-based
2622 Vector<uInt> scannos = scannocol.getColumn();
2623 uInt offset = min(scannos);
2624 scannos -= offset;
2625 scannocol.putColumn(scannos);
2626 uInt newscanno = max(scannos)+1;
2627 ++it;
2628
2629 // new IFNO
2630 uInt ifnoCounter = max(ifnocol.getColumn()) + 1;
2631
2632 // Here we assume that each IFNO has unique MOLECULE_ID
2633 // molIdMap:
2634 // KEY: IFNO
2635 // VALUE: MOLECULE_ID
2636 map<uInt, uInt> molIdMap;
2637 {
2638 TableIterator ifit(tout, "IFNO");
2639 while (!ifit.pastEnd()) {
2640 ROTableRow row(ifit.table());
2641 const TableRecord& rec = row.get(0);
2642 molIdMap[rec.asuInt("IFNO")] = rec.asuInt("MOLECULE_ID");
2643 ifit.next();
2644 }
2645 }
2646
2647 while ( it != in.end() ){
2648 // Check FREQUENCIES/BASEFRAME
2649 if ( out->frequencies().getFrame(true) != (*it)->frequencies().getFrame(true) ) {
2650 throw(AipsError("BASEFRAME is not identical"));
2651 }
2652
2653 if ( ! (*it)->conformant(*out) ) {
2654 // non conformant.
2655 LogIO os( LogOrigin( "STMath", "merge()", WHERE ) ) ;
2656 os << LogIO::SEVERE << "Can't merge scantables as header informations (any one of AntennaName, Equinox, and FluxUnit) differ." << LogIO::EXCEPTION ;
2657 }
2658 out->appendToHistoryTable((*it)->history());
2659 const Table& tab = (*it)->table();
2660
2661 Block<String> cols(3);
2662 cols[0] = String("FREQ_ID");
2663 cols[1] = String("MOLECULE_ID");
2664 cols[2] = String("FOCUS_ID");
2665
2666 TableIterator scanit(tab, "SCANNO");
2667 while (!scanit.pastEnd()) {
2668 TableIterator subit(scanit.table(), cols);
2669 while ( !subit.pastEnd() ) {
2670 uInt nrow = tout.nrow();
2671 Table thetab = subit.table();
2672 ROTableRow row(thetab);
2673 Vector<uInt> thecolvals(thetab.nrow());
2674 // The selected subset of table should have
2675 // the equal FREQ_ID, MOLECULE_ID, and FOCUS_ID values.
2676 const TableRecord& rec = row.get(0);
2677 tout.addRow(thetab.nrow());
2678 TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
2679
2680 Slicer slice(IPosition(1, nrow), IPosition(1, thetab.nrow()),
2681 Slicer::endIsLength);
2682
2683 // Set the proper FREQ_ID
2684 Double rv,rp,inc;
2685 (*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
2686 uInt id;
2687
2688 // default value is new unique IFNO
2689 uInt newifno = ifnoCounter;
2690 Bool isIfMerged = False;
2691 uInt nchan = rec.asArrayFloat("SPECTRA").shape()[0];
2692 //id = out->frequencies().addEntry(rp, rv, inc);
2693 if ( !out->frequencies().match(rp, rv, inc, freqTolInHz, id) ) {
2694 // add new entry to FREQUENCIES table
2695 id = out->frequencies().addEntry(rp, rv, inc);
2696
2697 // increment counter for IFNO
2698 ifnoCounter++;
2699 }
2700 else {
2701 // should renumber IFNO to be same as existing rows that have same FREQ_ID
2702 LogIO os(LogOrigin("STMath", "merge", WHERE));
2703 Table outFreqIdSelected = tout(tout.col("FREQ_ID") == id);
2704 TableIterator _iter(outFreqIdSelected, "IFNO");
2705 map<uInt, uInt> nchanMap;
2706 while (!_iter.pastEnd()) {
2707 const Table _table = _iter.table();
2708 ROTableRow _row(_table);
2709 const TableRecord &_rec = _row.get(0);
2710 uInt nchan = _rec.asArrayFloat("SPECTRA").shape()[0];
2711 if (nchanMap.find(nchan) != nchanMap.end()) {
2712 throw(AipsError("There are non-unique IFNOs assigned to spectra that have same FREQ_ID and same nchan. Something wrong."));
2713 }
2714 nchanMap[nchan] = _rec.asuInt("IFNO");
2715 _iter.next();
2716 }
2717
2718 os << LogIO::DEBUGGING << "nchanMap for " << id << ":" << LogIO::POST;
2719 for (map<uInt, uInt>::iterator i = nchanMap.begin(); i != nchanMap.end(); ++i) {
2720 os << LogIO::DEBUGGING << "nchanMap[" << i->first << "] = " << i->second << LogIO::POST;
2721 }
2722
2723 if (nchanMap.find(nchan) == nchanMap.end()) {
2724 // increment counter for IFNO
2725 ifnoCounter++;
2726 }
2727 else {
2728 // renumber IFNO to be same as existing value that corresponds to nchan
2729 newifno = nchanMap[nchan];
2730 isIfMerged = True;
2731 }
2732 os << LogIO::DEBUGGING << "newifno = " << newifno << LogIO::POST;
2733 }
2734 thecolvals = id;
2735 freqidcol.putColumnRange(slice, thecolvals);
2736
2737 thecolvals = newifno;
2738 ifnocol.putColumnRange(slice, thecolvals);
2739
2740 // Set the proper MOLECULE_ID
2741 Vector<String> name,fname;Vector<Double> rf;
2742 (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
2743 id = out->molecules().addEntry(rf, name, fname);
2744 if (molIdMap.find(newifno) == molIdMap.end()) {
2745 // add new entry to molIdMap
2746 molIdMap[newifno] = id;
2747 }
2748 if (isIfMerged) {
2749 id = molIdMap[newifno];
2750 }
2751 thecolvals = id;
2752 molidcol.putColumnRange(slice, thecolvals);
2753
2754 // Set the proper FOCUS_ID
2755 Float fpa,frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
2756 (*it)->focus().getEntry(fpa, fax, ftan, frot, fhand, fmount,fuser,
2757 fxy, fxyp, rec.asuInt("FOCUS_ID"));
2758 id = out->focus().addEntry(fpa, fax, ftan, frot, fhand, fmount,fuser,
2759 fxy, fxyp);
2760 thecolvals = id;
2761 focusidcol.putColumnRange(slice, thecolvals);
2762
2763 // Set the proper SCANNO
2764 thecolvals = newscanno;
2765 scannocol.putColumnRange(slice, thecolvals);
2766
2767 ++subit;
2768 }
2769 ++newscanno;
2770 ++scanit;
2771 }
2772 ++it;
2773 }
2774 return out;
2775}
2776
2777CountedPtr< Scantable >
2778 STMath::invertPhase( const CountedPtr < Scantable >& in )
2779{
2780 return applyToPol(in, &STPol::invertPhase, Float(0.0));
2781}
2782
2783CountedPtr< Scantable >
2784 STMath::rotateXYPhase( const CountedPtr < Scantable >& in, float phase )
2785{
2786 return applyToPol(in, &STPol::rotatePhase, Float(phase));
2787}
2788
2789CountedPtr< Scantable >
2790 STMath::rotateLinPolPhase( const CountedPtr < Scantable >& in, float phase )
2791{
2792 return applyToPol(in, &STPol::rotateLinPolPhase, Float(phase));
2793}
2794
2795CountedPtr< Scantable > STMath::applyToPol( const CountedPtr<Scantable>& in,
2796 STPol::polOperation fptr,
2797 Float phase )
2798{
2799 CountedPtr< Scantable > out = getScantable(in, false);
2800 Table& tout = out->table();
2801 Block<String> cols(4);
2802 cols[0] = String("SCANNO");
2803 cols[1] = String("BEAMNO");
2804 cols[2] = String("IFNO");
2805 cols[3] = String("CYCLENO");
2806 TableIterator iter(tout, cols);
2807 CountedPtr<STPol> stpol = STPol::getPolClass(out->factories_,
2808 out->getPolType() );
2809 while (!iter.pastEnd()) {
2810 Table t = iter.table();
2811 ArrayColumn<Float> speccol(t, "SPECTRA");
2812 ScalarColumn<uInt> focidcol(t, "FOCUS_ID");
2813 Matrix<Float> pols(speccol.getColumn());
2814 try {
2815 stpol->setSpectra(pols);
2816 Float fang,fhand;
2817 fang = in->focusTable_.getTotalAngle(focidcol(0));
2818 fhand = in->focusTable_.getFeedHand(focidcol(0));
2819 stpol->setPhaseCorrections(fang, fhand);
2820 // use a member function pointer in STPol. This only works on
2821 // the STPol pointer itself, not the Counted Pointer so
2822 // derefernce it.
2823 (&(*(stpol))->*fptr)(phase);
2824 speccol.putColumn(stpol->getSpectra());
2825 } catch (AipsError& e) {
2826 //delete stpol;stpol=0;
2827 throw(e);
2828 }
2829 ++iter;
2830 }
2831 //delete stpol;stpol=0;
2832 return out;
2833}
2834
2835CountedPtr< Scantable >
2836 STMath::swapPolarisations( const CountedPtr< Scantable > & in )
2837{
2838 CountedPtr< Scantable > out = getScantable(in, false);
2839 Table& tout = out->table();
2840 Table t0 = tout(tout.col("POLNO") == 0);
2841 Table t1 = tout(tout.col("POLNO") == 1);
2842 if ( t0.nrow() != t1.nrow() )
2843 throw(AipsError("Inconsistent number of polarisations"));
2844 ArrayColumn<Float> speccol0(t0, "SPECTRA");
2845 ArrayColumn<uChar> flagcol0(t0, "FLAGTRA");
2846 ArrayColumn<Float> speccol1(t1, "SPECTRA");
2847 ArrayColumn<uChar> flagcol1(t1, "FLAGTRA");
2848 Matrix<Float> s0 = speccol0.getColumn();
2849 Matrix<uChar> f0 = flagcol0.getColumn();
2850 speccol0.putColumn(speccol1.getColumn());
2851 flagcol0.putColumn(flagcol1.getColumn());
2852 speccol1.putColumn(s0);
2853 flagcol1.putColumn(f0);
2854 return out;
2855}
2856
2857CountedPtr< Scantable >
2858 STMath::averagePolarisations( const CountedPtr< Scantable > & in,
2859 const std::vector<bool>& mask,
2860 const std::string& weight )
2861{
2862 if (in->npol() < 2 )
2863 throw(AipsError("averagePolarisations can only be applied to two or more"
2864 "polarisations"));
2865 bool insitu = insitu_;
2866 setInsitu(false);
2867 CountedPtr< Scantable > pols = getScantable(in, true);
2868 setInsitu(insitu);
2869 Table& tout = pols->table();
2870 std::string taql = "SELECT FROM $1 WHERE POLNO IN [0,1]";
2871 Table tab = tableCommand(taql, in->table());
2872 if (tab.nrow() == 0 )
2873 throw(AipsError("Could not find any rows with POLNO==0 and POLNO==1"));
2874 TableCopy::copyRows(tout, tab);
2875 TableVector<uInt> vec(tout, "POLNO");
2876 vec = 0;
2877 pols->table_.rwKeywordSet().define("nPol", Int(1));
2878 pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
2879 //pols->table_.rwKeywordSet().define("POLTYPE", in->getPolType());
2880 std::vector<CountedPtr<Scantable> > vpols;
2881 vpols.push_back(pols);
2882 CountedPtr< Scantable > out = average(vpols, mask, weight, "SCAN");
2883 return out;
2884}
2885
2886CountedPtr< Scantable >
2887 STMath::averageBeams( const CountedPtr< Scantable > & in,
2888 const std::vector<bool>& mask,
2889 const std::string& weight )
2890{
2891 bool insitu = insitu_;
2892 setInsitu(false);
2893 CountedPtr< Scantable > beams = getScantable(in, false);
2894 setInsitu(insitu);
2895 Table& tout = beams->table();
2896 // give all rows the same BEAMNO
2897 TableVector<uInt> vec(tout, "BEAMNO");
2898 vec = 0;
2899 beams->table_.rwKeywordSet().define("nBeam", Int(1));
2900 std::vector<CountedPtr<Scantable> > vbeams;
2901 vbeams.push_back(beams);
2902 CountedPtr< Scantable > out = average(vbeams, mask, weight, "SCAN");
2903 return out;
2904}
2905
2906
2907CountedPtr< Scantable >
2908 asap::STMath::frequencyAlign( const CountedPtr< Scantable > & in,
2909 const std::string & refTime,
2910 const std::string & method)
2911{
2912 LogIO os( casa::LogOrigin("STMath", "frequencyAlign()", WHERE));
2913 // clone as this is not working insitu
2914 bool insitu = insitu_;
2915 setInsitu(false);
2916 CountedPtr< Scantable > out = getScantable(in, false);
2917 setInsitu(insitu);
2918 Table& tout = out->table();
2919 // Get reference Epoch to time of first row or given String
2920 Unit DAY(String("d"));
2921 MEpoch::Ref epochRef(in->getTimeReference());
2922 MEpoch refEpoch;
2923 if (refTime.length()>0) {
2924 Quantum<Double> qt;
2925 if (MVTime::read(qt,refTime)) {
2926 MVEpoch mv(qt);
2927 refEpoch = MEpoch(mv, epochRef);
2928 } else {
2929 throw(AipsError("Invalid format for Epoch string"));
2930 }
2931 } else {
2932 refEpoch = in->timeCol_(0);
2933 }
2934 MPosition refPos = in->getAntennaPosition();
2935
2936 InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
2937 /*
2938 // Comment from MV.
2939 // the following code has been commented out because different FREQ_IDs have to be aligned together even
2940 // if the frame doesn't change. So far, lack of this check didn't cause any problems.
2941 // test if user frame is different to base frame
2942 if ( in->frequencies().getFrameString(true)
2943 == in->frequencies().getFrameString(false) ) {
2944 throw(AipsError("Can't convert as no output frame has been set"
2945 " (use set_freqframe) or it is aligned already."));
2946 }
2947 */
2948 MFrequency::Types system = in->frequencies().getFrame();
2949 MVTime mvt(refEpoch.getValue());
2950 String epochout = mvt.string(MVTime::YMD) + String(" (") + refEpoch.getRefString() + String(")");
2951 os << "Aligning at reference Epoch " << epochout
2952 << " in frame " << MFrequency::showType(system) << LogIO::POST;
2953 // set up the iterator
2954 Block<String> cols(4);
2955 // select by constant direction
2956 cols[0] = String("SRCNAME");
2957 cols[1] = String("BEAMNO");
2958 // select by IF ( no of channels varies over this )
2959 cols[2] = String("IFNO");
2960 // select by restfrequency
2961 cols[3] = String("MOLECULE_ID");
2962 TableIterator iter(tout, cols);
2963 while ( !iter.pastEnd() ) {
2964 Table t = iter.table();
2965 ROScalarColumn<String> snCol(t, "SRCNAME");
2966 os << "Aligning to position of source '" << snCol(0) << "'" << LogIO::POST;
2967 MDirection::ROScalarColumn dirCol(t, "DIRECTION");
2968 TableIterator fiter(t, "FREQ_ID");
2969 // determine nchan from the first row. This should work as
2970 // we are iterating over BEAMNO and IFNO
2971 // we should have constant direction
2972
2973 ROArrayColumn<Float> sCol(t, "SPECTRA");
2974 const MDirection direction = dirCol(0);
2975 const uInt nchan = sCol(0).nelements();
2976
2977 // skip operations if there is nothing to align
2978 if (fiter.pastEnd()) {
2979 continue;
2980 }
2981
2982 Table ftab = fiter.table();
2983 // align all frequency ids with respect to the first encountered id
2984 ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
2985 // get the SpectralCoordinate for the freqid, which we are iterating over
2986 SpectralCoordinate sC = \
2987 in->frequencies().getSpectralCoordinate(freqidCol(0));
2988 FrequencyAligner<Float> fa( sC, nchan, refEpoch,
2989 direction, refPos, system );
2990 // realign the SpectralCoordinate and put into the output Scantable
2991 Vector<String> units(1);
2992 units = String("Hz");
2993 Bool linear=True;
2994 SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
2995 sc2.setWorldAxisUnits(units);
2996 const uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
2997 sc2.referenceValue()[0],
2998 sc2.increment()[0]);
2999 while ( !fiter.pastEnd() ) {
3000
3001 ftab = fiter.table();
3002 // spectral coordinate for the current FREQ_ID
3003 ScalarColumn<uInt> freqidCol2(ftab, "FREQ_ID");
3004 sC = in->frequencies().getSpectralCoordinate(freqidCol2(0));
3005 // create the "global" abcissa for alignment with same FREQ_ID
3006 Vector<Double> abc(nchan);
3007 for (uInt i=0; i<nchan; i++) {
3008 Double w;
3009 sC.toWorld(w,Double(i));
3010 abc[i] = w;
3011 }
3012 TableVector<uInt> tvec(ftab, "FREQ_ID");
3013 // assign new frequency id to all rows
3014 tvec = id;
3015 // cache abcissa for same time stamps, so iterate over those
3016 TableIterator timeiter(ftab, "TIME");
3017 while ( !timeiter.pastEnd() ) {
3018 Table tab = timeiter.table();
3019 ArrayColumn<Float> specCol(tab, "SPECTRA");
3020 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
3021 MEpoch::ROScalarColumn timeCol(tab, "TIME");
3022 // use align abcissa cache after the first row
3023 // these rows should be just be POLNO
3024 bool first = true;
3025 for (int i=0; i<int(tab.nrow()); ++i) {
3026 // input values
3027 Vector<uChar> flag = flagCol(i);
3028 Vector<Bool> mask(flag.shape());
3029 Vector<Float> specOut, spec;
3030 spec = specCol(i);
3031 Vector<Bool> maskOut;Vector<uChar> flagOut;
3032 convertArray(mask, flag);
3033 // alignment
3034 Bool ok = fa.align(specOut, maskOut, abc, spec,
3035 mask, timeCol(i), !first,
3036 interp, False);
3037 (void) ok; // unused stop compiler nagging
3038 // back into scantable
3039 flagOut.resize(maskOut.nelements());
3040 convertArray(flagOut, maskOut);
3041 flagCol.put(i, flagOut);
3042 specCol.put(i, specOut);
3043 // start abcissa caching
3044 first = false;
3045 }
3046 // next timestamp
3047 ++timeiter;
3048 }
3049 // next FREQ_ID
3050 ++fiter;
3051 }
3052 // next aligner
3053 ++iter;
3054 }
3055 // set this afterwards to ensure we are doing insitu correctly.
3056 out->frequencies().setFrame(system, true);
3057 return out;
3058}
3059
3060CountedPtr<Scantable>
3061 asap::STMath::convertPolarisation( const CountedPtr<Scantable>& in,
3062 const std::string & newtype )
3063{
3064 if (in->npol() != 2 && in->npol() != 4)
3065 throw(AipsError("Can only convert two or four polarisations."));
3066 if ( in->getPolType() == newtype )
3067 throw(AipsError("No need to convert."));
3068 if ( ! in->selector_.empty() )
3069 throw(AipsError("Can only convert whole scantable. Unset the selection."));
3070 bool insitu = insitu_;
3071 setInsitu(false);
3072 CountedPtr< Scantable > out = getScantable(in, true);
3073 setInsitu(insitu);
3074 Table& tout = out->table();
3075 tout.rwKeywordSet().define("POLTYPE", String(newtype));
3076
3077 Block<String> cols(4);
3078 cols[0] = "SCANNO";
3079 cols[1] = "CYCLENO";
3080 cols[2] = "BEAMNO";
3081 cols[3] = "IFNO";
3082 TableIterator it(in->originalTable_, cols);
3083 String basetype = in->getPolType();
3084 STPol* stpol = STPol::getPolClass(in->factories_, basetype);
3085 try {
3086 while ( !it.pastEnd() ) {
3087 Table tab = it.table();
3088 uInt row = tab.rowNumbers()[0];
3089 stpol->setSpectra(in->getPolMatrix(row));
3090 Float fang,fhand;
3091 fang = in->focusTable_.getTotalAngle(in->mfocusidCol_(row));
3092 fhand = in->focusTable_.getFeedHand(in->mfocusidCol_(row));
3093 stpol->setPhaseCorrections(fang, fhand);
3094 Int npolout = 0;
3095 for (uInt i=0; i<tab.nrow(); ++i) {
3096 Vector<Float> outvec = stpol->getSpectrum(i, newtype);
3097 if ( outvec.nelements() > 0 ) {
3098 tout.addRow();
3099 TableCopy::copyRows(tout, tab, tout.nrow()-1, 0, 1);
3100 ArrayColumn<Float> sCol(tout,"SPECTRA");
3101 ScalarColumn<uInt> pCol(tout,"POLNO");
3102 sCol.put(tout.nrow()-1 ,outvec);
3103 pCol.put(tout.nrow()-1 ,uInt(npolout));
3104 npolout++;
3105 }
3106 }
3107 tout.rwKeywordSet().define("nPol", npolout);
3108 ++it;
3109 }
3110 } catch (AipsError& e) {
3111 delete stpol;
3112 throw(e);
3113 }
3114 delete stpol;
3115 return out;
3116}
3117
3118CountedPtr< Scantable >
3119 asap::STMath::mxExtract( const CountedPtr< Scantable > & in,
3120 const std::string & scantype )
3121{
3122 bool insitu = insitu_;
3123 setInsitu(false);
3124 CountedPtr< Scantable > out = getScantable(in, true);
3125 setInsitu(insitu);
3126 Table& tout = out->table();
3127 std::string taql = "SELECT FROM $1 WHERE BEAMNO != REFBEAMNO";
3128 if (scantype == "on") {
3129 taql = "SELECT FROM $1 WHERE BEAMNO == REFBEAMNO";
3130 }
3131 Table tab = tableCommand(taql, in->table());
3132 TableCopy::copyRows(tout, tab);
3133 if (scantype == "on") {
3134 // re-index SCANNO to 0
3135 TableVector<uInt> vec(tout, "SCANNO");
3136 vec = 0;
3137 }
3138 return out;
3139}
3140
3141std::vector<float>
3142 asap::STMath::fft( const casa::CountedPtr< Scantable > & in,
3143 const std::vector<int>& whichrow,
3144 bool getRealImag )
3145{
3146 std::vector<float> res;
3147 Table tab = in->table();
3148 std::vector<bool> mask;
3149
3150 if (whichrow.size() < 1) { // for all rows (by default)
3151 int nrow = int(tab.nrow());
3152 for (int i = 0; i < nrow; ++i) {
3153 res = in->execFFT(i, mask, getRealImag);
3154 }
3155 } else { // for specified rows
3156 for (uInt i = 0; i < whichrow.size(); ++i) {
3157 res = in->execFFT(i, mask, getRealImag);
3158 }
3159 }
3160
3161 return res;
3162}
3163
3164
3165CountedPtr<Scantable>
3166 asap::STMath::lagFlag( const CountedPtr<Scantable>& in,
3167 double start, double end,
3168 const std::string& mode )
3169{
3170 CountedPtr<Scantable> out = getScantable(in, false);
3171 Table& tout = out->table();
3172 TableIterator iter(tout, "FREQ_ID");
3173 FFTServer<Float,Complex> ffts;
3174
3175 while ( !iter.pastEnd() ) {
3176 Table tab = iter.table();
3177 Double rp,rv,inc;
3178 ROTableRow row(tab);
3179 const TableRecord& rec = row.get(0);
3180 uInt freqid = rec.asuInt("FREQ_ID");
3181 out->frequencies().getEntry(rp, rv, inc, freqid);
3182 ArrayColumn<Float> specCol(tab, "SPECTRA");
3183 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
3184
3185 for (int i=0; i<int(tab.nrow()); ++i) {
3186 Vector<Float> spec = specCol(i);
3187 Vector<uChar> flag = flagCol(i);
3188 std::vector<bool> mask;
3189 for (uInt j = 0; j < flag.nelements(); ++j) {
3190 mask.push_back(!(flag[j]>0));
3191 }
3192 mathutil::doZeroOrderInterpolation(spec, mask);
3193
3194 Vector<Complex> lags;
3195 ffts.fft0(lags, spec);
3196
3197 Int lag0(start+0.5);
3198 Int lag1(end+0.5);
3199 if (mode == "frequency") {
3200 lag0 = Int(spec.nelements()*abs(inc)/(start)+0.5);
3201 lag1 = Int(spec.nelements()*abs(inc)/(end)+0.5);
3202 }
3203 Int lstart = max(0, lag0);
3204 Int lend = min(Int(lags.nelements()-1), lag1);
3205 if (lstart == lend) {
3206 lags[lstart] = Complex(0.0);
3207 } else {
3208 if (lstart > lend) {
3209 Int tmp = lend;
3210 lend = lstart;
3211 lstart = tmp;
3212 }
3213 for (int j=lstart; j <=lend ;++j) {
3214 lags[j] = Complex(0.0);
3215 }
3216 }
3217
3218 ffts.fft0(spec, lags);
3219
3220 specCol.put(i, spec);
3221 }
3222 ++iter;
3223 }
3224 return out;
3225}
3226
3227// Averaging spectra with different channel/resolution
3228CountedPtr<Scantable>
3229STMath::new_average( const std::vector<CountedPtr<Scantable> >& in,
3230 const bool& compel,
3231 const std::vector<bool>& mask,
3232 const std::string& weight,
3233 const std::string& avmode )
3234 throw ( casa::AipsError )
3235{
3236 LogIO os( LogOrigin( "STMath", "new_average()", WHERE ) ) ;
3237 if ( avmode == "SCAN" && in.size() != 1 )
3238 throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
3239 "Use merge first."));
3240
3241 CountedPtr<Scantable> out ; // processed result
3242 if ( compel ) {
3243 std::vector< CountedPtr<Scantable> > newin ; // input for average process
3244 uInt insize = in.size() ; // number of input scantables
3245
3246 // setup newin
3247 bool oldInsitu = insitu_ ;
3248 setInsitu( false ) ;
3249 newin.resize( insize ) ;
3250 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3251 newin[itable] = getScantable( in[itable], false ) ;
3252 }
3253 setInsitu( oldInsitu ) ;
3254
3255 // warning
3256 os << "Average spectra with different spectral resolution" << LogIO::POST ;
3257
3258 // temporarily set coordinfo
3259 vector<string> oldinfo( insize ) ;
3260 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3261 vector<string> coordinfo = in[itable]->getCoordInfo() ;
3262 oldinfo[itable] = coordinfo[0] ;
3263 coordinfo[0] = "Hz" ;
3264 newin[itable]->setCoordInfo( coordinfo ) ;
3265 }
3266
3267 ostringstream oss ;
3268
3269 // check IF frequency coverage
3270 // freqid: list of FREQ_ID, which is used, in each table
3271 // iffreq: list of minimum and maximum frequency for each FREQ_ID in
3272 // each table
3273 // freqid[insize][numIF]
3274 // freqid: [[id00, id01, ...],
3275 // [id10, id11, ...],
3276 // ...
3277 // [idn0, idn1, ...]]
3278 // iffreq[insize][numIF*2]
3279 // iffreq: [[min_id00, max_id00, min_id01, max_id01, ...],
3280 // [min_id10, max_id10, min_id11, max_id11, ...],
3281 // ...
3282 // [min_idn0, max_idn0, min_idn1, max_idn1, ...]]
3283 //os << "Check IF settings in each table" << LogIO::POST ;
3284 vector< vector<uInt> > freqid( insize );
3285 vector< vector<double> > iffreq( insize ) ;
3286 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3287 Vector<uInt> freqIds = newin[itable]->mfreqidCol_.getColumn() ;
3288 vector<uInt> uniqueFreqId = newin[itable]->getNumbers(newin[itable]->mfreqidCol_) ;
3289 for ( vector<uInt>::iterator i = uniqueFreqId.begin() ;
3290 i != uniqueFreqId.end() ; i++ ) {
3291 //os << "itable = " << itable << ": IF " << id << " is included in the list" << LogIO::POST ;
3292 uInt target = 0 ;
3293 while ( freqIds[target] != *i )
3294 target++ ;
3295 vector<double> abcissa = newin[itable]->getAbcissa( target ) ;
3296 freqid[itable].push_back( *i ) ;
3297 double incr = abs( abcissa[1] - abcissa[0] ) ;
3298 iffreq[itable].push_back( (*min_element(abcissa.begin(),abcissa.end()))-0.5*incr ) ;
3299 iffreq[itable].push_back( (*max_element(abcissa.begin(),abcissa.end()))+0.5*incr ) ;
3300 }
3301 }
3302
3303 // debug
3304// os << "IF settings summary:" << endl ;
3305// for ( uInt i = 0 ; i < freqid.size() ; i++ ) {
3306// os << " Table" << i << endl ;
3307// for ( uInt j = 0 ; j < freqid[i].size() ; j++ ) {
3308// os << " id = " << freqid[i][j] << " (min,max) = (" << iffreq[i][2*j] << "," << iffreq[i][2*j+1] << ")" << endl ;
3309// }
3310// }
3311// os << endl ;
3312// os.post() ;
3313
3314 // IF grouping based on their frequency coverage
3315 // ifgrp: number of member in each IF group
3316 // ifgrp[numgrp]
3317 // ifgrp: [n0, n1, ...]
3318 //os << "IF grouping based on their frequency coverage" << LogIO::POST ;
3319
3320 // parameter for IF grouping
3321 // groupmode = OR retrieve all region
3322 // AND only retrieve overlaped region
3323 //string groupmode = "AND" ;
3324 string groupmode = "OR" ;
3325 uInt sizecr = 0 ;
3326 if ( groupmode == "AND" )
3327 sizecr = 1 ;
3328 else if ( groupmode == "OR" )
3329 sizecr = 0 ;
3330
3331 vector<double> sortedfreq ;
3332 for ( uInt i = 0 ; i < iffreq.size() ; i++ ) {
3333 for ( uInt j = 0 ; j < iffreq[i].size() ; j++ ) {
3334 if ( count( sortedfreq.begin(), sortedfreq.end(), iffreq[i][j] ) == 0 )
3335 sortedfreq.push_back( iffreq[i][j] ) ;
3336 }
3337 }
3338 sort( sortedfreq.begin(), sortedfreq.end() ) ;
3339 vector<uInt> ifgrp( sortedfreq.size()-1, 0 ) ;
3340 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3341 for ( uInt iif = 0 ; iif < freqid[itable].size() ; iif++ ) {
3342 double range0 = iffreq[itable][2*iif] ;
3343 double range1 = iffreq[itable][2*iif+1] ;
3344 for ( uInt j = 0 ; j < sortedfreq.size()-1 ; j++ ) {
3345 double fmin = max( range0, sortedfreq[j] ) ;
3346 double fmax = min( range1, sortedfreq[j+1] ) ;
3347 if ( fmin < fmax ) {
3348 ifgrp[j]++ ;
3349 }
3350 }
3351 }
3352 }
3353
3354 // Grouping continuous IF groups (without frequency gap)
3355 // freqgrp: list of IF group indexes in each frequency group
3356 // freqgrp[numgrp][nummember]
3357 // freqgrp: [[ifgrp00, ifgrp01, ifgrp02, ...],
3358 // [ifgrp10, ifgrp11, ifgrp12, ...],
3359 // ...
3360 // [ifgrpn0, ifgrpn1, ifgrpn2, ...]]
3361 // grprange[2*numgrp]
3362 // grprange: [fmin0,fmax0,fmin1,fmax1,...]
3363 vector< vector<uInt> > freqgrp ;
3364 vector<double> grprange ;
3365 vector<uInt> grpedge ;
3366 for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
3367 if ( ifgrp[igrp] <= sizecr ) {
3368 grpedge.push_back( igrp ) ;
3369 }
3370 }
3371 grpedge.push_back( ifgrp.size() ) ;
3372 uInt itmp = 0 ;
3373 for ( uInt i = 0 ; i < grpedge.size() ; i++ ) {
3374 int n = grpedge[i] - itmp ;
3375 if ( n > 0 ) {
3376 vector<uInt> members( n ) ;
3377 for ( int j = 0 ; j < n ; j++ ) {
3378 members[j] = itmp+j ;
3379 }
3380 freqgrp.push_back( members ) ;
3381 grprange.push_back( sortedfreq[itmp] ) ;
3382 grprange.push_back( sortedfreq[grpedge[i]] ) ;
3383 }
3384 itmp += n + 1 ;
3385 }
3386
3387 // print frequency group
3388 oss.str("") ;
3389 oss << "Frequency Group summary: " << endl ;
3390 oss << " GROUP_ID: [FREQ_MIN, FREQ_MAX]" << endl ;
3391 for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
3392 oss << " GROUP " << setw( 2 ) << i << ": [" << grprange[2*i] << "," << grprange[2*i+1] << "]" ;
3393 oss << endl ;
3394 }
3395 oss << endl ;
3396 os << oss.str() << LogIO::POST ;
3397
3398 // groups: list of frequency group index whose frequency range overlaps
3399 // with that of each table and IF
3400 // groups[numtable][numIF]
3401 // groups: [[grpx, grpy,...],
3402 // [grpa, grpb,...],
3403 // ...
3404 // [grpk, grpm,...]]
3405 vector< vector<uInt> > groups( insize ) ;
3406 for ( uInt i = 0 ; i < insize ; i++ ) {
3407 groups[i].resize( freqid[i].size() ) ;
3408 }
3409 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3410 for ( uInt ifreq = 0 ; ifreq < freqid[itable].size() ; ifreq++ ) {
3411 double minf = iffreq[itable][2*ifreq] ;
3412 uInt groupid ;
3413 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3414 vector<uInt> memberlist = freqgrp[igrp] ;
3415 if ( (minf >= grprange[2*igrp]) && (minf <= grprange[2*igrp+1]) ) {
3416 groupid = igrp ;
3417 break ;
3418 }
3419 }
3420 groups[itable][ifreq] = groupid ;
3421 }
3422 }
3423
3424
3425 // print membership
3426 oss.str("") ;
3427 for ( uInt i = 0 ; i < insize ; i++ ) {
3428 oss << "Table " << i << endl ;
3429 for ( uInt j = 0 ; j < groups[i].size() ; j++ ) {
3430 oss << " FREQ_ID " << setw( 2 ) << freqid[i][j] << ": " ;
3431 oss << setw( 2 ) << groups[i][j] ;
3432 oss << endl ;
3433 }
3434 }
3435 os << oss.str() << LogIO::POST ;
3436
3437 // reset SCANNO and IFNO/FREQ_ID: IF is reset by the result of sortation
3438 //os << "All IF number is set to IF group index" << LogIO::POST ;
3439 // reset SCANNO only when avmode != "SCAN"
3440 if ( avmode != "SCAN" ) {
3441 os << "All scan number is set to 0" << LogIO::POST ;
3442 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3443 uInt nrow = newin[itable]->nrow() ;
3444 Vector<uInt> resetScan( nrow, 0 ) ;
3445 newin[itable]->scanCol_.putColumn( resetScan ) ;
3446 }
3447 }
3448
3449 // reset spectra and flagtra: align spectral resolution
3450 //os << "Align spectral resolution" << LogIO::POST ;
3451 // gmaxdnu: the coarsest frequency resolution in the frequency group
3452 // gminfreq: lower frequency edge of the frequency group
3453 // gnchan: number of channels for the frequency group
3454 vector<double> gmaxdnu( freqgrp.size(), 0.0 ) ;
3455 vector<double> gminfreq( freqgrp.size() ) ;
3456 vector<double> gnchan( freqgrp.size() ) ;
3457 for ( uInt i = 0 ; i < insize ; i++ ) {
3458 vector<uInt> members = groups[i] ;
3459 for ( uInt j = 0 ; j < members.size() ; j++ ) {
3460 uInt groupid = members[j] ;
3461 Double rp,rv,ic ;
3462 newin[i]->frequencies().getEntry( rp, rv, ic, j ) ;
3463 if ( abs(ic) > abs(gmaxdnu[groupid]) )
3464 gmaxdnu[groupid] = ic ;
3465 }
3466 }
3467 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3468 gminfreq[igrp] = grprange[2*igrp] ;
3469 double maxfreq = grprange[2*igrp+1] ;
3470 gnchan[igrp] = (int)(abs((maxfreq-gminfreq[igrp])/gmaxdnu[igrp])+0.9) ;
3471 }
3472
3473 // regrid spectral data and update frequency info
3474 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3475 Vector<uInt> oldFreqId = newin[itable]->mfreqidCol_.getColumn() ;
3476 Vector<uInt> newFreqId( oldFreqId.nelements() ) ;
3477
3478 // update MAIN
3479 for ( uInt irow = 0 ; irow < newin[itable]->nrow() ; irow++ ) {
3480 uInt groupid = groups[itable][oldFreqId[irow]] ;
3481 newFreqId[irow] = groupid ;
3482 newin[itable]->regridChannel( gnchan[groupid],
3483 gmaxdnu[groupid],
3484 gminfreq[groupid],
3485 irow ) ;
3486 }
3487 newin[itable]->mfreqidCol_.putColumn( newFreqId ) ;
3488 newin[itable]->ifCol_.putColumn( newFreqId ) ;
3489
3490 // update FREQUENCIES
3491 Table tab = newin[itable]->frequencies().table() ;
3492 ScalarColumn<uInt> fIdCol( tab, "ID" ) ;
3493 ScalarColumn<Double> fRefPixCol( tab, "REFPIX" ) ;
3494 ScalarColumn<Double> fRefValCol( tab, "REFVAL" ) ;
3495 ScalarColumn<Double> fIncrCol( tab, "INCREMENT" ) ;
3496 if ( freqgrp.size() > tab.nrow() ) {
3497 tab.addRow( freqgrp.size()-tab.nrow() ) ;
3498 }
3499 for ( uInt irow = 0 ; irow < freqgrp.size() ; irow++ ) {
3500 Double refval = gminfreq[irow] + 0.5 * abs(gmaxdnu[irow]) ;
3501 Double refpix = (gmaxdnu[irow] > 0.0) ? 0 : gnchan[irow]-1 ;
3502 Double increment = gmaxdnu[irow] ;
3503 fIdCol.put( irow, irow ) ;
3504 fRefPixCol.put( irow, refpix ) ;
3505 fRefValCol.put( irow, refval ) ;
3506 fIncrCol.put( irow, increment ) ;
3507 }
3508 }
3509
3510 // set back coordinfo
3511 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3512 vector<string> coordinfo = newin[itable]->getCoordInfo() ;
3513 coordinfo[0] = oldinfo[itable] ;
3514 newin[itable]->setCoordInfo( coordinfo ) ;
3515 }
3516
3517 // average
3518 out = average( newin, mask, weight, avmode ) ;
3519 }
3520 else {
3521 // simple average
3522 out = average( in, mask, weight, avmode ) ;
3523 }
3524
3525 return out;
3526}
3527
3528CountedPtr<Scantable> STMath::cwcal( const CountedPtr<Scantable>& s,
3529 const String calmode,
3530 const String antname )
3531{
3532 // frequency switch
3533 if ( calmode == "fs" ) {
3534 return cwcalfs( s, antname ) ;
3535 }
3536 else {
3537 vector<bool> masks = s->getMask( 0 ) ;
3538 vector<int> types ;
3539
3540 // save original table selection
3541 Table torg = s->table_ ;
3542
3543 // sky scan
3544 bool insitu = insitu_ ;
3545 insitu_ = false ;
3546 // share calibration scans before average with out
3547 CountedPtr<Scantable> out = getScantable( s, true ) ;
3548 insitu_ = insitu ;
3549 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::SKY ) ;
3550 out->attach() ;
3551 CountedPtr<Scantable> asky = averageWithinSession( out,
3552 masks,
3553 "TINT" ) ;
3554 // hot scan
3555 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::HOT ) ;
3556 out->attach() ;
3557 CountedPtr<Scantable> ahot = averageWithinSession( out,
3558 masks,
3559 "TINT" ) ;
3560 // cold scan
3561 CountedPtr<Scantable> acold ;
3562// out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::COLD ) ;
3563// out->attach() ;
3564// CountedPtr<Scantable> acold = averageWithinSession( out,
3565// masks,
3566// "TINT" ) ;
3567
3568 // off scan
3569 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSOFF ) ;
3570 out->attach() ;
3571 CountedPtr<Scantable> aoff = averageWithinSession( out,
3572 masks,
3573 "TINT" ) ;
3574
3575 // on scan
3576 s->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSON ) ;
3577 s->attach() ;
3578 out->table_ = out->originalTable_ ;
3579 out->attach() ;
3580 out->table().addRow( s->nrow() ) ;
3581 copyRows( out->table(), s->table(), 0, 0, s->nrow(), False, True, False ) ;
3582
3583 // process each on scan
3584 STSelector sel ;
3585 vector<string> cols( 3 ) ;
3586 cols[0] = "BEAMNO" ;
3587 cols[1] = "POLNO" ;
3588 cols[2] = "IFNO" ;
3589 STIdxIter2 iter(out, cols);
3590 while ( !iter.pastEnd() ) {
3591 Record ids = iter.currentValue();
3592 stringstream ss ;
3593 ss << "SELECT FROM $1 WHERE "
3594 << "BEAMNO==" << ids.asuInt("BEAMNO") << "&&"
3595 << "POLNO==" << ids.asuInt("POLNO") << "&&"
3596 << "IFNO==" << ids.asuInt("IFNO") ;
3597 //cout << "TaQL string: " << ss.str() << endl ;
3598 sel.setTaQL( ss.str() ) ;
3599 aoff->setSelection( sel ) ;
3600 ahot->setSelection( sel ) ;
3601 asky->setSelection( sel ) ;
3602 Vector<uInt> rows = iter.getRows( SHARE ) ;
3603 // out should be an exact copy of s except that SPECTRA column is empty
3604 calibrateCW( out, s, aoff, asky, ahot, acold, rows, antname ) ;
3605 aoff->unsetSelection() ;
3606 ahot->unsetSelection() ;
3607 asky->unsetSelection() ;
3608 sel.reset() ;
3609 iter.next() ;
3610 }
3611 s->table_ = torg ;
3612 s->attach() ;
3613
3614 // flux unit
3615 out->setFluxUnit( "K" ) ;
3616
3617 return out ;
3618 }
3619}
3620
3621CountedPtr<Scantable> STMath::almacal( const CountedPtr<Scantable>& s,
3622 const String calmode )
3623{
3624 // frequency switch
3625 if ( calmode == "fs" ) {
3626 return almacalfs( s ) ;
3627 }
3628 else {
3629// double t0, t1 ;
3630// t0 = mathutil::gettimeofday_sec() ;
3631 vector<bool> masks = s->getMask( 0 ) ;
3632
3633 // save original table selection
3634 Table torg = s->table_ ;
3635
3636 // off scan
3637 // TODO 2010/01/08 TN
3638 // Grouping by time should be needed before averaging.
3639 // Each group must have own unique SCANNO (should be renumbered).
3640 // See PIPELINE/SDCalibration.py
3641 bool insitu = insitu_ ;
3642 insitu_ = false ;
3643 // share off scan before average with out
3644 CountedPtr<Scantable> out = getScantable( s, true ) ;
3645 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSOFF ) ;
3646 out->attach() ;
3647 insitu_ = insitu ;
3648 CountedPtr<Scantable> aoff = averageWithinSession( out,
3649 masks,
3650 "TINT" ) ;
3651
3652 // on scan
3653// t0 = mathutil::gettimeofday_sec() ;
3654 s->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSON ) ;
3655 s->attach() ;
3656 out->table_ = out->originalTable_ ;
3657 out->attach() ;
3658 out->table().addRow( s->nrow() ) ;
3659 copyRows( out->table(), s->table(), 0, 0, s->nrow(), False ) ;
3660// t1 = mathutil::gettimeofday_sec() ;
3661// cout << "elapsed time for preparing output table: " << t1-t0 << " sec" << endl ;
3662
3663 // process each on scan
3664// t0 = mathutil::gettimeofday_sec() ;
3665
3666 // using STIdxIter2
3667 vector<string> cols( 3 ) ;
3668 cols[0] = "BEAMNO" ;
3669 cols[1] = "POLNO" ;
3670 cols[2] = "IFNO" ;
3671 STIdxIter2 iter( out, cols ) ;
3672 STSelector sel ;
3673 while ( !iter.pastEnd() ) {
3674 Record ids = iter.currentValue() ;
3675 stringstream ss ;
3676 ss << "SELECT FROM $1 WHERE "
3677 << "BEAMNO==" << ids.asuInt(cols[0]) << "&&"
3678 << "POLNO==" << ids.asuInt(cols[1]) << "&&"
3679 << "IFNO==" << ids.asuInt(cols[2]) ;
3680 //cout << "TaQL string: " << ss.str() << endl ;
3681 sel.setTaQL( ss.str() ) ;
3682 aoff->setSelection( sel ) ;
3683 Vector<uInt> rows = iter.getRows( SHARE ) ;
3684 // out should be an exact copy of s except that SPECTRA column is empty
3685 calibrateALMA( out, s, aoff, rows ) ;
3686 aoff->unsetSelection() ;
3687 sel.reset() ;
3688 iter.next() ;
3689 }
3690 s->table_ = torg ;
3691 s->attach() ;
3692
3693// t1 = mathutil::gettimeofday_sec() ;
3694// cout << "elapsed time for calibration: " << t1-t0 << " sec" << endl ;
3695
3696 // flux unit
3697 out->setFluxUnit( "K" ) ;
3698
3699 return out ;
3700 }
3701}
3702
3703CountedPtr<Scantable> STMath::cwcalfs( const CountedPtr<Scantable>& s,
3704 const String antname )
3705{
3706 vector<int> types ;
3707
3708 // APEX calibration mode
3709 int apexcalmode = 1 ;
3710
3711 if ( antname.find( "APEX" ) != string::npos ) {
3712 // check if off scan exists or not
3713 STSelector sel = STSelector() ;
3714 //sel.setName( offstr1 ) ;
3715 types.push_back( SrcType::FLOOFF ) ;
3716 sel.setTypes( types ) ;
3717 try {
3718 s->setSelection( sel ) ;
3719 }
3720 catch ( AipsError &e ) {
3721 apexcalmode = 0 ;
3722 }
3723 sel.reset() ;
3724 }
3725 s->unsetSelection() ;
3726 types.clear() ;
3727
3728 vector<bool> masks = s->getMask( 0 ) ;
3729 CountedPtr<Scantable> ssig, sref ;
3730 //CountedPtr<Scantable> out ;
3731 bool insitu = insitu_ ;
3732 insitu_ = False ;
3733 CountedPtr<Scantable> out = getScantable( s, true ) ;
3734 insitu_ = insitu ;
3735
3736 if ( antname.find( "APEX" ) != string::npos ) {
3737 // APEX calibration
3738 // sky scan
3739 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOSKY ) ;
3740 out->attach() ;
3741 CountedPtr<Scantable> askylo = averageWithinSession( out,
3742 masks,
3743 "TINT" ) ;
3744 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHISKY ) ;
3745 out->attach() ;
3746 CountedPtr<Scantable> askyhi = averageWithinSession( out,
3747 masks,
3748 "TINT" ) ;
3749
3750 // hot scan
3751 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOHOT ) ;
3752 out->attach() ;
3753 CountedPtr<Scantable> ahotlo = averageWithinSession( out,
3754 masks,
3755 "TINT" ) ;
3756 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHIHOT ) ;
3757 out->attach() ;
3758 CountedPtr<Scantable> ahothi = averageWithinSession( out,
3759 masks,
3760 "TINT" ) ;
3761
3762 // cold scan
3763 CountedPtr<Scantable> acoldlo, acoldhi ;
3764// out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOCOLD ) ;
3765// out->attach() ;
3766// CountedPtr<Scantable> acoldlo = averageWithinSession( out,
3767// masks,
3768// "TINT" ) ;
3769// out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHICOLD ) ;
3770// out->attach() ;
3771// CountedPtr<Scantable> acoldhi = averageWithinSession( out,
3772// masks,
3773// "TINT" ) ;
3774
3775 // ref scan
3776 insitu_ = false ;
3777 sref = getScantable( s, true ) ;
3778 CountedPtr<Scantable> rref = getScantable( s, true ) ;
3779 insitu_ = insitu ;
3780 rref->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FSLO ) ;
3781 rref->attach() ;
3782 copyRows( sref->table_, rref->table_, 0, 0, rref->nrow(), False, True, False ) ;
3783
3784 // sig scan
3785 insitu_ = false ;
3786 ssig = getScantable( s, true ) ;
3787 CountedPtr<Scantable> rsig = getScantable( s, true ) ;
3788 insitu_ = insitu ;
3789 rsig->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FSHI ) ;
3790 rsig->attach() ;
3791 copyRows( ssig->table_, rsig->table_, 0, 0, rsig->nrow(), False, True, False ) ;
3792
3793 if ( apexcalmode == 0 ) {
3794 // using STIdxIter2
3795 vector<string> cols( 3 ) ;
3796 cols[0] = "BEAMNO" ;
3797 cols[1] = "POLNO" ;
3798 cols[2] = "IFNO" ;
3799 STIdxIter2 iter(ssig, cols) ;
3800 STSelector sel ;
3801 vector< CountedPtr<Scantable> > on( 2 ) ;
3802 on[0] = rsig ;
3803 on[1] = rref ;
3804 vector< CountedPtr<Scantable> > sky( 2 ) ;
3805 sky[0] = askylo ;
3806 sky[1] = askyhi ;
3807 vector< CountedPtr<Scantable> > hot( 2 ) ;
3808 hot[0] = ahotlo ;
3809 hot[1] = ahothi ;
3810 vector< CountedPtr<Scantable> > cold( 2 ) ;
3811 while ( !iter.pastEnd() ) {
3812 Record ids = iter.currentValue() ;
3813 stringstream ss ;
3814 ss << "SELECT FROM $1 WHERE "
3815 << "BEAMNO==" << ids.asuInt("BEAMNO") << "&&"
3816 << "POLNO==" << ids.asuInt("POLNO") << "&&"
3817 << "IFNO==" << ids.asuInt("IFNO") ;
3818 //cout << "TaQL string: " << ss.str() << endl ;
3819 sel.setTaQL( ss.str() ) ;
3820 sky[0]->setSelection( sel ) ;
3821 sky[1]->setSelection( sel ) ;
3822 hot[0]->setSelection( sel ) ;
3823 hot[1]->setSelection( sel ) ;
3824 Vector<uInt> rows = iter.getRows( SHARE ) ;
3825 calibrateAPEXFS( ssig, sref, on, sky, hot, cold, rows ) ;
3826 sky[0]->unsetSelection() ;
3827 sky[1]->unsetSelection() ;
3828 hot[0]->unsetSelection() ;
3829 hot[1]->unsetSelection() ;
3830 sel.reset() ;
3831 iter.next() ;
3832 }
3833
3834 }
3835 else if ( apexcalmode == 1 ) {
3836 // APEX fs data with off scan
3837 // off scan
3838 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOOFF ) ;
3839 out->attach() ;
3840 CountedPtr<Scantable> aofflo = averageWithinSession( out,
3841 masks,
3842 "TINT" ) ;
3843 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHIOFF ) ;
3844 out->attach() ;
3845 CountedPtr<Scantable> aoffhi = averageWithinSession( out,
3846 masks,
3847 "TINT" ) ;
3848
3849 // process each sig and ref scan
3850// STSelector sel ;
3851 vector<string> cols( 3 ) ;
3852 cols[0] = "BEAMNO" ;
3853 cols[1] = "POLNO" ;
3854 cols[2] = "IFNO" ;
3855 STIdxIter2 iter(ssig, cols);
3856 STSelector sel ;
3857 while ( !iter.pastEnd() ) {
3858 Record ids = iter.currentValue() ;
3859 stringstream ss ;
3860 ss << "SELECT FROM $1 WHERE "
3861 << "BEAMNO==" << ids.asuInt("BEAMNO") << "&&"
3862 << "POLNO==" << ids.asuInt("POLNO") << "&&"
3863 << "IFNO==" << ids.asuInt("IFNO") ;
3864 //cout << "TaQL string: " << ss.str() << endl ;
3865 sel.setTaQL( ss.str() ) ;
3866 aofflo->setSelection( sel ) ;
3867 ahotlo->setSelection( sel ) ;
3868 askylo->setSelection( sel ) ;
3869 Vector<uInt> rows = iter.getRows( SHARE ) ;
3870 calibrateCW( ssig, rsig, aofflo, askylo, ahotlo, acoldlo, rows, antname ) ;
3871 aofflo->unsetSelection() ;
3872 ahotlo->unsetSelection() ;
3873 askylo->unsetSelection() ;
3874 sel.reset() ;
3875 iter.next() ;
3876 }
3877 iter = STIdxIter2(sref, cols);
3878 while ( !iter.pastEnd() ) {
3879 Record ids = iter.currentValue() ;
3880 stringstream ss ;
3881 ss << "SELECT FROM $1 WHERE "
3882 << "BEAMNO==" << ids.asuInt("BEAMNO") << "&&"
3883 << "POLNO==" << ids.asuInt("POLNO") << "&&"
3884 << "IFNO==" << ids.asuInt("IFNO") ;
3885 //cout << "TaQL string: " << ss.str() << endl ;
3886 sel.setTaQL( ss.str() ) ;
3887 aoffhi->setSelection( sel ) ;
3888 ahothi->setSelection( sel ) ;
3889 askyhi->setSelection( sel ) ;
3890 Vector<uInt> rows = iter.getRows( SHARE ) ;
3891 calibrateCW( sref, rref, aoffhi, askyhi, ahothi, acoldhi, rows, antname ) ;
3892 aoffhi->unsetSelection() ;
3893 ahothi->unsetSelection() ;
3894 askyhi->unsetSelection() ;
3895 sel.reset() ;
3896 iter.next() ;
3897 }
3898 }
3899 }
3900 else {
3901 // non-APEX fs data
3902 // sky scan
3903 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::SKY ) ;
3904 out->attach() ;
3905 CountedPtr<Scantable> asky = averageWithinSession( out,
3906 masks,
3907 "TINT" ) ;
3908 STSelector sel = STSelector() ;
3909
3910 // hot scan
3911 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::HOT ) ;
3912 out->attach() ;
3913 CountedPtr<Scantable> ahot = averageWithinSession( out,
3914 masks,
3915 "TINT" ) ;
3916
3917 // cold scan
3918 CountedPtr<Scantable> acold ;
3919// out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::COLD ) ;
3920// out->attach() ;
3921// CountedPtr<Scantable> acold = averageWithinSession( out,
3922// masks,
3923// "TINT" ) ;
3924
3925 // ref scan
3926 bool insitu = insitu_ ;
3927 insitu_ = false ;
3928 sref = getScantable( s, true ) ;
3929 CountedPtr<Scantable> rref = getScantable( s, true ) ;
3930 insitu_ = insitu ;
3931 rref->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSOFF ) ;
3932 rref->attach() ;
3933 copyRows( sref->table_, rref->table_, 0, 0, rref->nrow(), False, True, False ) ;
3934
3935 // sig scan
3936 insitu_ = false ;
3937 ssig = getScantable( s, true ) ;
3938 CountedPtr<Scantable> rsig = getScantable( s, true ) ;
3939 insitu_ = insitu ;
3940 rsig->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSON ) ;
3941 rsig->attach() ;
3942 copyRows( ssig->table_, rsig->table_, 0, 0, rsig->nrow(), False, True, False ) ;
3943
3944 // process each sig and ref scan
3945 vector<string> cols( 3 ) ;
3946 cols[0] = "BEAMNO" ;
3947 cols[1] = "POLNO" ;
3948 cols[2] = "IFNO" ;
3949 STIdxIter2 iter(ssig, cols);
3950 while ( !iter.pastEnd() ) {
3951 Record ids = iter.currentValue() ;
3952 stringstream ss ;
3953 ss << "SELECT FROM $1 WHERE "
3954 << "BEAMNO==" << ids.asuInt("BEAMNO") << "&&"
3955 << "POLNO==" << ids.asuInt("POLNO") << "&&"
3956 << "IFNO==" << ids.asuInt("IFNO") ;
3957 //cout << "TaQL string: " << ss.str() << endl ;
3958 sel.setTaQL( ss.str() ) ;
3959 ahot->setSelection( sel ) ;
3960 asky->setSelection( sel ) ;
3961 Vector<uInt> rows = iter.getRows( SHARE ) ;
3962 // out should be an exact copy of s except that SPECTRA column is empty
3963 calibrateFS( ssig, sref, rsig, rref, asky, ahot, acold, rows ) ;
3964 ahot->unsetSelection() ;
3965 asky->unsetSelection() ;
3966 sel.reset() ;
3967 iter.next() ;
3968 }
3969 }
3970
3971 // do folding if necessary
3972 Table sigtab = ssig->table() ;
3973 Table reftab = sref->table() ;
3974 ScalarColumn<uInt> reffidCol ;
3975 Int nchan = (Int)ssig->nchan() ;
3976 reffidCol.attach( reftab, "FREQ_ID" ) ;
3977 Vector<uInt> sfids = ssig->mfreqidCol_.getColumn() ;
3978 Vector<uInt> rfids = sref->mfreqidCol_.getColumn() ;
3979 vector<uInt> sfids_unique ;
3980 vector<uInt> rfids_unique ;
3981 vector<uInt> sifno_unique ;
3982 vector<uInt> rifno_unique ;
3983 for ( uInt i = 0 ; i < sfids.nelements() ; i++ ) {
3984 if ( count( sfids_unique.begin(), sfids_unique.end(), sfids[i] ) == 0 ) {
3985 sfids_unique.push_back( sfids[i] ) ;
3986 sifno_unique.push_back( ssig->getIF( i ) ) ;
3987 }
3988 if ( count( rfids_unique.begin(), rfids_unique.end(), rfids[i] ) == 0 ) {
3989 rfids_unique.push_back( rfids[i] ) ;
3990 rifno_unique.push_back( sref->getIF( i ) ) ;
3991 }
3992 }
3993 double refpix_sig, refval_sig, increment_sig ;
3994 double refpix_ref, refval_ref, increment_ref ;
3995 vector< CountedPtr<Scantable> > tmp( sfids_unique.size() ) ;
3996 for ( uInt i = 0 ; i < sfids_unique.size() ; i++ ) {
3997 ssig->frequencies().getEntry( refpix_sig, refval_sig, increment_sig, sfids_unique[i] ) ;
3998 sref->frequencies().getEntry( refpix_ref, refval_ref, increment_ref, rfids_unique[i] ) ;
3999 if ( refpix_sig == refpix_ref ) {
4000 double foffset = refval_ref - refval_sig ;
4001 int choffset = static_cast<int>(foffset/increment_sig) ;
4002 double doffset = foffset / increment_sig ;
4003 if ( abs(choffset) >= nchan ) {
4004 LogIO os( LogOrigin( "STMath", "cwcalfs", WHERE ) ) ;
4005 os << "FREQ_ID=[" << sfids_unique[i] << "," << rfids_unique[i] << "]: out-band frequency switching, no folding" << LogIO::POST ;
4006 os << "Just return signal data" << LogIO::POST ;
4007 //std::vector< CountedPtr<Scantable> > tabs ;
4008 //tabs.push_back( ssig ) ;
4009 //tabs.push_back( sref ) ;
4010 //out = merge( tabs ) ;
4011 tmp[i] = ssig ;
4012 }
4013 else {
4014 STSelector sel = STSelector() ;
4015 vector<int> v( 1, sifno_unique[i] ) ;
4016 sel.setIFs( v ) ;
4017 ssig->setSelection( sel ) ;
4018 sel.reset() ;
4019 v[0] = rifno_unique[i] ;
4020 sel.setIFs( v ) ;
4021 sref->setSelection( sel ) ;
4022 sel.reset() ;
4023 if ( antname.find( "APEX" ) != string::npos ) {
4024 tmp[i] = dofold( ssig, sref, 0.5*doffset, -0.5*doffset ) ;
4025 //tmp[i] = dofold( ssig, sref, doffset ) ;
4026 }
4027 else {
4028 tmp[i] = dofold( ssig, sref, doffset ) ;
4029 }
4030 ssig->unsetSelection() ;
4031 sref->unsetSelection() ;
4032 }
4033 }
4034 }
4035
4036 if ( tmp.size() > 1 ) {
4037 out = merge( tmp ) ;
4038 }
4039 else {
4040 out = tmp[0] ;
4041 }
4042
4043 // flux unit
4044 out->setFluxUnit( "K" ) ;
4045
4046 return out ;
4047}
4048
4049CountedPtr<Scantable> STMath::almacalfs( const CountedPtr<Scantable>& s )
4050{
4051 (void) s; //currently unused
4052 CountedPtr<Scantable> out ;
4053
4054 return out ;
4055}
4056
4057void STMath::calibrateCW( CountedPtr<Scantable> &out,
4058 const CountedPtr<Scantable>& on,
4059 const CountedPtr<Scantable>& off,
4060 const CountedPtr<Scantable>& sky,
4061 const CountedPtr<Scantable>& hot,
4062 const CountedPtr<Scantable>& cold,
4063 const Vector<uInt> &rows,
4064 const String &antname )
4065{
4066 // 2012/05/22 TN
4067 // Assume that out has empty SPECTRA column
4068
4069 // if rows is empty, just return
4070 if ( rows.nelements() == 0 )
4071 return ;
4072 ROScalarColumn<Double> timeCol( off->table(), "TIME" ) ;
4073 Vector<Double> timeOff = timeCol.getColumn() ;
4074 timeCol.attach( sky->table(), "TIME" ) ;
4075 Vector<Double> timeSky = timeCol.getColumn() ;
4076 timeCol.attach( hot->table(), "TIME" ) ;
4077 Vector<Double> timeHot = timeCol.getColumn() ;
4078 timeCol.attach( on->table(), "TIME" ) ;
4079 ROArrayColumn<Float> arrayFloatCol( off->table(), "SPECTRA" ) ;
4080 SpectralData offspectra(arrayFloatCol.getColumn());
4081 arrayFloatCol.attach( sky->table(), "SPECTRA" ) ;
4082 SpectralData skyspectra(arrayFloatCol.getColumn());
4083 arrayFloatCol.attach( hot->table(), "SPECTRA" ) ;
4084 SpectralData hotspectra(arrayFloatCol.getColumn());
4085 TcalData tcaldata(sky);
4086 TsysData tsysdata(sky);
4087 unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ;
4088 // I know that the data is contiguous
4089 const uInt *p = rows.data() ;
4090 vector<int> ids( 2 ) ;
4091 Block<uInt> flagchan( spsize ) ;
4092 uInt nflag = 0 ;
4093 for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
4094 double reftime = timeCol.asdouble(*p) ;
4095 ids = getRowIdFromTime( reftime, timeOff ) ;
4096 Vector<Float> spoff = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeOff, ids, offspectra, "linear");
4097 ids = getRowIdFromTime( reftime, timeSky ) ;
4098 Vector<Float> spsky = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSky, ids, skyspectra, "linear");
4099 Vector<Float> tcal = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSky, ids, tcaldata, "linear");
4100 Vector<Float> tsys = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSky, ids, tsysdata, "linear");
4101 ids = getRowIdFromTime( reftime, timeHot ) ;
4102 Vector<Float> sphot = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHot, ids, hotspectra, "linear");
4103 Vector<Float> spec = on->specCol_( *p ) ;
4104 if ( antname.find( "APEX" ) != String::npos ) {
4105 // using gain array
4106 for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
4107 if ( spoff[j] == 0.0 || (sphot[j]-spsky[j]) == 0.0 ) {
4108 spec[j] = 0.0 ;
4109 flagchan[nflag++] = j ;
4110 }
4111 else {
4112 spec[j] = ( ( spec[j] - spoff[j] ) / spoff[j] )
4113 * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
4114 }
4115 }
4116 }
4117 else {
4118 // Chopper-Wheel calibration (Ulich & Haas 1976)
4119 for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
4120 if ( (sphot[j]-spsky[j]) == 0.0 ) {
4121 spec[j] = 0.0 ;
4122 flagchan[nflag++] = j ;
4123 }
4124 else {
4125 spec[j] = ( spec[j] - spoff[j] ) / ( sphot[j] - spsky[j] ) * tcal[j] ;
4126 }
4127 }
4128 }
4129 out->specCol_.put( *p, spec ) ;
4130 out->tsysCol_.put( *p, tsys ) ;
4131 if ( nflag > 0 ) {
4132 Vector<uChar> fl = out->flagsCol_( *p ) ;
4133 for ( unsigned int j = 0 ; j < nflag ; j++ ) {
4134 fl[flagchan[j]] = (uChar)True ;
4135 }
4136 out->flagsCol_.put( *p, fl ) ;
4137 }
4138 nflag = 0 ;
4139 p++ ;
4140 }
4141}
4142
4143void STMath::calibrateALMA( CountedPtr<Scantable>& out,
4144 const CountedPtr<Scantable>& on,
4145 const CountedPtr<Scantable>& off,
4146 const Vector<uInt>& rows )
4147{
4148 // 2012/05/22 TN
4149 // Assume that out has empty SPECTRA column
4150
4151 // if rows is empty, just return
4152 if ( rows.nelements() == 0 )
4153 return ;
4154 ROScalarColumn<Double> timeCol( off->table(), "TIME" ) ;
4155 Vector<Double> timeVec = timeCol.getColumn() ;
4156 timeCol.attach( on->table(), "TIME" ) ;
4157 ROArrayColumn<Float> arrayFloatCol( off->table(), "SPECTRA" ) ;
4158 SpectralData offspectra(arrayFloatCol.getColumn());
4159 unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ;
4160 // I know that the data is contiguous
4161 const uInt *p = rows.data() ;
4162 vector<int> ids( 2 ) ;
4163 Block<uInt> flagchan( spsize ) ;
4164 uInt nflag = 0 ;
4165 for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
4166 double reftime = timeCol.asdouble(*p) ;
4167 ids = getRowIdFromTime( reftime, timeVec ) ;
4168 Vector<Float> spoff = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeVec, ids, offspectra, "linear");
4169 Vector<Float> spec = on->specCol_( *p ) ;
4170 Vector<Float> tsys = on->tsysCol_( *p ) ;
4171 // ALMA Calibration
4172 //
4173 // Ta* = Tsys * ( ON - OFF ) / OFF
4174 //
4175 // 2010/01/07 Takeshi Nakazato
4176 unsigned int tsyssize = tsys.nelements() ;
4177 for ( unsigned int j = 0 ; j < spsize ; j++ ) {
4178 if ( spoff[j] == 0.0 ) {
4179 spec[j] = 0.0 ;
4180 flagchan[nflag++] = j ;
4181 }
4182 else {
4183 spec[j] = ( spec[j] - spoff[j] ) / spoff[j] ;
4184 }
4185 if ( tsyssize == spsize )
4186 spec[j] *= tsys[j] ;
4187 else
4188 spec[j] *= tsys[0] ;
4189 }
4190 out->specCol_.put( *p, spec ) ;
4191 if ( nflag > 0 ) {
4192 Vector<uChar> fl = out->flagsCol_( *p ) ;
4193 for ( unsigned int j = 0 ; j < nflag ; j++ ) {
4194 fl[flagchan[j]] = (uChar)True ;
4195 }
4196 out->flagsCol_.put( *p, fl ) ;
4197 }
4198 nflag = 0 ;
4199 p++ ;
4200 }
4201}
4202
4203void STMath::calibrateAPEXFS( CountedPtr<Scantable> &sig,
4204 CountedPtr<Scantable> &ref,
4205 const vector< CountedPtr<Scantable> >& on,
4206 const vector< CountedPtr<Scantable> >& sky,
4207 const vector< CountedPtr<Scantable> >& hot,
4208 const vector< CountedPtr<Scantable> >& cold,
4209 const Vector<uInt> &rows )
4210{
4211 // if rows is empty, just return
4212 if ( rows.nelements() == 0 )
4213 return ;
4214 ROScalarColumn<Double> timeCol( sky[0]->table(), "TIME" ) ;
4215 Vector<Double> timeSkyS = timeCol.getColumn() ;
4216 timeCol.attach( sky[1]->table(), "TIME" ) ;
4217 Vector<Double> timeSkyR = timeCol.getColumn() ;
4218 timeCol.attach( hot[0]->table(), "TIME" ) ;
4219 Vector<Double> timeHotS = timeCol.getColumn() ;
4220 timeCol.attach( hot[1]->table(), "TIME" ) ;
4221 Vector<Double> timeHotR = timeCol.getColumn() ;
4222 timeCol.attach( sig->table(), "TIME" ) ;
4223 ROScalarColumn<Double> timeCol2( ref->table(), "TIME" ) ;
4224 ROArrayColumn<Float> arrayFloatCol( sky[0]->table(), "SPECTRA" ) ;
4225 SpectralData skyspectraS(arrayFloatCol.getColumn());
4226 arrayFloatCol.attach( sky[1]->table(), "SPECTRA" ) ;
4227 SpectralData skyspectraR(arrayFloatCol.getColumn());
4228 arrayFloatCol.attach( hot[0]->table(), "SPECTRA" ) ;
4229 SpectralData hotspectraS(arrayFloatCol.getColumn());
4230 arrayFloatCol.attach( hot[1]->table(), "SPECTRA" ) ;
4231 SpectralData hotspectraR(arrayFloatCol.getColumn());
4232 TcalData tcaldataS(sky[0]);
4233 TsysData tsysdataS(sky[0]);
4234 TcalData tcaldataR(sky[1]);
4235 TsysData tsysdataR(sky[1]);
4236 unsigned int spsize = sig->nchan( sig->getIF(rows[0]) ) ;
4237 Vector<Float> spec( spsize ) ;
4238 // I know that the data is contiguous
4239 const uInt *p = rows.data() ;
4240 vector<int> ids( 2 ) ;
4241 Block<uInt> flagchan( spsize ) ;
4242 uInt nflag = 0 ;
4243 for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
4244 double reftime = timeCol.asdouble(*p) ;
4245 ids = getRowIdFromTime( reftime, timeSkyS ) ;
4246 Vector<Float> spskyS = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSkyS, ids, skyspectraS, "linear");
4247 Vector<Float> tcalS = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSkyS, ids, tcaldataS, "linear");
4248 Vector<Float> tsysS = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSkyS, ids, tsysdataS, "linear");
4249 ids = getRowIdFromTime( reftime, timeHotS ) ;
4250 Vector<Float> sphotS = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHotS, ids, hotspectraS, "linear");
4251 reftime = timeCol2.asdouble(*p) ;
4252 ids = getRowIdFromTime( reftime, timeSkyR ) ;
4253 Vector<Float> spskyR = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSkyR, ids, skyspectraR, "linear");
4254 Vector<Float> tcalR = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSkyR, ids, tcaldataR, "linear");
4255 Vector<Float> tsysR = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSkyR, ids, tsysdataR, "linear");
4256 ids = getRowIdFromTime( reftime, timeHotR ) ;
4257 Vector<Float> sphotR = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHotR, ids, hotspectraR, "linear");
4258 Vector<Float> spsig = on[0]->specCol_( *p ) ;
4259 Vector<Float> spref = on[1]->specCol_( *p ) ;
4260 for ( unsigned int j = 0 ; j < spsize ; j++ ) {
4261 if ( (sphotS[j]-spskyS[j]) == 0.0 || (sphotR[j]-spskyR[j]) == 0.0 ) {
4262 spec[j] = 0.0 ;
4263 flagchan[nflag++] = j ;
4264 }
4265 else {
4266 spec[j] = tcalS[j] * spsig[j] / ( sphotS[j] - spskyS[j] )
4267 - tcalR[j] * spref[j] / ( sphotR[j] - spskyR[j] ) ;
4268 }
4269 }
4270 sig->specCol_.put( *p, spec ) ;
4271 sig->tsysCol_.put( *p, tsysS ) ;
4272 spec *= (Float)-1.0 ;
4273 ref->specCol_.put( *p, spec ) ;
4274 ref->tsysCol_.put( *p, tsysR ) ;
4275 if ( nflag > 0 ) {
4276 Vector<uChar> flsig = sig->flagsCol_( *p ) ;
4277 Vector<uChar> flref = ref->flagsCol_( *p ) ;
4278 for ( unsigned int j = 0 ; j < nflag ; j++ ) {
4279 flsig[flagchan[j]] = (uChar)True ;
4280 flref[flagchan[j]] = (uChar)True ;
4281 }
4282 sig->flagsCol_.put( *p, flsig ) ;
4283 ref->flagsCol_.put( *p, flref ) ;
4284 }
4285 nflag = 0 ;
4286 p++ ;
4287 }
4288}
4289
4290void STMath::calibrateFS( CountedPtr<Scantable> &sig,
4291 CountedPtr<Scantable> &ref,
4292 const CountedPtr<Scantable>& rsig,
4293 const CountedPtr<Scantable>& rref,
4294 const CountedPtr<Scantable>& sky,
4295 const CountedPtr<Scantable>& hot,
4296 const CountedPtr<Scantable>& cold,
4297 const Vector<uInt> &rows )
4298{
4299 // if rows is empty, just return
4300 if ( rows.nelements() == 0 )
4301 return ;
4302 ROScalarColumn<Double> timeCol( sky->table(), "TIME" ) ;
4303 Vector<Double> timeSky = timeCol.getColumn() ;
4304 timeCol.attach( hot->table(), "TIME" ) ;
4305 Vector<Double> timeHot = timeCol.getColumn() ;
4306 timeCol.attach( sig->table(), "TIME" ) ;
4307 ROScalarColumn<Double> timeCol2( ref->table(), "TIME" ) ;
4308 ROArrayColumn<Float> arrayFloatCol( sky->table(), "SPECTRA" ) ;
4309 SpectralData skyspectra(arrayFloatCol.getColumn());
4310 arrayFloatCol.attach( hot->table(), "SPECTRA" ) ;
4311 SpectralData hotspectra(arrayFloatCol.getColumn());
4312 TcalData tcaldata(sky);
4313 TsysData tsysdata(sky);
4314 unsigned int spsize = sig->nchan( sig->getIF(rows[0]) ) ;
4315 Vector<Float> spec( spsize ) ;
4316 // I know that the data is contiguous
4317 const uInt *p = rows.data() ;
4318 vector<int> ids( 2 ) ;
4319 Block<uInt> flagchan( spsize ) ;
4320 uInt nflag = 0 ;
4321 for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
4322 double reftime = timeCol.asdouble(*p) ;
4323 ids = getRowIdFromTime( reftime, timeSky ) ;
4324 Vector<Float> spsky = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSky, ids, skyspectra, "linear");
4325 Vector<Float> tcal = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSky, ids, tcaldata, "linear");
4326 Vector<Float> tsys = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSky, ids, tsysdata, "linear");
4327 ids = getRowIdFromTime( reftime, timeHot ) ;
4328 Vector<Float> sphot = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHot, ids, hotspectra, "linear");
4329 Vector<Float> spsig = rsig->specCol_( *p ) ;
4330 Vector<Float> spref = rref->specCol_( *p ) ;
4331 // using gain array
4332 for ( unsigned int j = 0 ; j < spsize ; j++ ) {
4333 if ( spref[j] == 0.0 || (sphot[j]-spsky[j]) == 0.0 ) {
4334 spec[j] = 0.0 ;
4335 flagchan[nflag++] = j ;
4336 }
4337 else {
4338 spec[j] = ( ( spsig[j] - spref[j] ) / spref[j] )
4339 * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
4340 }
4341 }
4342 sig->specCol_.put( *p, spec ) ;
4343 sig->tsysCol_.put( *p, tsys ) ;
4344 if ( nflag > 0 ) {
4345 Vector<uChar> fl = sig->flagsCol_( *p ) ;
4346 for ( unsigned int j = 0 ; j < nflag ; j++ ) {
4347 fl[flagchan[j]] = (uChar)True ;
4348 }
4349 sig->flagsCol_.put( *p, fl ) ;
4350 }
4351 nflag = 0 ;
4352
4353 reftime = timeCol2.asdouble(*p) ;
4354 ids = getRowIdFromTime( reftime, timeSky ) ;
4355 spsky = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeSky, ids, skyspectra, "linear");
4356 tcal = SimpleInterpolationHelper<TcalData>::GetFromTime(reftime, timeSky, ids, tcaldata, "linear");
4357 tsys = SimpleInterpolationHelper<TsysData>::GetFromTime(reftime, timeSky, ids, tsysdata, "linear");
4358 ids = getRowIdFromTime( reftime, timeHot ) ;
4359 sphot = SimpleInterpolationHelper<SpectralData>::GetFromTime(reftime, timeHot, ids, hotspectra, "linear");
4360 // using gain array
4361 for ( unsigned int j = 0 ; j < spsize ; j++ ) {
4362 if ( spsig[j] == 0.0 || (sphot[j]-spsky[j]) == 0.0 ) {
4363 spec[j] = 0.0 ;
4364 flagchan[nflag++] = j ;
4365 }
4366 else {
4367 spec[j] = ( ( spref[j] - spsig[j] ) / spsig[j] )
4368 * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
4369 }
4370 }
4371 ref->specCol_.put( *p, spec ) ;
4372 ref->tsysCol_.put( *p, tsys ) ;
4373 if ( nflag > 0 ) {
4374 Vector<uChar> fl = ref->flagsCol_( *p ) ;
4375 for ( unsigned int j = 0 ; j < nflag ; j++ ) {
4376 fl[flagchan[j]] = (uChar)True ;
4377 }
4378 ref->flagsCol_.put( *p, fl ) ;
4379 }
4380 nflag = 0 ;
4381 p++ ;
4382 }
4383}
4384
4385void STMath::copyRows( Table &out,
4386 const Table &in,
4387 uInt startout,
4388 uInt startin,
4389 uInt nrow,
4390 Bool copySpectra,
4391 Bool copyFlagtra,
4392 Bool copyTsys )
4393{
4394 uInt nexclude = 0 ;
4395 Block<String> excludeColsBlock( 3 ) ;
4396 if ( !copySpectra ) {
4397 excludeColsBlock[nexclude] = "SPECTRA" ;
4398 nexclude++ ;
4399 }
4400 if ( !copyFlagtra ) {
4401 excludeColsBlock[nexclude] = "FLAGTRA" ;
4402 nexclude++ ;
4403 }
4404 if ( !copyTsys ) {
4405 excludeColsBlock[nexclude] = "TSYS" ;
4406 nexclude++ ;
4407 }
4408 // if ( nexclude < 3 ) {
4409 // excludeCols.resize( nexclude, True ) ;
4410 // }
4411 Vector<String> excludeCols( IPosition(1,nexclude),
4412 excludeColsBlock.storage(),
4413 SHARE ) ;
4414// cout << "excludeCols=" << excludeCols << endl ;
4415 TableRow rowout( out, excludeCols, True ) ;
4416 ROTableRow rowin( in, excludeCols, True ) ;
4417 uInt rin = startin ;
4418 uInt rout = startout ;
4419 for ( uInt i = 0 ; i < nrow ; i++ ) {
4420 rowin.get( rin ) ;
4421 rowout.putMatchingFields( rout, rowin.record() ) ;
4422 rin++ ;
4423 rout++ ;
4424 }
4425}
4426
4427CountedPtr<Scantable> STMath::averageWithinSession( CountedPtr<Scantable> &s,
4428 vector<bool> &mask,
4429 string weight )
4430{
4431 // prepare output table
4432 bool insitu = insitu_ ;
4433 insitu_ = false ;
4434 CountedPtr<Scantable> a = getScantable( s, true ) ;
4435 insitu_ = insitu ;
4436 Table &atab = a->table() ;
4437 ScalarColumn<Double> timeColOut( atab, "TIME" ) ;
4438
4439 if ( s->nrow() == 0 )
4440 return a ;
4441
4442 // setup RowAccumulator
4443 WeightType wtype = stringToWeight( weight ) ;
4444 RowAccumulator acc( wtype ) ;
4445 Vector<Bool> cmask( mask ) ;
4446 acc.setUserMask( cmask ) ;
4447
4448 vector<string> cols( 3 ) ;
4449 cols[0] = "IFNO" ;
4450 cols[1] = "POLNO" ;
4451 cols[2] = "BEAMNO" ;
4452 STIdxIter2 iter( s, cols ) ;
4453
4454 Table ttab = s->table() ;
4455 ROScalarColumn<Double> *timeCol = new ROScalarColumn<Double>( ttab, "TIME" ) ;
4456 Vector<Double> timeVec = timeCol->getColumn() ;
4457 delete timeCol ;
4458 Vector<Double> interval = s->integrCol_.getColumn() ;
4459 uInt nrow = timeVec.nelements() ;
4460 uInt outrow = 0 ;
4461
4462 while( !iter.pastEnd() ) {
4463
4464 Vector<uInt> rows = iter.getRows( SHARE ) ;
4465 uInt len = rows.nelements() ;
4466
4467 if ( len == 0 ) {
4468 iter.next() ;
4469 continue ;
4470 }
4471
4472 uInt nchan = s->nchan(s->getIF(rows[0])) ;
4473 Vector<uChar> flag( nchan ) ;
4474 Vector<Bool> bflag( nchan ) ;
4475 Vector<Float> spec( nchan ) ;
4476 Vector<Float> tsys( nchan ) ;
4477
4478 Vector<Double> timeSep( len-1 ) ;
4479 for ( uInt i = 0 ; i < len-1 ; i++ ) {
4480 timeSep[i] = timeVec[rows[i+1]] - timeVec[rows[i]] ;
4481 }
4482
4483 uInt irow ;
4484 uInt jrow ;
4485 for ( uInt i = 0 ; i < len-1 ; i++ ) {
4486 irow = rows[i] ;
4487 jrow = rows[i+1] ;
4488 // accumulate data
4489 s->flagsCol_.get( irow, flag ) ;
4490 convertArray( bflag, flag ) ;
4491 s->specCol_.get( irow, spec ) ;
4492 tsys.assign( s->tsysCol_( irow ) ) ;
4493 if ( !allEQ(bflag,True) )
4494 acc.add( spec, !bflag, tsys, interval[irow], timeVec[irow] ) ;
4495 double gap = 2.0 * 86400.0 * timeSep[i] / ( interval[jrow] + interval[irow] ) ;
4496 //cout << "gap[" << i << "]=" << setw(5) << gap << endl ;
4497 if ( gap > 1.1 ) {
4498 //cout << "detected gap between " << i << " and " << i+1 << endl ;
4499 // put data to output table
4500 // reset RowAccumulator
4501 if ( acc.state() ) {
4502 atab.addRow() ;
4503 copyRows( atab, ttab, outrow, irow, 1, False, False, False ) ;
4504 acc.replaceNaN() ;
4505 const Vector<Bool> &msk = acc.getMask() ;
4506 convertArray( flag, !msk ) ;
4507 for (uInt k = 0; k < nchan; ++k) {
4508 uChar userFlag = 1 << 7;
4509 if (msk[k]==True) userFlag = 0 << 7;
4510 flag(k) = userFlag;
4511 }
4512 a->flagsCol_.put( outrow, flag ) ;
4513 a->specCol_.put( outrow, acc.getSpectrum() ) ;
4514 a->tsysCol_.put( outrow, acc.getTsys() ) ;
4515 a->integrCol_.put( outrow, acc.getInterval() ) ;
4516 timeColOut.put( outrow, acc.getTime() ) ;
4517 a->cycleCol_.put( outrow, 0 ) ;
4518 }
4519 acc.reset() ;
4520 outrow++ ;
4521 }
4522 }
4523
4524 // accumulate and add last data
4525 irow = rows[len-1] ;
4526 s->flagsCol_.get( irow, flag ) ;
4527 convertArray( bflag, flag ) ;
4528 s->specCol_.get( irow, spec ) ;
4529 tsys.assign( s->tsysCol_( irow ) ) ;
4530 if (!allEQ(bflag,True) )
4531 acc.add( spec, !bflag, tsys, interval[irow], timeVec[irow] ) ;
4532 if ( acc.state() ) {
4533 atab.addRow() ;
4534 copyRows( atab, ttab, outrow, irow, 1, False, False, False ) ;
4535 acc.replaceNaN() ;
4536 const Vector<Bool> &msk = acc.getMask() ;
4537 convertArray( flag, !msk ) ;
4538 for (uInt k = 0; k < nchan; ++k) {
4539 uChar userFlag = 1 << 7;
4540 if (msk[k]==True) userFlag = 0 << 7;
4541 flag(k) = userFlag;
4542 }
4543 a->flagsCol_.put( outrow, flag ) ;
4544 a->specCol_.put( outrow, acc.getSpectrum() ) ;
4545 a->tsysCol_.put( outrow, acc.getTsys() ) ;
4546 a->integrCol_.put( outrow, acc.getInterval() ) ;
4547 timeColOut.put( outrow, acc.getTime() ) ;
4548 a->cycleCol_.put( outrow, 0 ) ;
4549 }
4550 acc.reset() ;
4551 outrow++ ;
4552
4553 iter.next() ;
4554 }
4555
4556 return a ;
4557}
Note: See TracBrowser for help on using the repository browser.