source: branches/hpc34/src/STMath.cpp@ 2667

Last change on this file since 2667 was 2601, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs: sdsmooth unit test

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Minor speedup of STMath::smooth.
Create Vector outside the loop.


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