source: trunk/src/STMath.cpp@ 3017

Last change on this file since 3017 was 3010, checked in by WataruKawasaki, 10 years ago

New Development: No

JIRA Issue: Yes CAS-6599

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): sd

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