source: branches/casa-release-4_3-test02/src/STMath.cpp@ 3113

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

New Development: No

JIRA Issue: Yes CAS-6582

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...

STMath::average is changed so that empty averaged result (data to be
averaged are all flagged) is kept with all channel flags of 128 and
row flag of 1.


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