source: branches/hpc33/src/STMath.cpp@ 2450

Last change on this file since 2450 was 2449, checked in by KohjiNakamura, 13 years ago

probably faster TableExprNode

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