source: trunk/src/STMath.cpp@ 2873

Last change on this file since 2873 was 2829, checked in by Kana Sugimoto, 11 years ago

New Development: No

JIRA Issue: Yes

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): asap.scantable.average_time, asap.asapmath.average_time,

and all SD tasks that does averaging.

Description: STMath::average function now takes FLAGROW of scantable

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