source: trunk/src/STMath.cpp@ 2897

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

New Development: No (IMPORTANT bug fix)

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): scantable.smooth

Description: Fixed a bug that flag is ignored in smoothing.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 163.2 KB
Line 
1//
2// C++ Implementation: STMath
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12
13#include <sstream>
14
15#include <casa/iomanip.h>
16#include <casa/Arrays/MaskArrLogi.h>
17#include <casa/Arrays/MaskArrMath.h>
18#include <casa/Arrays/ArrayLogical.h>
19#include <casa/Arrays/ArrayMath.h>
20#include <casa/Arrays/Slice.h>
21#include <casa/Arrays/Slicer.h>
22#include <casa/BasicSL/String.h>
23#include <casa/Containers/Block.h>
24#include <casa/Containers/RecordField.h>
25#include <casa/Exceptions/Error.h>
26#include <casa/Logging/LogIO.h>
27
28#include <coordinates/Coordinates/CoordinateSystem.h>
29#include <coordinates/Coordinates/CoordinateUtil.h>
30#include <coordinates/Coordinates/FrequencyAligner.h>
31#include <coordinates/Coordinates/SpectralCoordinate.h>
32
33#include <lattices/Lattices/LatticeUtilities.h>
34
35#include <scimath/Functionals/Polynomial.h>
36#include <scimath/Mathematics/Convolver.h>
37#include <scimath/Mathematics/VectorKernel.h>
38
39#include <tables/Tables/ExprNode.h>
40#include <tables/Tables/ReadAsciiTable.h>
41#include <tables/Tables/TableCopy.h>
42#include <tables/Tables/TableIter.h>
43#include <tables/Tables/TableParse.h>
44#include <tables/Tables/TableRecord.h>
45#include <tables/Tables/TableRow.h>
46#include <tables/Tables/TableVector.h>
47#include <tables/Tables/TabVecMath.h>
48
49#include <atnf/PKSIO/SrcType.h>
50
51#include "RowAccumulator.h"
52#include "STAttr.h"
53#include "STMath.h"
54#include "STSelector.h"
55#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 mathutil::replaceMaskByZero(spec, !mask);
2530 //std::vector<bool> vmask;
2531 //(!mask).tovector(vmask);
2532 //mathutil::doZeroOrderInterpolation(spec, vmask);
2533 conv.linearConv(specout, spec);
2534 specCol.put(i, specout);
2535 }
2536 ++iter;
2537 }
2538 return out;
2539}
2540
2541CountedPtr< Scantable >
2542 STMath::merge( const std::vector< CountedPtr < Scantable > >& in )
2543{
2544 if ( in.size() < 2 ) {
2545 throw(AipsError("Need at least two scantables to perform a merge."));
2546 }
2547 std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
2548 bool insitu = insitu_;
2549 setInsitu(false);
2550 CountedPtr< Scantable > out = getScantable(*it, false);
2551 setInsitu(insitu);
2552 Table& tout = out->table();
2553 ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
2554 ScalarColumn<uInt> scannocol(tout,"SCANNO"), focusidcol(tout,"FOCUS_ID");
2555 // Renumber SCANNO to be 0-based
2556 Vector<uInt> scannos = scannocol.getColumn();
2557 uInt offset = min(scannos);
2558 scannos -= offset;
2559 scannocol.putColumn(scannos);
2560 uInt newscanno = max(scannos)+1;
2561 ++it;
2562 while ( it != in.end() ){
2563 if ( ! (*it)->conformant(*out) ) {
2564 // non conformant.
2565 LogIO os( LogOrigin( "STMath", "merge()", WHERE ) ) ;
2566 os << LogIO::SEVERE << "Can't merge scantables as header informations (any one of AntennaName, Equinox, and FluxUnit) differ." << LogIO::EXCEPTION ;
2567 }
2568 out->appendToHistoryTable((*it)->history());
2569 const Table& tab = (*it)->table();
2570
2571 Block<String> cols(3);
2572 cols[0] = String("FREQ_ID");
2573 cols[1] = String("MOLECULE_ID");
2574 cols[2] = String("FOCUS_ID");
2575
2576 TableIterator scanit(tab, "SCANNO");
2577 while (!scanit.pastEnd()) {
2578 ScalarColumn<uInt> thescannocol(scanit.table(),"SCANNO");
2579 Vector<uInt> thescannos(thescannocol.nrow(),newscanno);
2580 thescannocol.putColumn(thescannos);
2581 TableIterator subit(scanit.table(), cols);
2582 while ( !subit.pastEnd() ) {
2583 uInt nrow = tout.nrow();
2584 Table thetab = subit.table();
2585 ROTableRow row(thetab);
2586 Vector<uInt> thecolvals(thetab.nrow());
2587 ScalarColumn<uInt> thefreqidcol(thetab,"FREQ_ID");
2588 ScalarColumn<uInt> themolidcol(thetab, "MOLECULE_ID");
2589 ScalarColumn<uInt> thefocusidcol(thetab,"FOCUS_ID");
2590 // The selected subset of table should have
2591 // the equal FREQ_ID, MOLECULE_ID, and FOCUS_ID values.
2592 const TableRecord& rec = row.get(0);
2593 // Set the proper FREQ_ID
2594 Double rv,rp,inc;
2595 (*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
2596 uInt id;
2597 id = out->frequencies().addEntry(rp, rv, inc);
2598 thecolvals = id;
2599 thefreqidcol.putColumn(thecolvals);
2600 // Set the proper MOLECULE_ID
2601 Vector<String> name,fname;Vector<Double> rf;
2602 (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
2603 id = out->molecules().addEntry(rf, name, fname);
2604 thecolvals = id;
2605 themolidcol.putColumn(thecolvals);
2606 // Set the proper FOCUS_ID
2607 Float fpa,frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
2608 (*it)->focus().getEntry(fpa, fax, ftan, frot, fhand, fmount,fuser,
2609 fxy, fxyp, rec.asuInt("FOCUS_ID"));
2610 id = out->focus().addEntry(fpa, fax, ftan, frot, fhand, fmount,fuser,
2611 fxy, fxyp);
2612 thecolvals = id;
2613 thefocusidcol.putColumn(thecolvals);
2614
2615 tout.addRow(thetab.nrow());
2616 TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
2617
2618 ++subit;
2619 }
2620 ++newscanno;
2621 ++scanit;
2622 }
2623 ++it;
2624 }
2625 return out;
2626}
2627
2628CountedPtr< Scantable >
2629 STMath::invertPhase( const CountedPtr < Scantable >& in )
2630{
2631 return applyToPol(in, &STPol::invertPhase, Float(0.0));
2632}
2633
2634CountedPtr< Scantable >
2635 STMath::rotateXYPhase( const CountedPtr < Scantable >& in, float phase )
2636{
2637 return applyToPol(in, &STPol::rotatePhase, Float(phase));
2638}
2639
2640CountedPtr< Scantable >
2641 STMath::rotateLinPolPhase( const CountedPtr < Scantable >& in, float phase )
2642{
2643 return applyToPol(in, &STPol::rotateLinPolPhase, Float(phase));
2644}
2645
2646CountedPtr< Scantable > STMath::applyToPol( const CountedPtr<Scantable>& in,
2647 STPol::polOperation fptr,
2648 Float phase )
2649{
2650 CountedPtr< Scantable > out = getScantable(in, false);
2651 Table& tout = out->table();
2652 Block<String> cols(4);
2653 cols[0] = String("SCANNO");
2654 cols[1] = String("BEAMNO");
2655 cols[2] = String("IFNO");
2656 cols[3] = String("CYCLENO");
2657 TableIterator iter(tout, cols);
2658 CountedPtr<STPol> stpol = STPol::getPolClass(out->factories_,
2659 out->getPolType() );
2660 while (!iter.pastEnd()) {
2661 Table t = iter.table();
2662 ArrayColumn<Float> speccol(t, "SPECTRA");
2663 ScalarColumn<uInt> focidcol(t, "FOCUS_ID");
2664 Matrix<Float> pols(speccol.getColumn());
2665 try {
2666 stpol->setSpectra(pols);
2667 Float fang,fhand;
2668 fang = in->focusTable_.getTotalAngle(focidcol(0));
2669 fhand = in->focusTable_.getFeedHand(focidcol(0));
2670 stpol->setPhaseCorrections(fang, fhand);
2671 // use a member function pointer in STPol. This only works on
2672 // the STPol pointer itself, not the Counted Pointer so
2673 // derefernce it.
2674 (&(*(stpol))->*fptr)(phase);
2675 speccol.putColumn(stpol->getSpectra());
2676 } catch (AipsError& e) {
2677 //delete stpol;stpol=0;
2678 throw(e);
2679 }
2680 ++iter;
2681 }
2682 //delete stpol;stpol=0;
2683 return out;
2684}
2685
2686CountedPtr< Scantable >
2687 STMath::swapPolarisations( const CountedPtr< Scantable > & in )
2688{
2689 CountedPtr< Scantable > out = getScantable(in, false);
2690 Table& tout = out->table();
2691 Table t0 = tout(tout.col("POLNO") == 0);
2692 Table t1 = tout(tout.col("POLNO") == 1);
2693 if ( t0.nrow() != t1.nrow() )
2694 throw(AipsError("Inconsistent number of polarisations"));
2695 ArrayColumn<Float> speccol0(t0, "SPECTRA");
2696 ArrayColumn<uChar> flagcol0(t0, "FLAGTRA");
2697 ArrayColumn<Float> speccol1(t1, "SPECTRA");
2698 ArrayColumn<uChar> flagcol1(t1, "FLAGTRA");
2699 Matrix<Float> s0 = speccol0.getColumn();
2700 Matrix<uChar> f0 = flagcol0.getColumn();
2701 speccol0.putColumn(speccol1.getColumn());
2702 flagcol0.putColumn(flagcol1.getColumn());
2703 speccol1.putColumn(s0);
2704 flagcol1.putColumn(f0);
2705 return out;
2706}
2707
2708CountedPtr< Scantable >
2709 STMath::averagePolarisations( const CountedPtr< Scantable > & in,
2710 const std::vector<bool>& mask,
2711 const std::string& weight )
2712{
2713 if (in->npol() < 2 )
2714 throw(AipsError("averagePolarisations can only be applied to two or more"
2715 "polarisations"));
2716 bool insitu = insitu_;
2717 setInsitu(false);
2718 CountedPtr< Scantable > pols = getScantable(in, true);
2719 setInsitu(insitu);
2720 Table& tout = pols->table();
2721 std::string taql = "SELECT FROM $1 WHERE POLNO IN [0,1]";
2722 Table tab = tableCommand(taql, in->table());
2723 if (tab.nrow() == 0 )
2724 throw(AipsError("Could not find any rows with POLNO==0 and POLNO==1"));
2725 TableCopy::copyRows(tout, tab);
2726 TableVector<uInt> vec(tout, "POLNO");
2727 vec = 0;
2728 pols->table_.rwKeywordSet().define("nPol", Int(1));
2729 pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
2730 //pols->table_.rwKeywordSet().define("POLTYPE", in->getPolType());
2731 std::vector<CountedPtr<Scantable> > vpols;
2732 vpols.push_back(pols);
2733 CountedPtr< Scantable > out = average(vpols, mask, weight, "SCAN");
2734 return out;
2735}
2736
2737CountedPtr< Scantable >
2738 STMath::averageBeams( const CountedPtr< Scantable > & in,
2739 const std::vector<bool>& mask,
2740 const std::string& weight )
2741{
2742 bool insitu = insitu_;
2743 setInsitu(false);
2744 CountedPtr< Scantable > beams = getScantable(in, false);
2745 setInsitu(insitu);
2746 Table& tout = beams->table();
2747 // give all rows the same BEAMNO
2748 TableVector<uInt> vec(tout, "BEAMNO");
2749 vec = 0;
2750 beams->table_.rwKeywordSet().define("nBeam", Int(1));
2751 std::vector<CountedPtr<Scantable> > vbeams;
2752 vbeams.push_back(beams);
2753 CountedPtr< Scantable > out = average(vbeams, mask, weight, "SCAN");
2754 return out;
2755}
2756
2757
2758CountedPtr< Scantable >
2759 asap::STMath::frequencyAlign( const CountedPtr< Scantable > & in,
2760 const std::string & refTime,
2761 const std::string & method)
2762{
2763 LogIO os( casa::LogOrigin("STMath", "frequencyAlign()", WHERE));
2764 // clone as this is not working insitu
2765 bool insitu = insitu_;
2766 setInsitu(false);
2767 CountedPtr< Scantable > out = getScantable(in, false);
2768 setInsitu(insitu);
2769 Table& tout = out->table();
2770 // Get reference Epoch to time of first row or given String
2771 Unit DAY(String("d"));
2772 MEpoch::Ref epochRef(in->getTimeReference());
2773 MEpoch refEpoch;
2774 if (refTime.length()>0) {
2775 Quantum<Double> qt;
2776 if (MVTime::read(qt,refTime)) {
2777 MVEpoch mv(qt);
2778 refEpoch = MEpoch(mv, epochRef);
2779 } else {
2780 throw(AipsError("Invalid format for Epoch string"));
2781 }
2782 } else {
2783 refEpoch = in->timeCol_(0);
2784 }
2785 MPosition refPos = in->getAntennaPosition();
2786
2787 InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
2788 /*
2789 // Comment from MV.
2790 // the following code has been commented out because different FREQ_IDs have to be aligned together even
2791 // if the frame doesn't change. So far, lack of this check didn't cause any problems.
2792 // test if user frame is different to base frame
2793 if ( in->frequencies().getFrameString(true)
2794 == in->frequencies().getFrameString(false) ) {
2795 throw(AipsError("Can't convert as no output frame has been set"
2796 " (use set_freqframe) or it is aligned already."));
2797 }
2798 */
2799 MFrequency::Types system = in->frequencies().getFrame();
2800 MVTime mvt(refEpoch.getValue());
2801 String epochout = mvt.string(MVTime::YMD) + String(" (") + refEpoch.getRefString() + String(")");
2802 os << "Aligning at reference Epoch " << epochout
2803 << " in frame " << MFrequency::showType(system) << LogIO::POST;
2804 // set up the iterator
2805 Block<String> cols(4);
2806 // select by constant direction
2807 cols[0] = String("SRCNAME");
2808 cols[1] = String("BEAMNO");
2809 // select by IF ( no of channels varies over this )
2810 cols[2] = String("IFNO");
2811 // select by restfrequency
2812 cols[3] = String("MOLECULE_ID");
2813 TableIterator iter(tout, cols);
2814 while ( !iter.pastEnd() ) {
2815 Table t = iter.table();
2816 ROScalarColumn<String> snCol(t, "SRCNAME");
2817 os << "Aligning to position of source '" << snCol(0) << "'" << LogIO::POST;
2818 MDirection::ROScalarColumn dirCol(t, "DIRECTION");
2819 TableIterator fiter(t, "FREQ_ID");
2820 // determine nchan from the first row. This should work as
2821 // we are iterating over BEAMNO and IFNO
2822 // we should have constant direction
2823
2824 ROArrayColumn<Float> sCol(t, "SPECTRA");
2825 const MDirection direction = dirCol(0);
2826 const uInt nchan = sCol(0).nelements();
2827
2828 // skip operations if there is nothing to align
2829 if (fiter.pastEnd()) {
2830 continue;
2831 }
2832
2833 Table ftab = fiter.table();
2834 // align all frequency ids with respect to the first encountered id
2835 ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
2836 // get the SpectralCoordinate for the freqid, which we are iterating over
2837 SpectralCoordinate sC = \
2838 in->frequencies().getSpectralCoordinate(freqidCol(0));
2839 FrequencyAligner<Float> fa( sC, nchan, refEpoch,
2840 direction, refPos, system );
2841 // realign the SpectralCoordinate and put into the output Scantable
2842 Vector<String> units(1);
2843 units = String("Hz");
2844 Bool linear=True;
2845 SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
2846 sc2.setWorldAxisUnits(units);
2847 const uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
2848 sc2.referenceValue()[0],
2849 sc2.increment()[0]);
2850 while ( !fiter.pastEnd() ) {
2851
2852 ftab = fiter.table();
2853 // spectral coordinate for the current FREQ_ID
2854 ScalarColumn<uInt> freqidCol2(ftab, "FREQ_ID");
2855 sC = in->frequencies().getSpectralCoordinate(freqidCol2(0));
2856 // create the "global" abcissa for alignment with same FREQ_ID
2857 Vector<Double> abc(nchan);
2858 for (uInt i=0; i<nchan; i++) {
2859 Double w;
2860 sC.toWorld(w,Double(i));
2861 abc[i] = w;
2862 }
2863 TableVector<uInt> tvec(ftab, "FREQ_ID");
2864 // assign new frequency id to all rows
2865 tvec = id;
2866 // cache abcissa for same time stamps, so iterate over those
2867 TableIterator timeiter(ftab, "TIME");
2868 while ( !timeiter.pastEnd() ) {
2869 Table tab = timeiter.table();
2870 ArrayColumn<Float> specCol(tab, "SPECTRA");
2871 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2872 MEpoch::ROScalarColumn timeCol(tab, "TIME");
2873 // use align abcissa cache after the first row
2874 // these rows should be just be POLNO
2875 bool first = true;
2876 for (int i=0; i<int(tab.nrow()); ++i) {
2877 // input values
2878 Vector<uChar> flag = flagCol(i);
2879 Vector<Bool> mask(flag.shape());
2880 Vector<Float> specOut, spec;
2881 spec = specCol(i);
2882 Vector<Bool> maskOut;Vector<uChar> flagOut;
2883 convertArray(mask, flag);
2884 // alignment
2885 Bool ok = fa.align(specOut, maskOut, abc, spec,
2886 mask, timeCol(i), !first,
2887 interp, False);
2888 (void) ok; // unused stop compiler nagging
2889 // back into scantable
2890 flagOut.resize(maskOut.nelements());
2891 convertArray(flagOut, maskOut);
2892 flagCol.put(i, flagOut);
2893 specCol.put(i, specOut);
2894 // start abcissa caching
2895 first = false;
2896 }
2897 // next timestamp
2898 ++timeiter;
2899 }
2900 // next FREQ_ID
2901 ++fiter;
2902 }
2903 // next aligner
2904 ++iter;
2905 }
2906 // set this afterwards to ensure we are doing insitu correctly.
2907 out->frequencies().setFrame(system, true);
2908 return out;
2909}
2910
2911CountedPtr<Scantable>
2912 asap::STMath::convertPolarisation( const CountedPtr<Scantable>& in,
2913 const std::string & newtype )
2914{
2915 if (in->npol() != 2 && in->npol() != 4)
2916 throw(AipsError("Can only convert two or four polarisations."));
2917 if ( in->getPolType() == newtype )
2918 throw(AipsError("No need to convert."));
2919 if ( ! in->selector_.empty() )
2920 throw(AipsError("Can only convert whole scantable. Unset the selection."));
2921 bool insitu = insitu_;
2922 setInsitu(false);
2923 CountedPtr< Scantable > out = getScantable(in, true);
2924 setInsitu(insitu);
2925 Table& tout = out->table();
2926 tout.rwKeywordSet().define("POLTYPE", String(newtype));
2927
2928 Block<String> cols(4);
2929 cols[0] = "SCANNO";
2930 cols[1] = "CYCLENO";
2931 cols[2] = "BEAMNO";
2932 cols[3] = "IFNO";
2933 TableIterator it(in->originalTable_, cols);
2934 String basetype = in->getPolType();
2935 STPol* stpol = STPol::getPolClass(in->factories_, basetype);
2936 try {
2937 while ( !it.pastEnd() ) {
2938 Table tab = it.table();
2939 uInt row = tab.rowNumbers()[0];
2940 stpol->setSpectra(in->getPolMatrix(row));
2941 Float fang,fhand;
2942 fang = in->focusTable_.getTotalAngle(in->mfocusidCol_(row));
2943 fhand = in->focusTable_.getFeedHand(in->mfocusidCol_(row));
2944 stpol->setPhaseCorrections(fang, fhand);
2945 Int npolout = 0;
2946 for (uInt i=0; i<tab.nrow(); ++i) {
2947 Vector<Float> outvec = stpol->getSpectrum(i, newtype);
2948 if ( outvec.nelements() > 0 ) {
2949 tout.addRow();
2950 TableCopy::copyRows(tout, tab, tout.nrow()-1, 0, 1);
2951 ArrayColumn<Float> sCol(tout,"SPECTRA");
2952 ScalarColumn<uInt> pCol(tout,"POLNO");
2953 sCol.put(tout.nrow()-1 ,outvec);
2954 pCol.put(tout.nrow()-1 ,uInt(npolout));
2955 npolout++;
2956 }
2957 }
2958 tout.rwKeywordSet().define("nPol", npolout);
2959 ++it;
2960 }
2961 } catch (AipsError& e) {
2962 delete stpol;
2963 throw(e);
2964 }
2965 delete stpol;
2966 return out;
2967}
2968
2969CountedPtr< Scantable >
2970 asap::STMath::mxExtract( const CountedPtr< Scantable > & in,
2971 const std::string & scantype )
2972{
2973 bool insitu = insitu_;
2974 setInsitu(false);
2975 CountedPtr< Scantable > out = getScantable(in, true);
2976 setInsitu(insitu);
2977 Table& tout = out->table();
2978 std::string taql = "SELECT FROM $1 WHERE BEAMNO != REFBEAMNO";
2979 if (scantype == "on") {
2980 taql = "SELECT FROM $1 WHERE BEAMNO == REFBEAMNO";
2981 }
2982 Table tab = tableCommand(taql, in->table());
2983 TableCopy::copyRows(tout, tab);
2984 if (scantype == "on") {
2985 // re-index SCANNO to 0
2986 TableVector<uInt> vec(tout, "SCANNO");
2987 vec = 0;
2988 }
2989 return out;
2990}
2991
2992std::vector<float>
2993 asap::STMath::fft( const casa::CountedPtr< Scantable > & in,
2994 const std::vector<int>& whichrow,
2995 bool getRealImag )
2996{
2997 std::vector<float> res;
2998 Table tab = in->table();
2999 std::vector<bool> mask;
3000
3001 if (whichrow.size() < 1) { // for all rows (by default)
3002 int nrow = int(tab.nrow());
3003 for (int i = 0; i < nrow; ++i) {
3004 res = in->execFFT(i, mask, getRealImag);
3005 }
3006 } else { // for specified rows
3007 for (uInt i = 0; i < whichrow.size(); ++i) {
3008 res = in->execFFT(i, mask, getRealImag);
3009 }
3010 }
3011
3012 return res;
3013}
3014
3015
3016CountedPtr<Scantable>
3017 asap::STMath::lagFlag( const CountedPtr<Scantable>& in,
3018 double start, double end,
3019 const std::string& mode )
3020{
3021 CountedPtr<Scantable> out = getScantable(in, false);
3022 Table& tout = out->table();
3023 TableIterator iter(tout, "FREQ_ID");
3024 FFTServer<Float,Complex> ffts;
3025
3026 while ( !iter.pastEnd() ) {
3027 Table tab = iter.table();
3028 Double rp,rv,inc;
3029 ROTableRow row(tab);
3030 const TableRecord& rec = row.get(0);
3031 uInt freqid = rec.asuInt("FREQ_ID");
3032 out->frequencies().getEntry(rp, rv, inc, freqid);
3033 ArrayColumn<Float> specCol(tab, "SPECTRA");
3034 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
3035
3036 for (int i=0; i<int(tab.nrow()); ++i) {
3037 Vector<Float> spec = specCol(i);
3038 Vector<uChar> flag = flagCol(i);
3039 std::vector<bool> mask;
3040 for (uInt j = 0; j < flag.nelements(); ++j) {
3041 mask.push_back(!(flag[j]>0));
3042 }
3043 mathutil::doZeroOrderInterpolation(spec, mask);
3044
3045 Vector<Complex> lags;
3046 ffts.fft0(lags, spec);
3047
3048 Int lag0(start+0.5);
3049 Int lag1(end+0.5);
3050 if (mode == "frequency") {
3051 lag0 = Int(spec.nelements()*abs(inc)/(start)+0.5);
3052 lag1 = Int(spec.nelements()*abs(inc)/(end)+0.5);
3053 }
3054 Int lstart = max(0, lag0);
3055 Int lend = min(Int(lags.nelements()-1), lag1);
3056 if (lstart == lend) {
3057 lags[lstart] = Complex(0.0);
3058 } else {
3059 if (lstart > lend) {
3060 Int tmp = lend;
3061 lend = lstart;
3062 lstart = tmp;
3063 }
3064 for (int j=lstart; j <=lend ;++j) {
3065 lags[j] = Complex(0.0);
3066 }
3067 }
3068
3069 ffts.fft0(spec, lags);
3070
3071 specCol.put(i, spec);
3072 }
3073 ++iter;
3074 }
3075 return out;
3076}
3077
3078// Averaging spectra with different channel/resolution
3079CountedPtr<Scantable>
3080STMath::new_average( const std::vector<CountedPtr<Scantable> >& in,
3081 const bool& compel,
3082 const std::vector<bool>& mask,
3083 const std::string& weight,
3084 const std::string& avmode )
3085 throw ( casa::AipsError )
3086{
3087 LogIO os( LogOrigin( "STMath", "new_average()", WHERE ) ) ;
3088 if ( avmode == "SCAN" && in.size() != 1 )
3089 throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
3090 "Use merge first."));
3091
3092 CountedPtr<Scantable> out ; // processed result
3093 if ( compel ) {
3094 std::vector< CountedPtr<Scantable> > newin ; // input for average process
3095 uInt insize = in.size() ; // number of input scantables
3096
3097 // setup newin
3098 bool oldInsitu = insitu_ ;
3099 setInsitu( false ) ;
3100 newin.resize( insize ) ;
3101 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3102 newin[itable] = getScantable( in[itable], false ) ;
3103 }
3104 setInsitu( oldInsitu ) ;
3105
3106 // warning
3107 os << "Average spectra with different spectral resolution" << LogIO::POST ;
3108
3109 // temporarily set coordinfo
3110 vector<string> oldinfo( insize ) ;
3111 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3112 vector<string> coordinfo = in[itable]->getCoordInfo() ;
3113 oldinfo[itable] = coordinfo[0] ;
3114 coordinfo[0] = "Hz" ;
3115 newin[itable]->setCoordInfo( coordinfo ) ;
3116 }
3117
3118 ostringstream oss ;
3119
3120 // check IF frequency coverage
3121 // freqid: list of FREQ_ID, which is used, in each table
3122 // iffreq: list of minimum and maximum frequency for each FREQ_ID in
3123 // each table
3124 // freqid[insize][numIF]
3125 // freqid: [[id00, id01, ...],
3126 // [id10, id11, ...],
3127 // ...
3128 // [idn0, idn1, ...]]
3129 // iffreq[insize][numIF*2]
3130 // iffreq: [[min_id00, max_id00, min_id01, max_id01, ...],
3131 // [min_id10, max_id10, min_id11, max_id11, ...],
3132 // ...
3133 // [min_idn0, max_idn0, min_idn1, max_idn1, ...]]
3134 //os << "Check IF settings in each table" << LogIO::POST ;
3135 vector< vector<uInt> > freqid( insize );
3136 vector< vector<double> > iffreq( insize ) ;
3137 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3138 Vector<uInt> freqIds = newin[itable]->mfreqidCol_.getColumn() ;
3139 vector<uInt> uniqueFreqId = newin[itable]->getNumbers(newin[itable]->mfreqidCol_) ;
3140 for ( vector<uInt>::iterator i = uniqueFreqId.begin() ;
3141 i != uniqueFreqId.end() ; i++ ) {
3142 //os << "itable = " << itable << ": IF " << id << " is included in the list" << LogIO::POST ;
3143 uInt target = 0 ;
3144 while ( freqIds[target] != *i )
3145 target++ ;
3146 vector<double> abcissa = newin[itable]->getAbcissa( target ) ;
3147 freqid[itable].push_back( *i ) ;
3148 double incr = abs( abcissa[1] - abcissa[0] ) ;
3149 iffreq[itable].push_back( (*min_element(abcissa.begin(),abcissa.end()))-0.5*incr ) ;
3150 iffreq[itable].push_back( (*max_element(abcissa.begin(),abcissa.end()))+0.5*incr ) ;
3151 }
3152 }
3153
3154 // debug
3155// os << "IF settings summary:" << endl ;
3156// for ( uInt i = 0 ; i < freqid.size() ; i++ ) {
3157// os << " Table" << i << endl ;
3158// for ( uInt j = 0 ; j < freqid[i].size() ; j++ ) {
3159// os << " id = " << freqid[i][j] << " (min,max) = (" << iffreq[i][2*j] << "," << iffreq[i][2*j+1] << ")" << endl ;
3160// }
3161// }
3162// os << endl ;
3163// os.post() ;
3164
3165 // IF grouping based on their frequency coverage
3166 // ifgrp: number of member in each IF group
3167 // ifgrp[numgrp]
3168 // ifgrp: [n0, n1, ...]
3169 //os << "IF grouping based on their frequency coverage" << LogIO::POST ;
3170
3171 // parameter for IF grouping
3172 // groupmode = OR retrieve all region
3173 // AND only retrieve overlaped region
3174 //string groupmode = "AND" ;
3175 string groupmode = "OR" ;
3176 uInt sizecr = 0 ;
3177 if ( groupmode == "AND" )
3178 sizecr = 1 ;
3179 else if ( groupmode == "OR" )
3180 sizecr = 0 ;
3181
3182 vector<double> sortedfreq ;
3183 for ( uInt i = 0 ; i < iffreq.size() ; i++ ) {
3184 for ( uInt j = 0 ; j < iffreq[i].size() ; j++ ) {
3185 if ( count( sortedfreq.begin(), sortedfreq.end(), iffreq[i][j] ) == 0 )
3186 sortedfreq.push_back( iffreq[i][j] ) ;
3187 }
3188 }
3189 sort( sortedfreq.begin(), sortedfreq.end() ) ;
3190 vector<uInt> ifgrp( sortedfreq.size()-1, 0 ) ;
3191 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3192 for ( uInt iif = 0 ; iif < freqid[itable].size() ; iif++ ) {
3193 double range0 = iffreq[itable][2*iif] ;
3194 double range1 = iffreq[itable][2*iif+1] ;
3195 for ( uInt j = 0 ; j < sortedfreq.size()-1 ; j++ ) {
3196 double fmin = max( range0, sortedfreq[j] ) ;
3197 double fmax = min( range1, sortedfreq[j+1] ) ;
3198 if ( fmin < fmax ) {
3199 ifgrp[j]++ ;
3200 }
3201 }
3202 }
3203 }
3204
3205 // Grouping continuous IF groups (without frequency gap)
3206 // freqgrp: list of IF group indexes in each frequency group
3207 // freqgrp[numgrp][nummember]
3208 // freqgrp: [[ifgrp00, ifgrp01, ifgrp02, ...],
3209 // [ifgrp10, ifgrp11, ifgrp12, ...],
3210 // ...
3211 // [ifgrpn0, ifgrpn1, ifgrpn2, ...]]
3212 // grprange[2*numgrp]
3213 // grprange: [fmin0,fmax0,fmin1,fmax1,...]
3214 vector< vector<uInt> > freqgrp ;
3215 vector<double> grprange ;
3216 vector<uInt> grpedge ;
3217 for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
3218 if ( ifgrp[igrp] <= sizecr ) {
3219 grpedge.push_back( igrp ) ;
3220 }
3221 }
3222 grpedge.push_back( ifgrp.size() ) ;
3223 uInt itmp = 0 ;
3224 for ( uInt i = 0 ; i < grpedge.size() ; i++ ) {
3225 int n = grpedge[i] - itmp ;
3226 if ( n > 0 ) {
3227 vector<uInt> members( n ) ;
3228 for ( int j = 0 ; j < n ; j++ ) {
3229 members[j] = itmp+j ;
3230 }
3231 freqgrp.push_back( members ) ;
3232 grprange.push_back( sortedfreq[itmp] ) ;
3233 grprange.push_back( sortedfreq[grpedge[i]] ) ;
3234 }
3235 itmp += n + 1 ;
3236 }
3237
3238 // print frequency group
3239 oss.str("") ;
3240 oss << "Frequency Group summary: " << endl ;
3241 oss << " GROUP_ID: [FREQ_MIN, FREQ_MAX]" << endl ;
3242 for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
3243 oss << " GROUP " << setw( 2 ) << i << ": [" << grprange[2*i] << "," << grprange[2*i+1] << "]" ;
3244 oss << endl ;
3245 }
3246 oss << endl ;
3247 os << oss.str() << LogIO::POST ;
3248
3249 // groups: list of frequency group index whose frequency range overlaps
3250 // with that of each table and IF
3251 // groups[numtable][numIF]
3252 // groups: [[grpx, grpy,...],
3253 // [grpa, grpb,...],
3254 // ...
3255 // [grpk, grpm,...]]
3256 vector< vector<uInt> > groups( insize ) ;
3257 for ( uInt i = 0 ; i < insize ; i++ ) {
3258 groups[i].resize( freqid[i].size() ) ;
3259 }
3260 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3261 for ( uInt ifreq = 0 ; ifreq < freqid[itable].size() ; ifreq++ ) {
3262 double minf = iffreq[itable][2*ifreq] ;
3263 uInt groupid ;
3264 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3265 vector<uInt> memberlist = freqgrp[igrp] ;
3266 if ( (minf >= grprange[2*igrp]) && (minf <= grprange[2*igrp+1]) ) {
3267 groupid = igrp ;
3268 break ;
3269 }
3270 }
3271 groups[itable][ifreq] = groupid ;
3272 }
3273 }
3274
3275
3276 // print membership
3277 oss.str("") ;
3278 for ( uInt i = 0 ; i < insize ; i++ ) {
3279 oss << "Table " << i << endl ;
3280 for ( uInt j = 0 ; j < groups[i].size() ; j++ ) {
3281 oss << " FREQ_ID " << setw( 2 ) << freqid[i][j] << ": " ;
3282 oss << setw( 2 ) << groups[i][j] ;
3283 oss << endl ;
3284 }
3285 }
3286 os << oss.str() << LogIO::POST ;
3287
3288 // reset SCANNO and IFNO/FREQ_ID: IF is reset by the result of sortation
3289 //os << "All IF number is set to IF group index" << LogIO::POST ;
3290 // reset SCANNO only when avmode != "SCAN"
3291 if ( avmode != "SCAN" ) {
3292 os << "All scan number is set to 0" << LogIO::POST ;
3293 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3294 uInt nrow = newin[itable]->nrow() ;
3295 Vector<uInt> resetScan( nrow, 0 ) ;
3296 newin[itable]->scanCol_.putColumn( resetScan ) ;
3297 }
3298 }
3299
3300 // reset spectra and flagtra: align spectral resolution
3301 //os << "Align spectral resolution" << LogIO::POST ;
3302 // gmaxdnu: the coarsest frequency resolution in the frequency group
3303 // gminfreq: lower frequency edge of the frequency group
3304 // gnchan: number of channels for the frequency group
3305 vector<double> gmaxdnu( freqgrp.size(), 0.0 ) ;
3306 vector<double> gminfreq( freqgrp.size() ) ;
3307 vector<double> gnchan( freqgrp.size() ) ;
3308 for ( uInt i = 0 ; i < insize ; i++ ) {
3309 vector<uInt> members = groups[i] ;
3310 for ( uInt j = 0 ; j < members.size() ; j++ ) {
3311 uInt groupid = members[j] ;
3312 Double rp,rv,ic ;
3313 newin[i]->frequencies().getEntry( rp, rv, ic, j ) ;
3314 if ( abs(ic) > abs(gmaxdnu[groupid]) )
3315 gmaxdnu[groupid] = ic ;
3316 }
3317 }
3318 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3319 gminfreq[igrp] = grprange[2*igrp] ;
3320 double maxfreq = grprange[2*igrp+1] ;
3321 gnchan[igrp] = (int)(abs((maxfreq-gminfreq[igrp])/gmaxdnu[igrp])+0.9) ;
3322 }
3323
3324 // regrid spectral data and update frequency info
3325 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3326 Vector<uInt> oldFreqId = newin[itable]->mfreqidCol_.getColumn() ;
3327 Vector<uInt> newFreqId( oldFreqId.nelements() ) ;
3328
3329 // update MAIN
3330 for ( uInt irow = 0 ; irow < newin[itable]->nrow() ; irow++ ) {
3331 uInt groupid = groups[itable][oldFreqId[irow]] ;
3332 newFreqId[irow] = groupid ;
3333 newin[itable]->regridChannel( gnchan[groupid],
3334 gmaxdnu[groupid],
3335 gminfreq[groupid],
3336 irow ) ;
3337 }
3338 newin[itable]->mfreqidCol_.putColumn( newFreqId ) ;
3339 newin[itable]->ifCol_.putColumn( newFreqId ) ;
3340
3341 // update FREQUENCIES
3342 Table tab = newin[itable]->frequencies().table() ;
3343 ScalarColumn<uInt> fIdCol( tab, "ID" ) ;
3344 ScalarColumn<Double> fRefPixCol( tab, "REFPIX" ) ;
3345 ScalarColumn<Double> fRefValCol( tab, "REFVAL" ) ;
3346 ScalarColumn<Double> fIncrCol( tab, "INCREMENT" ) ;
3347 if ( freqgrp.size() > tab.nrow() ) {
3348 tab.addRow( freqgrp.size()-tab.nrow() ) ;
3349 }
3350 for ( uInt irow = 0 ; irow < freqgrp.size() ; irow++ ) {
3351 Double refval = gminfreq[irow] + 0.5 * abs(gmaxdnu[irow]) ;
3352 Double refpix = (gmaxdnu[irow] > 0.0) ? 0 : gnchan[irow]-1 ;
3353 Double increment = gmaxdnu[irow] ;
3354 fIdCol.put( irow, irow ) ;
3355 fRefPixCol.put( irow, refpix ) ;
3356 fRefValCol.put( irow, refval ) ;
3357 fIncrCol.put( irow, increment ) ;
3358 }
3359 }
3360
3361 // set back coordinfo
3362 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3363 vector<string> coordinfo = newin[itable]->getCoordInfo() ;
3364 coordinfo[0] = oldinfo[itable] ;
3365 newin[itable]->setCoordInfo( coordinfo ) ;
3366 }
3367
3368 // average
3369 out = average( newin, mask, weight, avmode ) ;
3370 }
3371 else {
3372 // simple average
3373 out = average( in, mask, weight, avmode ) ;
3374 }
3375
3376 return out;
3377}
3378
3379CountedPtr<Scantable> STMath::cwcal( const CountedPtr<Scantable>& s,
3380 const String calmode,
3381 const String antname )
3382{
3383 // frequency switch
3384 if ( calmode == "fs" ) {
3385 return cwcalfs( s, antname ) ;
3386 }
3387 else {
3388 vector<bool> masks = s->getMask( 0 ) ;
3389 vector<int> types ;
3390
3391 // save original table selection
3392 Table torg = s->table_ ;
3393
3394 // sky scan
3395 bool insitu = insitu_ ;
3396 insitu_ = false ;
3397 // share calibration scans before average with out
3398 CountedPtr<Scantable> out = getScantable( s, true ) ;
3399 insitu_ = insitu ;
3400 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::SKY ) ;
3401 out->attach() ;
3402 CountedPtr<Scantable> asky = averageWithinSession( out,
3403 masks,
3404 "TINT" ) ;
3405 // hot scan
3406 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::HOT ) ;
3407 out->attach() ;
3408 CountedPtr<Scantable> ahot = averageWithinSession( out,
3409 masks,
3410 "TINT" ) ;
3411 // cold scan
3412 CountedPtr<Scantable> acold ;
3413// out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::COLD ) ;
3414// out->attach() ;
3415// CountedPtr<Scantable> acold = averageWithinSession( out,
3416// masks,
3417// "TINT" ) ;
3418
3419 // off scan
3420 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSOFF ) ;
3421 out->attach() ;
3422 CountedPtr<Scantable> aoff = averageWithinSession( out,
3423 masks,
3424 "TINT" ) ;
3425
3426 // on scan
3427 s->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSON ) ;
3428 s->attach() ;
3429 out->table_ = out->originalTable_ ;
3430 out->attach() ;
3431 out->table().addRow( s->nrow() ) ;
3432 copyRows( out->table(), s->table(), 0, 0, s->nrow(), False, True, False ) ;
3433
3434 // process each on scan
3435 STSelector sel ;
3436 vector<string> cols( 3 ) ;
3437 cols[0] = "BEAMNO" ;
3438 cols[1] = "POLNO" ;
3439 cols[2] = "IFNO" ;
3440 STIdxIter *iter = new STIdxIterAcc( out, cols ) ;
3441 while ( !iter->pastEnd() ) {
3442 Vector<uInt> ids = iter->current() ;
3443 stringstream ss ;
3444 ss << "SELECT FROM $1 WHERE "
3445 << "BEAMNO==" << ids[0] << "&&"
3446 << "POLNO==" << ids[1] << "&&"
3447 << "IFNO==" << ids[2] ;
3448 //cout << "TaQL string: " << ss.str() << endl ;
3449 sel.setTaQL( ss.str() ) ;
3450 aoff->setSelection( sel ) ;
3451 ahot->setSelection( sel ) ;
3452 asky->setSelection( sel ) ;
3453 Vector<uInt> rows = iter->getRows( SHARE ) ;
3454 // out should be an exact copy of s except that SPECTRA column is empty
3455 calibrateCW( out, s, aoff, asky, ahot, acold, rows, antname ) ;
3456 aoff->unsetSelection() ;
3457 ahot->unsetSelection() ;
3458 asky->unsetSelection() ;
3459 sel.reset() ;
3460 iter->next() ;
3461 }
3462 delete iter ;
3463 s->table_ = torg ;
3464 s->attach() ;
3465
3466 // flux unit
3467 out->setFluxUnit( "K" ) ;
3468
3469 return out ;
3470 }
3471}
3472
3473CountedPtr<Scantable> STMath::almacal( const CountedPtr<Scantable>& s,
3474 const String calmode )
3475{
3476 // frequency switch
3477 if ( calmode == "fs" ) {
3478 return almacalfs( s ) ;
3479 }
3480 else {
3481// double t0, t1 ;
3482// t0 = mathutil::gettimeofday_sec() ;
3483 vector<bool> masks = s->getMask( 0 ) ;
3484
3485 // save original table selection
3486 Table torg = s->table_ ;
3487
3488 // off scan
3489 // TODO 2010/01/08 TN
3490 // Grouping by time should be needed before averaging.
3491 // Each group must have own unique SCANNO (should be renumbered).
3492 // See PIPELINE/SDCalibration.py
3493 bool insitu = insitu_ ;
3494 insitu_ = false ;
3495 // share off scan before average with out
3496 CountedPtr<Scantable> out = getScantable( s, true ) ;
3497 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSOFF ) ;
3498 out->attach() ;
3499 insitu_ = insitu ;
3500 CountedPtr<Scantable> aoff = averageWithinSession( out,
3501 masks,
3502 "TINT" ) ;
3503
3504 // on scan
3505// t0 = mathutil::gettimeofday_sec() ;
3506 s->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSON ) ;
3507 s->attach() ;
3508 out->table_ = out->originalTable_ ;
3509 out->attach() ;
3510 out->table().addRow( s->nrow() ) ;
3511 copyRows( out->table(), s->table(), 0, 0, s->nrow(), False ) ;
3512// t1 = mathutil::gettimeofday_sec() ;
3513// cout << "elapsed time for preparing output table: " << t1-t0 << " sec" << endl ;
3514
3515 // process each on scan
3516// t0 = mathutil::gettimeofday_sec() ;
3517
3518 // using STIdxIterAcc
3519 vector<string> cols( 3 ) ;
3520 cols[0] = "BEAMNO" ;
3521 cols[1] = "POLNO" ;
3522 cols[2] = "IFNO" ;
3523 STIdxIter *iter = new STIdxIterAcc( out, cols ) ;
3524 STSelector sel ;
3525 while ( !iter->pastEnd() ) {
3526 Vector<uInt> ids = iter->current() ;
3527 stringstream ss ;
3528 ss << "SELECT FROM $1 WHERE "
3529 << "BEAMNO==" << ids[0] << "&&"
3530 << "POLNO==" << ids[1] << "&&"
3531 << "IFNO==" << ids[2] ;
3532 //cout << "TaQL string: " << ss.str() << endl ;
3533 sel.setTaQL( ss.str() ) ;
3534 aoff->setSelection( sel ) ;
3535 Vector<uInt> rows = iter->getRows( SHARE ) ;
3536 // out should be an exact copy of s except that SPECTRA column is empty
3537 calibrateALMA( out, s, aoff, rows ) ;
3538 aoff->unsetSelection() ;
3539 sel.reset() ;
3540 iter->next() ;
3541 }
3542 delete iter ;
3543 s->table_ = torg ;
3544 s->attach() ;
3545
3546// t1 = mathutil::gettimeofday_sec() ;
3547// cout << "elapsed time for calibration: " << t1-t0 << " sec" << endl ;
3548
3549 // flux unit
3550 out->setFluxUnit( "K" ) ;
3551
3552 return out ;
3553 }
3554}
3555
3556CountedPtr<Scantable> STMath::cwcalfs( const CountedPtr<Scantable>& s,
3557 const String antname )
3558{
3559 vector<int> types ;
3560
3561 // APEX calibration mode
3562 int apexcalmode = 1 ;
3563
3564 if ( antname.find( "APEX" ) != string::npos ) {
3565 // check if off scan exists or not
3566 STSelector sel = STSelector() ;
3567 //sel.setName( offstr1 ) ;
3568 types.push_back( SrcType::FLOOFF ) ;
3569 sel.setTypes( types ) ;
3570 try {
3571 s->setSelection( sel ) ;
3572 }
3573 catch ( AipsError &e ) {
3574 apexcalmode = 0 ;
3575 }
3576 sel.reset() ;
3577 }
3578 s->unsetSelection() ;
3579 types.clear() ;
3580
3581 vector<bool> masks = s->getMask( 0 ) ;
3582 CountedPtr<Scantable> ssig, sref ;
3583 //CountedPtr<Scantable> out ;
3584 bool insitu = insitu_ ;
3585 insitu_ = False ;
3586 CountedPtr<Scantable> out = getScantable( s, true ) ;
3587 insitu_ = insitu ;
3588
3589 if ( antname.find( "APEX" ) != string::npos ) {
3590 // APEX calibration
3591 // sky scan
3592 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOSKY ) ;
3593 out->attach() ;
3594 CountedPtr<Scantable> askylo = averageWithinSession( out,
3595 masks,
3596 "TINT" ) ;
3597 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHISKY ) ;
3598 out->attach() ;
3599 CountedPtr<Scantable> askyhi = averageWithinSession( out,
3600 masks,
3601 "TINT" ) ;
3602
3603 // hot scan
3604 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOHOT ) ;
3605 out->attach() ;
3606 CountedPtr<Scantable> ahotlo = averageWithinSession( out,
3607 masks,
3608 "TINT" ) ;
3609 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHIHOT ) ;
3610 out->attach() ;
3611 CountedPtr<Scantable> ahothi = averageWithinSession( out,
3612 masks,
3613 "TINT" ) ;
3614
3615 // cold scan
3616 CountedPtr<Scantable> acoldlo, acoldhi ;
3617// out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOCOLD ) ;
3618// out->attach() ;
3619// CountedPtr<Scantable> acoldlo = averageWithinSession( out,
3620// masks,
3621// "TINT" ) ;
3622// out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHICOLD ) ;
3623// out->attach() ;
3624// CountedPtr<Scantable> acoldhi = averageWithinSession( out,
3625// masks,
3626// "TINT" ) ;
3627
3628 // ref scan
3629 insitu_ = false ;
3630 sref = getScantable( s, true ) ;
3631 CountedPtr<Scantable> rref = getScantable( s, true ) ;
3632 insitu_ = insitu ;
3633 rref->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FSLO ) ;
3634 rref->attach() ;
3635 copyRows( sref->table_, rref->table_, 0, 0, rref->nrow(), False, True, False ) ;
3636
3637 // sig scan
3638 insitu_ = false ;
3639 ssig = getScantable( s, true ) ;
3640 CountedPtr<Scantable> rsig = getScantable( s, true ) ;
3641 insitu_ = insitu ;
3642 rsig->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FSHI ) ;
3643 rsig->attach() ;
3644 copyRows( ssig->table_, rsig->table_, 0, 0, rsig->nrow(), False, True, False ) ;
3645
3646 if ( apexcalmode == 0 ) {
3647 // using STIdxIterAcc
3648 vector<string> cols( 3 ) ;
3649 cols[0] = "BEAMNO" ;
3650 cols[1] = "POLNO" ;
3651 cols[2] = "IFNO" ;
3652 STIdxIter *iter = new STIdxIterAcc( ssig, cols ) ;
3653 STSelector sel ;
3654 vector< CountedPtr<Scantable> > on( 2 ) ;
3655 on[0] = rsig ;
3656 on[1] = rref ;
3657 vector< CountedPtr<Scantable> > sky( 2 ) ;
3658 sky[0] = askylo ;
3659 sky[1] = askyhi ;
3660 vector< CountedPtr<Scantable> > hot( 2 ) ;
3661 hot[0] = ahotlo ;
3662 hot[1] = ahothi ;
3663 vector< CountedPtr<Scantable> > cold( 2 ) ;
3664 while ( !iter->pastEnd() ) {
3665 Vector<uInt> ids = iter->current() ;
3666 stringstream ss ;
3667 ss << "SELECT FROM $1 WHERE "
3668 << "BEAMNO==" << ids[0] << "&&"
3669 << "POLNO==" << ids[1] << "&&"
3670 << "IFNO==" << ids[2] ;
3671 //cout << "TaQL string: " << ss.str() << endl ;
3672 sel.setTaQL( ss.str() ) ;
3673 sky[0]->setSelection( sel ) ;
3674 sky[1]->setSelection( sel ) ;
3675 hot[0]->setSelection( sel ) ;
3676 hot[1]->setSelection( sel ) ;
3677 Vector<uInt> rows = iter->getRows( SHARE ) ;
3678 calibrateAPEXFS( ssig, sref, on, sky, hot, cold, rows ) ;
3679 sky[0]->unsetSelection() ;
3680 sky[1]->unsetSelection() ;
3681 hot[0]->unsetSelection() ;
3682 hot[1]->unsetSelection() ;
3683 sel.reset() ;
3684 iter->next() ;
3685 }
3686 delete iter ;
3687
3688 }
3689 else if ( apexcalmode == 1 ) {
3690 // APEX fs data with off scan
3691 // off scan
3692 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FLOOFF ) ;
3693 out->attach() ;
3694 CountedPtr<Scantable> aofflo = averageWithinSession( out,
3695 masks,
3696 "TINT" ) ;
3697 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::FHIOFF ) ;
3698 out->attach() ;
3699 CountedPtr<Scantable> aoffhi = averageWithinSession( out,
3700 masks,
3701 "TINT" ) ;
3702
3703 // process each sig and ref scan
3704// STSelector sel ;
3705 vector<string> cols( 3 ) ;
3706 cols[0] = "BEAMNO" ;
3707 cols[1] = "POLNO" ;
3708 cols[2] = "IFNO" ;
3709 STIdxIter *iter = new STIdxIterAcc( ssig, cols ) ;
3710 STSelector sel ;
3711 while ( !iter->pastEnd() ) {
3712 Vector<uInt> ids = iter->current() ;
3713 stringstream ss ;
3714 ss << "SELECT FROM $1 WHERE "
3715 << "BEAMNO==" << ids[0] << "&&"
3716 << "POLNO==" << ids[1] << "&&"
3717 << "IFNO==" << ids[2] ;
3718 //cout << "TaQL string: " << ss.str() << endl ;
3719 sel.setTaQL( ss.str() ) ;
3720 aofflo->setSelection( sel ) ;
3721 ahotlo->setSelection( sel ) ;
3722 askylo->setSelection( sel ) ;
3723 Vector<uInt> rows = iter->getRows( SHARE ) ;
3724 calibrateCW( ssig, rsig, aofflo, askylo, ahotlo, acoldlo, rows, antname ) ;
3725 aofflo->unsetSelection() ;
3726 ahotlo->unsetSelection() ;
3727 askylo->unsetSelection() ;
3728 sel.reset() ;
3729 iter->next() ;
3730 }
3731 delete iter ;
3732 iter = new STIdxIterAcc( sref, cols ) ;
3733 while ( !iter->pastEnd() ) {
3734 Vector<uInt> ids = iter->current() ;
3735 stringstream ss ;
3736 ss << "SELECT FROM $1 WHERE "
3737 << "BEAMNO==" << ids[0] << "&&"
3738 << "POLNO==" << ids[1] << "&&"
3739 << "IFNO==" << ids[2] ;
3740 //cout << "TaQL string: " << ss.str() << endl ;
3741 sel.setTaQL( ss.str() ) ;
3742 aoffhi->setSelection( sel ) ;
3743 ahothi->setSelection( sel ) ;
3744 askyhi->setSelection( sel ) ;
3745 Vector<uInt> rows = iter->getRows( SHARE ) ;
3746 calibrateCW( sref, rref, aoffhi, askyhi, ahothi, acoldhi, rows, antname ) ;
3747 aoffhi->unsetSelection() ;
3748 ahothi->unsetSelection() ;
3749 askyhi->unsetSelection() ;
3750 sel.reset() ;
3751 iter->next() ;
3752 }
3753 delete iter ;
3754 }
3755 }
3756 else {
3757 // non-APEX fs data
3758 // sky scan
3759 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::SKY ) ;
3760 out->attach() ;
3761 CountedPtr<Scantable> asky = averageWithinSession( out,
3762 masks,
3763 "TINT" ) ;
3764 STSelector sel = STSelector() ;
3765
3766 // hot scan
3767 out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::HOT ) ;
3768 out->attach() ;
3769 CountedPtr<Scantable> ahot = averageWithinSession( out,
3770 masks,
3771 "TINT" ) ;
3772
3773 // cold scan
3774 CountedPtr<Scantable> acold ;
3775// out->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::COLD ) ;
3776// out->attach() ;
3777// CountedPtr<Scantable> acold = averageWithinSession( out,
3778// masks,
3779// "TINT" ) ;
3780
3781 // ref scan
3782 bool insitu = insitu_ ;
3783 insitu_ = false ;
3784 sref = getScantable( s, true ) ;
3785 CountedPtr<Scantable> rref = getScantable( s, true ) ;
3786 insitu_ = insitu ;
3787 rref->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSOFF ) ;
3788 rref->attach() ;
3789 copyRows( sref->table_, rref->table_, 0, 0, rref->nrow(), False, True, False ) ;
3790
3791 // sig scan
3792 insitu_ = false ;
3793 ssig = getScantable( s, true ) ;
3794 CountedPtr<Scantable> rsig = getScantable( s, true ) ;
3795 insitu_ = insitu ;
3796 rsig->table_ = s->table_( s->table_.col("SRCTYPE") == (Int)SrcType::PSON ) ;
3797 rsig->attach() ;
3798 copyRows( ssig->table_, rsig->table_, 0, 0, rsig->nrow(), False, True, False ) ;
3799
3800 // process each sig and ref scan
3801 vector<string> cols( 3 ) ;
3802 cols[0] = "BEAMNO" ;
3803 cols[1] = "POLNO" ;
3804 cols[2] = "IFNO" ;
3805 STIdxIter *iter = new STIdxIterAcc( ssig, cols ) ;
3806 while ( !iter->pastEnd() ) {
3807 Vector<uInt> ids = iter->current() ;
3808 stringstream ss ;
3809 ss << "SELECT FROM $1 WHERE "
3810 << "BEAMNO==" << ids[0] << "&&"
3811 << "POLNO==" << ids[1] << "&&"
3812 << "IFNO==" << ids[2] ;
3813 //cout << "TaQL string: " << ss.str() << endl ;
3814 sel.setTaQL( ss.str() ) ;
3815 ahot->setSelection( sel ) ;
3816 asky->setSelection( sel ) ;
3817 Vector<uInt> rows = iter->getRows( SHARE ) ;
3818 // out should be an exact copy of s except that SPECTRA column is empty
3819 calibrateFS( ssig, sref, rsig, rref, asky, ahot, acold, rows ) ;
3820 ahot->unsetSelection() ;
3821 asky->unsetSelection() ;
3822 sel.reset() ;
3823 iter->next() ;
3824 }
3825 delete iter ;
3826 }
3827
3828 // do folding if necessary
3829 Table sigtab = ssig->table() ;
3830 Table reftab = sref->table() ;
3831 ScalarColumn<uInt> reffidCol ;
3832 Int nchan = (Int)ssig->nchan() ;
3833 reffidCol.attach( reftab, "FREQ_ID" ) ;
3834 Vector<uInt> sfids = ssig->mfreqidCol_.getColumn() ;
3835 Vector<uInt> rfids = sref->mfreqidCol_.getColumn() ;
3836 vector<uInt> sfids_unique ;
3837 vector<uInt> rfids_unique ;
3838 vector<uInt> sifno_unique ;
3839 vector<uInt> rifno_unique ;
3840 for ( uInt i = 0 ; i < sfids.nelements() ; i++ ) {
3841 if ( count( sfids_unique.begin(), sfids_unique.end(), sfids[i] ) == 0 ) {
3842 sfids_unique.push_back( sfids[i] ) ;
3843 sifno_unique.push_back( ssig->getIF( i ) ) ;
3844 }
3845 if ( count( rfids_unique.begin(), rfids_unique.end(), rfids[i] ) == 0 ) {
3846 rfids_unique.push_back( rfids[i] ) ;
3847 rifno_unique.push_back( sref->getIF( i ) ) ;
3848 }
3849 }
3850 double refpix_sig, refval_sig, increment_sig ;
3851 double refpix_ref, refval_ref, increment_ref ;
3852 vector< CountedPtr<Scantable> > tmp( sfids_unique.size() ) ;
3853 for ( uInt i = 0 ; i < sfids_unique.size() ; i++ ) {
3854 ssig->frequencies().getEntry( refpix_sig, refval_sig, increment_sig, sfids_unique[i] ) ;
3855 sref->frequencies().getEntry( refpix_ref, refval_ref, increment_ref, rfids_unique[i] ) ;
3856 if ( refpix_sig == refpix_ref ) {
3857 double foffset = refval_ref - refval_sig ;
3858 int choffset = static_cast<int>(foffset/increment_sig) ;
3859 double doffset = foffset / increment_sig ;
3860 if ( abs(choffset) >= nchan ) {
3861 LogIO os( LogOrigin( "STMath", "cwcalfs", WHERE ) ) ;
3862 os << "FREQ_ID=[" << sfids_unique[i] << "," << rfids_unique[i] << "]: out-band frequency switching, no folding" << LogIO::POST ;
3863 os << "Just return signal data" << LogIO::POST ;
3864 //std::vector< CountedPtr<Scantable> > tabs ;
3865 //tabs.push_back( ssig ) ;
3866 //tabs.push_back( sref ) ;
3867 //out = merge( tabs ) ;
3868 tmp[i] = ssig ;
3869 }
3870 else {
3871 STSelector sel = STSelector() ;
3872 vector<int> v( 1, sifno_unique[i] ) ;
3873 sel.setIFs( v ) ;
3874 ssig->setSelection( sel ) ;
3875 sel.reset() ;
3876 v[0] = rifno_unique[i] ;
3877 sel.setIFs( v ) ;
3878 sref->setSelection( sel ) ;
3879 sel.reset() ;
3880 if ( antname.find( "APEX" ) != string::npos ) {
3881 tmp[i] = dofold( ssig, sref, 0.5*doffset, -0.5*doffset ) ;
3882 //tmp[i] = dofold( ssig, sref, doffset ) ;
3883 }
3884 else {
3885 tmp[i] = dofold( ssig, sref, doffset ) ;
3886 }
3887 ssig->unsetSelection() ;
3888 sref->unsetSelection() ;
3889 }
3890 }
3891 }
3892
3893 if ( tmp.size() > 1 ) {
3894 out = merge( tmp ) ;
3895 }
3896 else {
3897 out = tmp[0] ;
3898 }
3899
3900 // flux unit
3901 out->setFluxUnit( "K" ) ;
3902
3903 return out ;
3904}
3905
3906CountedPtr<Scantable> STMath::almacalfs( const CountedPtr<Scantable>& s )
3907{
3908 (void) s; //currently unused
3909 CountedPtr<Scantable> out ;
3910
3911 return out ;
3912}
3913
3914Vector<Float> STMath::getSpectrumFromTime( double reftime,
3915 const Vector<Double> &timeVec,
3916 const vector<int> &idx,
3917 const Matrix<Float>& spectra,
3918 string mode )
3919{
3920 LogIO os( LogOrigin( "STMath", "getSpectrumFromTime", WHERE ) ) ;
3921 Vector<Float> sp ;
3922 uInt ncol = spectra.ncolumn() ;
3923
3924 if ( ncol == 0 ) {
3925 os << LogIO::SEVERE << "No spectra in the input scantable. Return empty spectrum." << LogIO::POST ;
3926 return sp ;
3927 }
3928 else if ( ncol == 1 ) {
3929 //os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ;
3930 sp.reference( spectra.column( 0 ) ) ;
3931 return sp ;
3932 }
3933 else {
3934 if ( mode == "before" ) {
3935 int id = -1 ;
3936 if ( idx[0] != -1 ) {
3937 id = idx[0] ;
3938 }
3939 else if ( idx[1] != -1 ) {
3940 os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
3941 id = idx[1] ;
3942 }
3943 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
3944 sp.reference( spectra.column( id ) ) ;
3945 }
3946 else if ( mode == "after" ) {
3947 int id = -1 ;
3948 if ( idx[1] != -1 ) {
3949 id = idx[1] ;
3950 }
3951 else if ( idx[0] != -1 ) {
3952 os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
3953 id = idx[1] ;
3954 }
3955 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
3956 sp.reference( spectra.column( id ) ) ;
3957 }
3958 else if ( mode == "nearest" ) {
3959 int id = -1 ;
3960 if ( idx[0] == -1 ) {
3961 id = idx[1] ;
3962 }
3963 else if ( idx[1] == -1 ) {
3964 id = idx[0] ;
3965 }
3966 else if ( idx[0] == idx[1] ) {
3967 id = idx[0] ;
3968 }
3969 else {
3970 double t0 = timeVec[idx[0]] ;
3971 double t1 = timeVec[idx[1]] ;
3972 double tref = reftime ;
3973 if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
3974 id = idx[1] ;
3975 }
3976 else {
3977 id = idx[0] ;
3978 }
3979 }
3980 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
3981 sp.reference( spectra.column( id ) ) ;
3982 }
3983 else if ( mode == "linear" ) {
3984 if ( idx[0] == -1 ) {
3985 // use after
3986 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
3987 int id = idx[1] ;
3988 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
3989 sp.reference( spectra.column( id ) ) ;
3990 }
3991 else if ( idx[1] == -1 ) {
3992 // use before
3993 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
3994 int id = idx[0] ;
3995 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
3996 sp.reference( spectra.column( id ) ) ;
3997 }
3998 else if ( idx[0] == idx[1] ) {
3999 // use before
4000 //os << "No need to interporate." << LogIO::POST ;
4001 int id = idx[0] ;
4002 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4003 sp.reference( spectra.column( id ) ) ;
4004 }
4005 else {
4006 // do interpolation
4007 //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
4008 double t0 = timeVec[idx[0]] ;
4009 double t1 = timeVec[idx[1]] ;
4010 double tref = reftime ;
4011 sp = spectra.column( idx[0] ).copy() ;
4012 Vector<Float> sp1( spectra.column( idx[1] ) ) ;
4013 double tfactor = ( tref - t0 ) / ( t1 - t0 ) ;
4014 for ( unsigned int i = 0 ; i < sp.size() ; i++ ) {
4015 sp[i] = ( sp1[i] - sp[i] ) * tfactor + sp[i] ;
4016 }
4017 }
4018 }
4019 else {
4020 os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
4021 }
4022 return sp ;
4023 }
4024}
4025
4026vector<int> STMath::getRowIdFromTime( double reftime, const Vector<Double> &t )
4027{
4028// double reft = reftime ;
4029 double dtmin = 1.0e100 ;
4030 double dtmax = -1.0e100 ;
4031// vector<double> dt ;
4032 int just_before = -1 ;
4033 int just_after = -1 ;
4034 Vector<Double> dt = t - reftime ;
4035 for ( unsigned int i = 0 ; i < dt.size() ; i++ ) {
4036 if ( dt[i] > 0.0 ) {
4037 // after reftime
4038 if ( dt[i] < dtmin ) {
4039 just_after = i ;
4040 dtmin = dt[i] ;
4041 }
4042 }
4043 else if ( dt[i] < 0.0 ) {
4044 // before reftime
4045 if ( dt[i] > dtmax ) {
4046 just_before = i ;
4047 dtmax = dt[i] ;
4048 }
4049 }
4050 else {
4051 // just a reftime
4052 just_before = i ;
4053 just_after = i ;
4054 dtmax = 0 ;
4055 dtmin = 0 ;
4056 break ;
4057 }
4058 }
4059
4060 vector<int> v(2) ;
4061 v[0] = just_before ;
4062 v[1] = just_after ;
4063
4064 return v ;
4065}
4066
4067Vector<Float> STMath::getTcalFromTime( double reftime,
4068 const Vector<Double> &timeVec,
4069 const vector<int> &idx,
4070 const CountedPtr<Scantable>& s,
4071 string mode )
4072{
4073 LogIO os( LogOrigin( "STMath", "getTcalFromTime", WHERE ) ) ;
4074 STTcal tcalTable = s->tcal() ;
4075 String time ;
4076 Vector<Float> tcalval ;
4077 if ( s->nrow() == 0 ) {
4078 os << LogIO::SEVERE << "No row in the input scantable. Return empty tcal." << LogIO::POST ;
4079 return tcalval ;
4080 }
4081 else if ( s->nrow() == 1 ) {
4082 uInt tcalid = s->getTcalId( 0 ) ;
4083 //os << "use row " << 0 << " (tcalid = " << tcalid << ")" << LogIO::POST ;
4084 tcalTable.getEntry( time, tcalval, tcalid ) ;
4085 return tcalval ;
4086 }
4087 else {
4088 if ( mode == "before" ) {
4089 int id = -1 ;
4090 if ( idx[0] != -1 ) {
4091 id = idx[0] ;
4092 }
4093 else if ( idx[1] != -1 ) {
4094 os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
4095 id = idx[1] ;
4096 }
4097 uInt tcalid = s->getTcalId( id ) ;
4098 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
4099 tcalTable.getEntry( time, tcalval, tcalid ) ;
4100 }
4101 else if ( mode == "after" ) {
4102 int id = -1 ;
4103 if ( idx[1] != -1 ) {
4104 id = idx[1] ;
4105 }
4106 else if ( idx[0] != -1 ) {
4107 os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
4108 id = idx[1] ;
4109 }
4110 uInt tcalid = s->getTcalId( id ) ;
4111 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
4112 tcalTable.getEntry( time, tcalval, tcalid ) ;
4113 }
4114 else if ( mode == "nearest" ) {
4115 int id = -1 ;
4116 if ( idx[0] == -1 ) {
4117 id = idx[1] ;
4118 }
4119 else if ( idx[1] == -1 ) {
4120 id = idx[0] ;
4121 }
4122 else if ( idx[0] == idx[1] ) {
4123 id = idx[0] ;
4124 }
4125 else {
4126 double t0 = timeVec[idx[0]] ;
4127 double t1 = timeVec[idx[1]] ;
4128 if ( abs( t0 - reftime ) > abs( t1 - reftime ) ) {
4129 id = idx[1] ;
4130 }
4131 else {
4132 id = idx[0] ;
4133 }
4134 }
4135 uInt tcalid = s->getTcalId( id ) ;
4136 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
4137 tcalTable.getEntry( time, tcalval, tcalid ) ;
4138 }
4139 else if ( mode == "linear" ) {
4140 if ( idx[0] == -1 ) {
4141 // use after
4142 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
4143 int id = idx[1] ;
4144 uInt tcalid = s->getTcalId( id ) ;
4145 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
4146 tcalTable.getEntry( time, tcalval, tcalid ) ;
4147 }
4148 else if ( idx[1] == -1 ) {
4149 // use before
4150 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
4151 int id = idx[0] ;
4152 uInt tcalid = s->getTcalId( id ) ;
4153 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
4154 tcalTable.getEntry( time, tcalval, tcalid ) ;
4155 }
4156 else if ( idx[0] == idx[1] ) {
4157 // use before
4158 //os << "No need to interporate." << LogIO::POST ;
4159 int id = idx[0] ;
4160 uInt tcalid = s->getTcalId( id ) ;
4161 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
4162 tcalTable.getEntry( time, tcalval, tcalid ) ;
4163 }
4164 else {
4165 // do interpolation
4166 //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
4167 double t0 = timeVec[idx[0]] ;
4168 double t1 = timeVec[idx[1]] ;
4169 Vector<Float> tcal0 ;
4170 uInt tcalid0 = s->getTcalId( idx[0] ) ;
4171 uInt tcalid1 = s->getTcalId( idx[1] ) ;
4172 tcalTable.getEntry( time, tcal0, tcalid0 ) ;
4173 tcalTable.getEntry( time, tcalval, tcalid1 ) ;
4174 double tfactor = (reftime - t0) / (t1 - t0) ;
4175 for ( unsigned int i = 0 ; i < tcal0.size() ; i++ ) {
4176 tcalval[i] = ( tcalval[i] - tcal0[i] ) * tfactor + tcal0[i] ;
4177 }
4178 }
4179 }
4180 else {
4181 os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
4182 }
4183 return tcalval ;
4184 }
4185}
4186
4187Vector<Float> STMath::getTsysFromTime( double reftime,
4188 const Vector<Double> &timeVec,
4189 const vector<int> &idx,
4190 const CountedPtr<Scantable> &s,
4191 string mode )
4192{
4193 LogIO os( LogOrigin( "STMath", "getTsysFromTime", WHERE ) ) ;
4194 ArrayColumn<Float> tsysCol ;
4195 tsysCol.attach( s->table(), "TSYS" ) ;
4196 Vector<Float> tsysval ;
4197 if ( s->nrow() == 0 ) {
4198 os << LogIO::SEVERE << "No row in the input scantable. Return empty tsys." << LogIO::POST ;
4199 return tsysval ;
4200 }
4201 else if ( s->nrow() == 1 ) {
4202 //os << "use row " << 0 << LogIO::POST ;
4203 tsysval = tsysCol( 0 ) ;
4204 return tsysval ;
4205 }
4206 else {
4207 if ( mode == "before" ) {
4208 int id = -1 ;
4209 if ( idx[0] != -1 ) {
4210 id = idx[0] ;
4211 }
4212 else if ( idx[1] != -1 ) {
4213 os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
4214 id = idx[1] ;
4215 }
4216 //os << "use row " << id << LogIO::POST ;
4217 tsysval = tsysCol( id ) ;
4218 }
4219 else if ( mode == "after" ) {
4220 int id = -1 ;
4221 if ( idx[1] != -1 ) {
4222 id = idx[1] ;
4223 }
4224 else if ( idx[0] != -1 ) {
4225 os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
4226 id = idx[1] ;
4227 }
4228 //os << "use row " << id << LogIO::POST ;
4229 tsysval = tsysCol( id ) ;
4230 }
4231 else if ( mode == "nearest" ) {
4232 int id = -1 ;
4233 if ( idx[0] == -1 ) {
4234 id = idx[1] ;
4235 }
4236 else if ( idx[1] == -1 ) {
4237 id = idx[0] ;
4238 }
4239 else if ( idx[0] == idx[1] ) {
4240 id = idx[0] ;
4241 }
4242 else {
4243 double t0 = timeVec[idx[0]] ;
4244 double t1 = timeVec[idx[1]] ;
4245 if ( abs( t0 - reftime ) > abs( t1 - reftime ) ) {
4246 id = idx[1] ;
4247 }
4248 else {
4249 id = idx[0] ;
4250 }
4251 }
4252 //os << "use row " << id << LogIO::POST ;
4253 tsysval = tsysCol( id ) ;
4254 }
4255 else if ( mode == "linear" ) {
4256 if ( idx[0] == -1 ) {
4257 // use after
4258 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
4259 int id = idx[1] ;
4260 //os << "use row " << id << LogIO::POST ;
4261 tsysval = tsysCol( id ) ;
4262 }
4263 else if ( idx[1] == -1 ) {
4264 // use before
4265 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
4266 int id = idx[0] ;
4267 //os << "use row " << id << LogIO::POST ;
4268 tsysval = tsysCol( id ) ;
4269 }
4270 else if ( idx[0] == idx[1] ) {
4271 // use before
4272 //os << "No need to interporate." << LogIO::POST ;
4273 int id = idx[0] ;
4274 //os << "use row " << id << LogIO::POST ;
4275 tsysval = tsysCol( id ) ;
4276 }
4277 else {
4278 // do interpolation
4279 //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
4280 double t0 = timeVec[idx[0]] ;
4281 double t1 = timeVec[idx[1]] ;
4282 Vector<Float> tsys0 ;
4283 tsys0 = tsysCol( idx[0] ) ;
4284 tsysval = tsysCol( idx[1] ) ;
4285 double tfactor = (reftime - t0) / (t1 - t0) ;
4286 for ( unsigned int i = 0 ; i < tsys0.size() ; i++ ) {
4287 tsysval[i] = ( tsysval[i] - tsys0[i] ) * tfactor + tsys0[i] ;
4288 }
4289 }
4290 }
4291 else {
4292 os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
4293 }
4294 return tsysval ;
4295 }
4296}
4297
4298void STMath::calibrateCW( CountedPtr<Scantable> &out,
4299 const CountedPtr<Scantable>& on,
4300 const CountedPtr<Scantable>& off,
4301 const CountedPtr<Scantable>& sky,
4302 const CountedPtr<Scantable>& hot,
4303 const CountedPtr<Scantable>& cold,
4304 const Vector<uInt> &rows,
4305 const String &antname )
4306{
4307 // 2012/05/22 TN
4308 // Assume that out has empty SPECTRA column
4309
4310 // if rows is empty, just return
4311 if ( rows.nelements() == 0 )
4312 return ;
4313 ROScalarColumn<Double> timeCol( off->table(), "TIME" ) ;
4314 Vector<Double> timeOff = timeCol.getColumn() ;
4315 timeCol.attach( sky->table(), "TIME" ) ;
4316 Vector<Double> timeSky = timeCol.getColumn() ;
4317 timeCol.attach( hot->table(), "TIME" ) ;
4318 Vector<Double> timeHot = timeCol.getColumn() ;
4319 timeCol.attach( on->table(), "TIME" ) ;
4320 ROArrayColumn<Float> arrayFloatCol( off->table(), "SPECTRA" ) ;
4321 Matrix<Float> offspectra = arrayFloatCol.getColumn() ;
4322 arrayFloatCol.attach( sky->table(), "SPECTRA" ) ;
4323 Matrix<Float> skyspectra = arrayFloatCol.getColumn() ;
4324 arrayFloatCol.attach( hot->table(), "SPECTRA" ) ;
4325 Matrix<Float> hotspectra = arrayFloatCol.getColumn() ;
4326 unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ;
4327 // I know that the data is contiguous
4328 const uInt *p = rows.data() ;
4329 vector<int> ids( 2 ) ;
4330 Block<uInt> flagchan( spsize ) ;
4331 uInt nflag = 0 ;
4332 for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
4333 double reftime = timeCol.asdouble(*p) ;
4334 ids = getRowIdFromTime( reftime, timeOff ) ;
4335 Vector<Float> spoff = getSpectrumFromTime( reftime, timeOff, ids, offspectra, "linear" ) ;
4336 ids = getRowIdFromTime( reftime, timeSky ) ;
4337 Vector<Float> spsky = getSpectrumFromTime( reftime, timeSky, ids, skyspectra, "linear" ) ;
4338 Vector<Float> tcal = getTcalFromTime( reftime, timeSky, ids, sky, "linear" ) ;
4339 Vector<Float> tsys = getTsysFromTime( reftime, timeSky, ids, sky, "linear" ) ;
4340 ids = getRowIdFromTime( reftime, timeHot ) ;
4341 Vector<Float> sphot = getSpectrumFromTime( reftime, timeHot, ids, hotspectra, "linear" ) ;
4342 Vector<Float> spec = on->specCol_( *p ) ;
4343 if ( antname.find( "APEX" ) != String::npos ) {
4344 // using gain array
4345 for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
4346 if ( spoff[j] == 0.0 || (sphot[j]-spsky[j]) == 0.0 ) {
4347 spec[j] = 0.0 ;
4348 flagchan[nflag++] = j ;
4349 }
4350 else {
4351 spec[j] = ( ( spec[j] - spoff[j] ) / spoff[j] )
4352 * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
4353 }
4354 }
4355 }
4356 else {
4357 // Chopper-Wheel calibration (Ulich & Haas 1976)
4358 for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
4359 if ( (sphot[j]-spsky[j]) == 0.0 ) {
4360 spec[j] = 0.0 ;
4361 flagchan[nflag++] = j ;
4362 }
4363 else {
4364 spec[j] = ( spec[j] - spoff[j] ) / ( sphot[j] - spsky[j] ) * tcal[j] ;
4365 }
4366 }
4367 }
4368 out->specCol_.put( *p, spec ) ;
4369 out->tsysCol_.put( *p, tsys ) ;
4370 if ( nflag > 0 ) {
4371 Vector<uChar> fl = out->flagsCol_( *p ) ;
4372 for ( unsigned int j = 0 ; j < nflag ; j++ ) {
4373 fl[flagchan[j]] = (uChar)True ;
4374 }
4375 out->flagsCol_.put( *p, fl ) ;
4376 }
4377 nflag = 0 ;
4378 p++ ;
4379 }
4380}
4381
4382void STMath::calibrateALMA( CountedPtr<Scantable>& out,
4383 const CountedPtr<Scantable>& on,
4384 const CountedPtr<Scantable>& off,
4385 const Vector<uInt>& rows )
4386{
4387 // 2012/05/22 TN
4388 // Assume that out has empty SPECTRA column
4389
4390 // if rows is empty, just return
4391 if ( rows.nelements() == 0 )
4392 return ;
4393 ROScalarColumn<Double> timeCol( off->table(), "TIME" ) ;
4394 Vector<Double> timeVec = timeCol.getColumn() ;
4395 timeCol.attach( on->table(), "TIME" ) ;
4396 ROArrayColumn<Float> arrayFloatCol( off->table(), "SPECTRA" ) ;
4397 Matrix<Float> offspectra = arrayFloatCol.getColumn() ;
4398 unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ;
4399 // I know that the data is contiguous
4400 const uInt *p = rows.data() ;
4401 vector<int> ids( 2 ) ;
4402 Block<uInt> flagchan( spsize ) ;
4403 uInt nflag = 0 ;
4404 for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
4405 double reftime = timeCol.asdouble(*p) ;
4406 ids = getRowIdFromTime( reftime, timeVec ) ;
4407 Vector<Float> spoff = getSpectrumFromTime( reftime, timeVec, ids, offspectra, "linear" ) ;
4408 //Vector<Float> spoff = getSpectrumFromTime( reftime, timeVec, off, "linear" ) ;
4409 Vector<Float> spec = on->specCol_( *p ) ;
4410 Vector<Float> tsys = on->tsysCol_( *p ) ;
4411 // ALMA Calibration
4412 //
4413 // Ta* = Tsys * ( ON - OFF ) / OFF
4414 //
4415 // 2010/01/07 Takeshi Nakazato
4416 unsigned int tsyssize = tsys.nelements() ;
4417 for ( unsigned int j = 0 ; j < spsize ; j++ ) {
4418 if ( spoff[j] == 0.0 ) {
4419 spec[j] = 0.0 ;
4420 flagchan[nflag++] = j ;
4421 }
4422 else {
4423 spec[j] = ( spec[j] - spoff[j] ) / spoff[j] ;
4424 }
4425 if ( tsyssize == spsize )
4426 spec[j] *= tsys[j] ;
4427 else
4428 spec[j] *= tsys[0] ;
4429 }
4430 out->specCol_.put( *p, spec ) ;
4431 if ( nflag > 0 ) {
4432 Vector<uChar> fl = out->flagsCol_( *p ) ;
4433 for ( unsigned int j = 0 ; j < nflag ; j++ ) {
4434 fl[flagchan[j]] = (uChar)True ;
4435 }
4436 out->flagsCol_.put( *p, fl ) ;
4437 }
4438 nflag = 0 ;
4439 p++ ;
4440 }
4441}
4442
4443void STMath::calibrateAPEXFS( CountedPtr<Scantable> &sig,
4444 CountedPtr<Scantable> &ref,
4445 const vector< CountedPtr<Scantable> >& on,
4446 const vector< CountedPtr<Scantable> >& sky,
4447 const vector< CountedPtr<Scantable> >& hot,
4448 const vector< CountedPtr<Scantable> >& cold,
4449 const Vector<uInt> &rows )
4450{
4451 // if rows is empty, just return
4452 if ( rows.nelements() == 0 )
4453 return ;
4454 ROScalarColumn<Double> timeCol( sky[0]->table(), "TIME" ) ;
4455 Vector<Double> timeSkyS = timeCol.getColumn() ;
4456 timeCol.attach( sky[1]->table(), "TIME" ) ;
4457 Vector<Double> timeSkyR = timeCol.getColumn() ;
4458 timeCol.attach( hot[0]->table(), "TIME" ) ;
4459 Vector<Double> timeHotS = timeCol.getColumn() ;
4460 timeCol.attach( hot[1]->table(), "TIME" ) ;
4461 Vector<Double> timeHotR = timeCol.getColumn() ;
4462 timeCol.attach( sig->table(), "TIME" ) ;
4463 ROScalarColumn<Double> timeCol2( ref->table(), "TIME" ) ;
4464 ROArrayColumn<Float> arrayFloatCol( sky[0]->table(), "SPECTRA" ) ;
4465 Matrix<Float> skyspectraS = arrayFloatCol.getColumn() ;
4466 arrayFloatCol.attach( sky[1]->table(), "SPECTRA" ) ;
4467 Matrix<Float> skyspectraR = arrayFloatCol.getColumn() ;
4468 arrayFloatCol.attach( hot[0]->table(), "SPECTRA" ) ;
4469 Matrix<Float> hotspectraS = arrayFloatCol.getColumn() ;
4470 arrayFloatCol.attach( hot[1]->table(), "SPECTRA" ) ;
4471 Matrix<Float> hotspectraR = arrayFloatCol.getColumn() ;
4472 unsigned int spsize = sig->nchan( sig->getIF(rows[0]) ) ;
4473 Vector<Float> spec( spsize ) ;
4474 // I know that the data is contiguous
4475 const uInt *p = rows.data() ;
4476 vector<int> ids( 2 ) ;
4477 Block<uInt> flagchan( spsize ) ;
4478 uInt nflag = 0 ;
4479 for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
4480 double reftime = timeCol.asdouble(*p) ;
4481 ids = getRowIdFromTime( reftime, timeSkyS ) ;
4482 Vector<Float> spskyS = getSpectrumFromTime( reftime, timeSkyS, ids, skyspectraS, "linear" ) ;
4483 Vector<Float> tcalS = getTcalFromTime( reftime, timeSkyS, ids, sky[0], "linear" ) ;
4484 Vector<Float> tsysS = getTsysFromTime( reftime, timeSkyS, ids, sky[0], "linear" ) ;
4485 ids = getRowIdFromTime( reftime, timeHotS ) ;
4486 Vector<Float> sphotS = getSpectrumFromTime( reftime, timeHotS, ids, hotspectraS ) ;
4487 reftime = timeCol2.asdouble(*p) ;
4488 ids = getRowIdFromTime( reftime, timeSkyR ) ;
4489 Vector<Float> spskyR = getSpectrumFromTime( reftime, timeSkyR, ids, skyspectraR, "linear" ) ;
4490 Vector<Float> tcalR = getTcalFromTime( reftime, timeSkyR, ids, sky[1], "linear" ) ;
4491 Vector<Float> tsysR = getTsysFromTime( reftime, timeSkyR, ids, sky[1], "linear" ) ;
4492 ids = getRowIdFromTime( reftime, timeHotR ) ;
4493 Vector<Float> sphotR = getSpectrumFromTime( reftime, timeHotR, ids, hotspectraR ) ;
4494 Vector<Float> spsig = on[0]->specCol_( *p ) ;
4495 Vector<Float> spref = on[1]->specCol_( *p ) ;
4496 for ( unsigned int j = 0 ; j < spsize ; j++ ) {
4497 if ( (sphotS[j]-spskyS[j]) == 0.0 || (sphotR[j]-spskyR[j]) == 0.0 ) {
4498 spec[j] = 0.0 ;
4499 flagchan[nflag++] = j ;
4500 }
4501 else {
4502 spec[j] = tcalS[j] * spsig[j] / ( sphotS[j] - spskyS[j] )
4503 - tcalR[j] * spref[j] / ( sphotR[j] - spskyR[j] ) ;
4504 }
4505 }
4506 sig->specCol_.put( *p, spec ) ;
4507 sig->tsysCol_.put( *p, tsysS ) ;
4508 spec *= (Float)-1.0 ;
4509 ref->specCol_.put( *p, spec ) ;
4510 ref->tsysCol_.put( *p, tsysR ) ;
4511 if ( nflag > 0 ) {
4512 Vector<uChar> flsig = sig->flagsCol_( *p ) ;
4513 Vector<uChar> flref = ref->flagsCol_( *p ) ;
4514 for ( unsigned int j = 0 ; j < nflag ; j++ ) {
4515 flsig[flagchan[j]] = (uChar)True ;
4516 flref[flagchan[j]] = (uChar)True ;
4517 }
4518 sig->flagsCol_.put( *p, flsig ) ;
4519 ref->flagsCol_.put( *p, flref ) ;
4520 }
4521 nflag = 0 ;
4522 p++ ;
4523 }
4524}
4525
4526void STMath::calibrateFS( CountedPtr<Scantable> &sig,
4527 CountedPtr<Scantable> &ref,
4528 const CountedPtr<Scantable>& rsig,
4529 const CountedPtr<Scantable>& rref,
4530 const CountedPtr<Scantable>& sky,
4531 const CountedPtr<Scantable>& hot,
4532 const CountedPtr<Scantable>& cold,
4533 const Vector<uInt> &rows )
4534{
4535 // if rows is empty, just return
4536 if ( rows.nelements() == 0 )
4537 return ;
4538 ROScalarColumn<Double> timeCol( sky->table(), "TIME" ) ;
4539 Vector<Double> timeSky = timeCol.getColumn() ;
4540 timeCol.attach( hot->table(), "TIME" ) ;
4541 Vector<Double> timeHot = timeCol.getColumn() ;
4542 timeCol.attach( sig->table(), "TIME" ) ;
4543 ROScalarColumn<Double> timeCol2( ref->table(), "TIME" ) ;
4544 ROArrayColumn<Float> arrayFloatCol( sky->table(), "SPECTRA" ) ;
4545 Matrix<Float> skyspectra = arrayFloatCol.getColumn() ;
4546 arrayFloatCol.attach( hot->table(), "SPECTRA" ) ;
4547 Matrix<Float> hotspectra = arrayFloatCol.getColumn() ;
4548 unsigned int spsize = sig->nchan( sig->getIF(rows[0]) ) ;
4549 Vector<Float> spec( spsize ) ;
4550 // I know that the data is contiguous
4551 const uInt *p = rows.data() ;
4552 vector<int> ids( 2 ) ;
4553 Block<uInt> flagchan( spsize ) ;
4554 uInt nflag = 0 ;
4555 for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
4556 double reftime = timeCol.asdouble(*p) ;
4557 ids = getRowIdFromTime( reftime, timeSky ) ;
4558 Vector<Float> spsky = getSpectrumFromTime( reftime, timeSky, ids, skyspectra, "linear" ) ;
4559 Vector<Float> tcal = getTcalFromTime( reftime, timeSky, ids, sky, "linear" ) ;
4560 Vector<Float> tsys = getTsysFromTime( reftime, timeSky, ids, sky, "linear" ) ;
4561 ids = getRowIdFromTime( reftime, timeHot ) ;
4562 Vector<Float> sphot = getSpectrumFromTime( reftime, timeHot, ids, hotspectra ) ;
4563 Vector<Float> spsig = rsig->specCol_( *p ) ;
4564 Vector<Float> spref = rref->specCol_( *p ) ;
4565 // using gain array
4566 for ( unsigned int j = 0 ; j < spsize ; j++ ) {
4567 if ( spref[j] == 0.0 || (sphot[j]-spsky[j]) == 0.0 ) {
4568 spec[j] = 0.0 ;
4569 flagchan[nflag++] = j ;
4570 }
4571 else {
4572 spec[j] = ( ( spsig[j] - spref[j] ) / spref[j] )
4573 * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
4574 }
4575 }
4576 sig->specCol_.put( *p, spec ) ;
4577 sig->tsysCol_.put( *p, tsys ) ;
4578 if ( nflag > 0 ) {
4579 Vector<uChar> fl = sig->flagsCol_( *p ) ;
4580 for ( unsigned int j = 0 ; j < nflag ; j++ ) {
4581 fl[flagchan[j]] = (uChar)True ;
4582 }
4583 sig->flagsCol_.put( *p, fl ) ;
4584 }
4585 nflag = 0 ;
4586
4587 reftime = timeCol2.asdouble(*p) ;
4588 spsky = getSpectrumFromTime( reftime, timeSky, ids, skyspectra, "linear" ) ;
4589 tcal = getTcalFromTime( reftime, timeSky, ids, sky, "linear" ) ;
4590 tsys = getTsysFromTime( reftime, timeSky, ids, sky, "linear" ) ;
4591 ids = getRowIdFromTime( reftime, timeHot ) ;
4592 sphot = getSpectrumFromTime( reftime, timeHot, ids, hotspectra ) ;
4593 // using gain array
4594 for ( unsigned int j = 0 ; j < spsize ; j++ ) {
4595 if ( spsig[j] == 0.0 || (sphot[j]-spsky[j]) == 0.0 ) {
4596 spec[j] = 0.0 ;
4597 flagchan[nflag++] = j ;
4598 }
4599 else {
4600 spec[j] = ( ( spref[j] - spsig[j] ) / spsig[j] )
4601 * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
4602 }
4603 }
4604 ref->specCol_.put( *p, spec ) ;
4605 ref->tsysCol_.put( *p, tsys ) ;
4606 if ( nflag > 0 ) {
4607 Vector<uChar> fl = ref->flagsCol_( *p ) ;
4608 for ( unsigned int j = 0 ; j < nflag ; j++ ) {
4609 fl[flagchan[j]] = (uChar)True ;
4610 }
4611 ref->flagsCol_.put( *p, fl ) ;
4612 }
4613 nflag = 0 ;
4614 p++ ;
4615 }
4616}
4617
4618void STMath::copyRows( Table &out,
4619 const Table &in,
4620 uInt startout,
4621 uInt startin,
4622 uInt nrow,
4623 Bool copySpectra,
4624 Bool copyFlagtra,
4625 Bool copyTsys )
4626{
4627 uInt nexclude = 0 ;
4628 Block<String> excludeColsBlock( 3 ) ;
4629 if ( !copySpectra ) {
4630 excludeColsBlock[nexclude] = "SPECTRA" ;
4631 nexclude++ ;
4632 }
4633 if ( !copyFlagtra ) {
4634 excludeColsBlock[nexclude] = "FLAGTRA" ;
4635 nexclude++ ;
4636 }
4637 if ( !copyTsys ) {
4638 excludeColsBlock[nexclude] = "TSYS" ;
4639 nexclude++ ;
4640 }
4641 // if ( nexclude < 3 ) {
4642 // excludeCols.resize( nexclude, True ) ;
4643 // }
4644 Vector<String> excludeCols( IPosition(1,nexclude),
4645 excludeColsBlock.storage(),
4646 SHARE ) ;
4647// cout << "excludeCols=" << excludeCols << endl ;
4648 TableRow rowout( out, excludeCols, True ) ;
4649 ROTableRow rowin( in, excludeCols, True ) ;
4650 uInt rin = startin ;
4651 uInt rout = startout ;
4652 for ( uInt i = 0 ; i < nrow ; i++ ) {
4653 rowin.get( rin ) ;
4654 rowout.putMatchingFields( rout, rowin.record() ) ;
4655 rin++ ;
4656 rout++ ;
4657 }
4658}
4659
4660CountedPtr<Scantable> STMath::averageWithinSession( CountedPtr<Scantable> &s,
4661 vector<bool> &mask,
4662 string weight )
4663{
4664 // prepare output table
4665 bool insitu = insitu_ ;
4666 insitu_ = false ;
4667 CountedPtr<Scantable> a = getScantable( s, true ) ;
4668 insitu_ = insitu ;
4669 Table &atab = a->table() ;
4670 ScalarColumn<Double> timeColOut( atab, "TIME" ) ;
4671
4672 if ( s->nrow() == 0 )
4673 return a ;
4674
4675 // setup RowAccumulator
4676 WeightType wtype = stringToWeight( weight ) ;
4677 RowAccumulator acc( wtype ) ;
4678 Vector<Bool> cmask( mask ) ;
4679 acc.setUserMask( cmask ) ;
4680
4681 vector<string> cols( 3 ) ;
4682 cols[0] = "IFNO" ;
4683 cols[1] = "POLNO" ;
4684 cols[2] = "BEAMNO" ;
4685 STIdxIterAcc iter( s, cols ) ;
4686
4687 Table ttab = s->table() ;
4688 ROScalarColumn<Double> *timeCol = new ROScalarColumn<Double>( ttab, "TIME" ) ;
4689 Vector<Double> timeVec = timeCol->getColumn() ;
4690 delete timeCol ;
4691 Vector<Double> interval = s->integrCol_.getColumn() ;
4692 uInt nrow = timeVec.nelements() ;
4693 uInt outrow = 0 ;
4694
4695 while( !iter.pastEnd() ) {
4696
4697 Vector<uInt> rows = iter.getRows( SHARE ) ;
4698 uInt len = rows.nelements() ;
4699
4700 if ( len == 0 ) {
4701 iter.next() ;
4702 continue ;
4703 }
4704
4705 uInt nchan = s->nchan(s->getIF(rows[0])) ;
4706 Vector<uChar> flag( nchan ) ;
4707 Vector<Bool> bflag( nchan ) ;
4708 Vector<Float> spec( nchan ) ;
4709 Vector<Float> tsys( nchan ) ;
4710
4711 Vector<Double> timeSep( len-1 ) ;
4712 for ( uInt i = 0 ; i < len-1 ; i++ ) {
4713 timeSep[i] = timeVec[rows[i+1]] - timeVec[rows[i]] ;
4714 }
4715
4716 uInt irow ;
4717 uInt jrow ;
4718 for ( uInt i = 0 ; i < len-1 ; i++ ) {
4719 irow = rows[i] ;
4720 jrow = rows[i+1] ;
4721 // accumulate data
4722 s->flagsCol_.get( irow, flag ) ;
4723 convertArray( bflag, flag ) ;
4724 s->specCol_.get( irow, spec ) ;
4725 tsys.assign( s->tsysCol_( irow ) ) ;
4726 if ( !allEQ(bflag,True) )
4727 acc.add( spec, !bflag, tsys, interval[irow], timeVec[irow] ) ;
4728 double gap = 2.0 * 86400.0 * timeSep[i] / ( interval[jrow] + interval[irow] ) ;
4729 //cout << "gap[" << i << "]=" << setw(5) << gap << endl ;
4730 if ( gap > 1.1 ) {
4731 //cout << "detected gap between " << i << " and " << i+1 << endl ;
4732 // put data to output table
4733 // reset RowAccumulator
4734 if ( acc.state() ) {
4735 atab.addRow() ;
4736 copyRows( atab, ttab, outrow, irow, 1, False, False, False ) ;
4737 acc.replaceNaN() ;
4738 const Vector<Bool> &msk = acc.getMask() ;
4739 convertArray( flag, !msk ) ;
4740 for (uInt k = 0; k < nchan; ++k) {
4741 uChar userFlag = 1 << 7;
4742 if (msk[k]==True) userFlag = 0 << 7;
4743 flag(k) = userFlag;
4744 }
4745 a->flagsCol_.put( outrow, flag ) ;
4746 a->specCol_.put( outrow, acc.getSpectrum() ) ;
4747 a->tsysCol_.put( outrow, acc.getTsys() ) ;
4748 a->integrCol_.put( outrow, acc.getInterval() ) ;
4749 timeColOut.put( outrow, acc.getTime() ) ;
4750 a->cycleCol_.put( outrow, 0 ) ;
4751 }
4752 acc.reset() ;
4753 outrow++ ;
4754 }
4755 }
4756
4757 // accumulate and add last data
4758 irow = rows[len-1] ;
4759 s->flagsCol_.get( irow, flag ) ;
4760 convertArray( bflag, flag ) ;
4761 s->specCol_.get( irow, spec ) ;
4762 tsys.assign( s->tsysCol_( irow ) ) ;
4763 if (!allEQ(bflag,True) )
4764 acc.add( spec, !bflag, tsys, interval[irow], timeVec[irow] ) ;
4765 if ( acc.state() ) {
4766 atab.addRow() ;
4767 copyRows( atab, ttab, outrow, irow, 1, False, False, False ) ;
4768 acc.replaceNaN() ;
4769 const Vector<Bool> &msk = acc.getMask() ;
4770 convertArray( flag, !msk ) ;
4771 for (uInt k = 0; k < nchan; ++k) {
4772 uChar userFlag = 1 << 7;
4773 if (msk[k]==True) userFlag = 0 << 7;
4774 flag(k) = userFlag;
4775 }
4776 a->flagsCol_.put( outrow, flag ) ;
4777 a->specCol_.put( outrow, acc.getSpectrum() ) ;
4778 a->tsysCol_.put( outrow, acc.getTsys() ) ;
4779 a->integrCol_.put( outrow, acc.getInterval() ) ;
4780 timeColOut.put( outrow, acc.getTime() ) ;
4781 a->cycleCol_.put( outrow, 0 ) ;
4782 }
4783 acc.reset() ;
4784 outrow++ ;
4785
4786 iter.next() ;
4787 }
4788
4789 return a ;
4790}
Note: See TracBrowser for help on using the repository browser.