source: trunk/src/STMath.cpp@ 2984

Last change on this file since 2984 was 2978, checked in by WataruKawasaki, 10 years ago

New Development: No

JIRA Issue: Yes CAS-6584, CAS-6787

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs: test_sdcal

Put in Release Notes:

Module(s): sd

Description: (1) modified STMath::almacal(), STMath::cwcal(), and the relevant functions so that flag information are properly handled in sdcal's 'ps'(for ALMA and APEX), 'otf', and 'otfraster' modes.

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