source: trunk/src/STMath.cpp@ 2921

Last change on this file since 2921 was 2920, checked in by Takeshi Nakazato, 10 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here..

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