source: trunk/src/STMath.cpp@ 2953

Last change on this file since 2953 was 2952, checked in by WataruKawasaki, 10 years ago

New Development: No

JIRA Issue: Yes CAS-6598

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: add a parameter for some functions

Test Programs: test_sdscale

Put in Release Notes:

Module(s): sd

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