source: trunk/src/STMath.cpp@ 2977

Last change on this file since 2977 was 2974, checked in by WataruKawasaki, 10 years ago

New Development: No

JIRA Issue: Yes CAS-6594

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs: test_sdmath

Put in Release Notes: No

Module(s): sd

Description: modified STMath::binaryOperate() so that:

(1) If one of or both spectra are row-flagged,

operation will be skipped so the 'left'
spectrum values are not modified, while
the FLAGROW of the output spectrum will
be set flagged.

(2) Otherwise, operation is executed for all

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