source: trunk/src/STMath.cpp@ 3102

Last change on this file since 3102 was 3089, checked in by Kana Sugimoto, 9 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s):

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