source: trunk/src/STMath.cpp@ 3117

Last change on this file since 3117 was 3106, checked in by Takeshi Nakazato, 8 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes/No

Interface Changes: Yes/No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...


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