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

Last change on this file since 2552 was 2551, checked in by Takeshi Nakazato, 13 years ago

New Development: No

JIRA Issue: Yes CAS-4195

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Load whole (averaged) off spectra on memory instead to get one spectrum
when it is needed since those are not so heavy. It might reduce number of
call of ArrayColumn::get().
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 192.1 KB
Line 
1//
2// C++ Implementation: STMath
3//
4// Description:
5//
6//
7// Author: Malte Marquarding <asap@atnf.csiro.au>, (C) 2006
8//
9// Copyright: See COPYING file that comes with this distribution
10//
11//
12
13#include <sstream>
14
15#include <casa/iomanip.h>
16#include <casa/Arrays/MaskArrLogi.h>
17#include <casa/Arrays/MaskArrMath.h>
18#include <casa/Arrays/ArrayLogical.h>
19#include <casa/Arrays/ArrayMath.h>
20#include <casa/Arrays/Slice.h>
21#include <casa/Arrays/Slicer.h>
22#include <casa/BasicSL/String.h>
23#include <casa/Containers/Block.h>
24#include <casa/Containers/RecordField.h>
25#include <casa/Exceptions/Error.h>
26#include <casa/Logging/LogIO.h>
27
28#include <coordinates/Coordinates/CoordinateSystem.h>
29#include <coordinates/Coordinates/CoordinateUtil.h>
30#include <coordinates/Coordinates/FrequencyAligner.h>
31#include <coordinates/Coordinates/SpectralCoordinate.h>
32
33#include <lattices/Lattices/LatticeUtilities.h>
34
35#include <scimath/Functionals/Polynomial.h>
36#include <scimath/Mathematics/Convolver.h>
37#include <scimath/Mathematics/VectorKernel.h>
38
39#include <tables/Tables/ExprNode.h>
40#include <tables/Tables/ReadAsciiTable.h>
41#include <tables/Tables/TableCopy.h>
42#include <tables/Tables/TableIter.h>
43#include <tables/Tables/TableParse.h>
44#include <tables/Tables/TableRecord.h>
45#include <tables/Tables/TableRow.h>
46#include <tables/Tables/TableVector.h>
47#include <tables/Tables/TabVecMath.h>
48
49#include <atnf/PKSIO/SrcType.h>
50
51#include "RowAccumulator.h"
52#include "STAttr.h"
53#include "STMath.h"
54#include "STSelector.h"
55#include "Accelerator.h"
56#include "STIdxIter.h"
57
58using namespace casa;
59using namespace asap;
60
61
62// 2012/02/17 TN
63// Since STGrid is implemented, average doesn't consider direction
64// when accumulating
65// tolerance for direction comparison (rad)
66// #define TOL_OTF 1.0e-15
67// #define TOL_POINT 2.9088821e-4 // 1 arcmin
68
69STMath::STMath(bool insitu) :
70 insitu_(insitu)
71{
72}
73
74
75STMath::~STMath()
76{
77}
78
79CountedPtr<Scantable>
80STMath::average( const std::vector<CountedPtr<Scantable> >& in,
81 const std::vector<bool>& mask,
82 const std::string& weight,
83 const std::string& avmode)
84{
85// double t0, t1 ;
86// t0 = mathutil::gettimeofday_sec() ;
87
88 LogIO os( LogOrigin( "STMath", "average()", WHERE ) ) ;
89 if ( avmode == "SCAN" && in.size() != 1 )
90 throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
91 "Use merge first."));
92 WeightType wtype = stringToWeight(weight);
93
94 // 2012/02/17 TN
95 // Since STGrid is implemented, average doesn't consider direction
96 // when accumulating
97 // check if OTF observation
98// String obstype = in[0]->getHeader().obstype ;
99// Double tol = 0.0 ;
100// if ( (obstype.find( "OTF" ) != String::npos) || (obstype.find( "OBSERVE_TARGET" ) != String::npos) ) {
101// tol = TOL_OTF ;
102// }
103// else {
104// tol = TOL_POINT ;
105// }
106
107 // output
108 // clone as this is non insitu
109 bool insitu = insitu_;
110 setInsitu(false);
111 CountedPtr< Scantable > out = getScantable(in[0], true);
112 setInsitu(insitu);
113 std::vector<CountedPtr<Scantable> >::const_iterator stit = in.begin();
114 ++stit;
115 while ( stit != in.end() ) {
116 out->appendToHistoryTable((*stit)->history());
117 ++stit;
118 }
119
120 Table& tout = out->table();
121
122 /// @todo check if all scantables are conformant
123
124 ArrayColumn<Float> specColOut(tout,"SPECTRA");
125 ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
126 ArrayColumn<Float> tsysColOut(tout,"TSYS");
127 ScalarColumn<Double> mjdColOut(tout,"TIME");
128 ScalarColumn<Double> intColOut(tout,"INTERVAL");
129 ScalarColumn<uInt> cycColOut(tout,"CYCLENO");
130 ScalarColumn<uInt> scanColOut(tout,"SCANNO");
131
132 // set up the output table rows. These are based on the structure of the
133 // FIRST scantable in the vector
134 const Table& baset = in[0]->table();
135
136 RowAccumulator acc(wtype);
137 Vector<Bool> cmask(mask);
138 acc.setUserMask(cmask);
139// ROTableRow row(tout);
140 ROArrayColumn<Float> specCol, tsysCol;
141 ROArrayColumn<uChar> flagCol;
142 ROScalarColumn<Double> mjdCol, intCol;
143 ROScalarColumn<Int> scanIDCol;
144
145 //Vector<uInt> rowstodelete;
146 Block<uInt> rowstodelB( in[0]->nrow() ) ;
147 uInt nrowdel = 0 ;
148
149// Block<String> cols(3);
150 vector<string> cols(3) ;
151 cols[0] = String("BEAMNO");
152 cols[1] = String("IFNO");
153 cols[2] = String("POLNO");
154 if ( avmode == "SOURCE" ) {
155 cols.resize(4);
156 cols[3] = String("SRCNAME");
157 }
158 if ( avmode == "SCAN" && in.size() == 1) {
159 //cols.resize(4);
160 //cols[3] = String("SCANNO");
161 cols.resize(5);
162 cols[3] = String("SRCNAME");
163 cols[4] = String("SCANNO");
164 }
165 uInt outrowCount = 0;
166 // use STIdxIterExAcc instead of TableIterator
167 STIdxIterExAcc iter( in[0], cols ) ;
168// double t2 = 0 ;
169// double t3 = 0 ;
170// double t4 = 0 ;
171// double t5 = 0 ;
172// TableIterator iter(baset, cols);
173// int count = 0 ;
174 while (!iter.pastEnd()) {
175 Vector<uInt> rows = iter.getRows( SHARE ) ;
176 if ( rows.nelements() == 0 ) {
177 iter.next() ;
178 continue ;
179 }
180 Vector<uInt> current = iter.current() ;
181 String srcname = iter.getSrcName() ;
182 //Table subt = iter.table();
183 // copy the first row of this selection into the new table
184 tout.addRow();
185// t4 = mathutil::gettimeofday_sec() ;
186 // skip to copy SPECTRA, FLAGTRA, and TSYS since those heavy columns are
187 // overwritten in the following process
188 copyRows( tout, baset, outrowCount, rows[0], 1, False, False, False ) ;
189// t5 += mathutil::gettimeofday_sec() - t4 ;
190 // re-index to 0
191 if ( avmode != "SCAN" && avmode != "SOURCE" ) {
192 scanColOut.put(outrowCount, uInt(0));
193 }
194
195 // 2012/02/17 TN
196 // Since STGrid is implemented, average doesn't consider direction
197 // when accumulating
198// MDirection::ScalarColumn dircol ;
199// dircol.attach( subt, "DIRECTION" ) ;
200// Int length = subt.nrow() ;
201// vector< Vector<Double> > dirs ;
202// vector<int> indexes ;
203// for ( Int i = 0 ; i < length ; i++ ) {
204// Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
205// //os << << count++ << ": " ;
206// //os << "[" << t[0] << "," << t[1] << "]" << LogIO::POST ;
207// bool adddir = true ;
208// for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
209// //if ( allTrue( t == dirs[j] ) ) {
210// Double dx = t[0] - dirs[j][0] ;
211// Double dy = t[1] - dirs[j][1] ;
212// Double dd = sqrt( dx * dx + dy * dy ) ;
213// //if ( allNearAbs( t, dirs[j], tol ) ) {
214// if ( dd <= tol ) {
215// adddir = false ;
216// break ;
217// }
218// }
219// if ( adddir ) {
220// dirs.push_back( t ) ;
221// indexes.push_back( i ) ;
222// }
223// }
224// uInt rowNum = dirs.size() ;
225// tout.addRow( rowNum ) ;
226// for ( uInt i = 0 ; i < rowNum ; i++ ) {
227// TableCopy::copyRows( tout, subt, outrowCount+i, indexes[i], 1 ) ;
228// // re-index to 0
229// if ( avmode != "SCAN" && avmode != "SOURCE" ) {
230// scanColOut.put(outrowCount+i, uInt(0));
231// }
232// }
233// outrowCount += rowNum ;
234
235 // merge loop
236 uInt i = outrowCount ;
237 // in[0] is already selected by iterator
238 specCol.attach(baset,"SPECTRA");
239 flagCol.attach(baset,"FLAGTRA");
240 tsysCol.attach(baset,"TSYS");
241 intCol.attach(baset,"INTERVAL");
242 mjdCol.attach(baset,"TIME");
243 Vector<Float> spec,tsys;
244 Vector<uChar> flag;
245 Double inter,time;
246
247 for (uInt l = 0; l < rows.nelements(); ++l ) {
248 uInt k = rows[l] ;
249 flagCol.get(k, flag);
250 Vector<Bool> bflag(flag.shape());
251 convertArray(bflag, flag);
252 /*
253 if ( allEQ(bflag, True) ) {
254 continue;//don't accumulate
255 }
256 */
257 specCol.get(k, spec);
258 tsysCol.get(k, tsys);
259 intCol.get(k, inter);
260 mjdCol.get(k, time);
261 // spectrum has to be added last to enable weighting by the other values
262// t2 = mathutil::gettimeofday_sec() ;
263 acc.add(spec, !bflag, tsys, inter, time);
264// t3 += mathutil::gettimeofday_sec() - t2 ;
265
266 }
267
268
269 // in[0] is already selected by TableIterator so that index is
270 // started from 1
271 for ( int j=1; j < int(in.size()); ++j ) {
272 const Table& tin = in[j]->table();
273 //const TableRecord& rec = row.get(i);
274 ROScalarColumn<Double> tmp(tin, "TIME");
275 Double td;tmp.get(0,td);
276
277#if 1
278 static char const*const colNames1[] = { "IFNO", "BEAMNO", "POLNO" };
279 //uInt const values1[] = { rec.asuInt("IFNO"), rec.asuInt("BEAMNO"), rec.asuInt("POLNO") };
280 uInt const values1[] = { current[1], current[0], current[2] };
281 SingleTypeEqPredicate<uInt, 3> myPred(tin, colNames1, values1);
282 CustomTableExprNodeRep myNodeRep(tin, myPred);
283 myNodeRep.link(); // to avoid automatic delete when myExpr is destructed.
284 CustomTableExprNode myExpr(myNodeRep);
285 Table basesubt = tin(myExpr);
286#else
287// Table basesubt = tin( tin.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
288// && tin.col("IFNO") == Int(rec.asuInt("IFNO"))
289// && tin.col("POLNO") == Int(rec.asuInt("POLNO")) );
290 Table basesubt = tin( tin.col("BEAMNO") == current[0]
291 && tin.col("IFNO") == current[1]
292 && tin.col("POLNO") == current[2] );
293#endif
294 Table subt;
295 if ( avmode == "SOURCE") {
296// subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME"));
297 subt = basesubt( basesubt.col("SRCNAME") == srcname );
298
299 } else if (avmode == "SCAN") {
300// subt = basesubt( basesubt.col("SRCNAME") == rec.asString("SRCNAME")
301// && basesubt.col("SCANNO") == Int(rec.asuInt("SCANNO")) );
302 subt = basesubt( basesubt.col("SRCNAME") == srcname
303 && basesubt.col("SCANNO") == current[4] );
304 } else {
305 subt = basesubt;
306 }
307
308 // 2012/02/17 TN
309 // Since STGrid is implemented, average doesn't consider direction
310 // when accumulating
311// vector<uInt> removeRows ;
312// uInt nrsubt = subt.nrow() ;
313// for ( uInt irow = 0 ; irow < nrsubt ; irow++ ) {
314// //if ( !allTrue((subt.col("DIRECTION").getArrayDouble(TableExprId(irow)))==rec.asArrayDouble("DIRECTION")) ) {
315// Vector<Double> x0 = (subt.col("DIRECTION").getArrayDouble(TableExprId(irow))) ;
316// Vector<Double> x1 = rec.asArrayDouble("DIRECTION") ;
317// double dx = x0[0] - x1[0];
318// double dy = x0[1] - x1[1];
319// Double dd = sqrt( dx * dx + dy * dy ) ;
320// //if ( !allNearAbs((subt.col("DIRECTION").getArrayDouble(TableExprId(irow))), rec.asArrayDouble("DIRECTION"), tol ) ) {
321// if ( dd > tol ) {
322// removeRows.push_back( irow ) ;
323// }
324// }
325// if ( removeRows.size() != 0 ) {
326// subt.removeRow( removeRows ) ;
327// }
328
329// if ( nrsubt == removeRows.size() )
330// throw(AipsError("Averaging data is empty.")) ;
331
332 specCol.attach(subt,"SPECTRA");
333 flagCol.attach(subt,"FLAGTRA");
334 tsysCol.attach(subt,"TSYS");
335 intCol.attach(subt,"INTERVAL");
336 mjdCol.attach(subt,"TIME");
337 for (uInt k = 0; k < subt.nrow(); ++k ) {
338 flagCol.get(k, flag);
339 Vector<Bool> bflag(flag.shape());
340 convertArray(bflag, flag);
341 /*
342 if ( allEQ(bflag, True) ) {
343 continue;//don't accumulate
344 }
345 */
346 specCol.get(k, spec);
347 tsysCol.get(k, tsys);
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 for (uInt i=0; i < tout.nrow(); ++i ) {
2023 MaskedArray<Float> main = maskedArray(specCol(i), flagCol(i));
2024 MaskedArray<Float> maout;
2025 LatticeUtilities::bin(maout, main, 0, Int(width));
2026 /// @todo implement channel based tsys binning
2027 specCol.put(i, maout.getArray());
2028 flagCol.put(i, flagsFromMA(maout));
2029 // take only the first binned spectrum's length for the deprecated
2030 // global header item nChan
2031 if (i==0) tout.rwKeywordSet().define(String("nChan"),
2032 Int(maout.getArray().nelements()));
2033 }
2034 return out;
2035}
2036
2037CountedPtr< Scantable > STMath::resample( const CountedPtr< Scantable >& in,
2038 const std::string& method,
2039 float width )
2040//
2041// Should add the possibility of width being specified in km/s. This means
2042// that for each freqID (SpectralCoordinate) we will need to convert to an
2043// average channel width (say at the reference pixel). Then we would need
2044// to be careful to make sure each spectrum (of different freqID)
2045// is the same length.
2046//
2047{
2048 //InterpolateArray1D<Double,Float>::InterpolationMethod interp;
2049 Int interpMethod(stringToIMethod(method));
2050
2051 CountedPtr< Scantable > out = getScantable(in, false);
2052 Table& tout = out->table();
2053
2054// Resample SpectralCoordinates (one per freqID)
2055 out->frequencies().rescale(width, "RESAMPLE");
2056 TableIterator iter(tout, "IFNO");
2057 TableRow row(tout);
2058 while ( !iter.pastEnd() ) {
2059 Table tab = iter.table();
2060 ArrayColumn<Float> specCol(tab, "SPECTRA");
2061 //ArrayColumn<Float> tsysCol(tout, "TSYS");
2062 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2063 Vector<Float> spec;
2064 Vector<uChar> flag;
2065 specCol.get(0,spec); // the number of channels should be constant per IF
2066 uInt nChanIn = spec.nelements();
2067 Vector<Float> xIn(nChanIn); indgen(xIn);
2068 Int fac = Int(nChanIn/width);
2069 Vector<Float> xOut(fac+10); // 10 to be safe - resize later
2070 uInt k = 0;
2071 Float x = 0.0;
2072 while (x < Float(nChanIn) ) {
2073 xOut(k) = x;
2074 k++;
2075 x += width;
2076 }
2077 uInt nChanOut = k;
2078 xOut.resize(nChanOut, True);
2079 // process all rows for this IFNO
2080 Vector<Float> specOut;
2081 Vector<Bool> maskOut;
2082 Vector<uChar> flagOut;
2083 for (uInt i=0; i < tab.nrow(); ++i) {
2084 specCol.get(i, spec);
2085 flagCol.get(i, flag);
2086 Vector<Bool> mask(flag.nelements());
2087 convertArray(mask, flag);
2088
2089 IPosition shapeIn(spec.shape());
2090 //sh.nchan = nChanOut;
2091 InterpolateArray1D<Float,Float>::interpolate(specOut, maskOut, xOut,
2092 xIn, spec, mask,
2093 interpMethod, True, True);
2094 /// @todo do the same for channel based Tsys
2095 flagOut.resize(maskOut.nelements());
2096 convertArray(flagOut, maskOut);
2097 specCol.put(i, specOut);
2098 flagCol.put(i, flagOut);
2099 }
2100 ++iter;
2101 }
2102
2103 return out;
2104}
2105
2106STMath::imethod STMath::stringToIMethod(const std::string& in)
2107{
2108 static STMath::imap lookup;
2109
2110 // initialize the lookup table if necessary
2111 if ( lookup.empty() ) {
2112 lookup["nearest"] = InterpolateArray1D<Double,Float>::nearestNeighbour;
2113 lookup["linear"] = InterpolateArray1D<Double,Float>::linear;
2114 lookup["cubic"] = InterpolateArray1D<Double,Float>::cubic;
2115 lookup["spline"] = InterpolateArray1D<Double,Float>::spline;
2116 }
2117
2118 STMath::imap::const_iterator iter = lookup.find(in);
2119
2120 if ( lookup.end() == iter ) {
2121 std::string message = in;
2122 message += " is not a valid interpolation mode";
2123 throw(AipsError(message));
2124 }
2125 return iter->second;
2126}
2127
2128WeightType STMath::stringToWeight(const std::string& in)
2129{
2130 static std::map<std::string, WeightType> lookup;
2131
2132 // initialize the lookup table if necessary
2133 if ( lookup.empty() ) {
2134 lookup["NONE"] = asap::W_NONE;
2135 lookup["TINT"] = asap::W_TINT;
2136 lookup["TINTSYS"] = asap::W_TINTSYS;
2137 lookup["TSYS"] = asap::W_TSYS;
2138 lookup["VAR"] = asap::W_VAR;
2139 }
2140
2141 std::map<std::string, WeightType>::const_iterator iter = lookup.find(in);
2142
2143 if ( lookup.end() == iter ) {
2144 std::string message = in;
2145 message += " is not a valid weighting mode";
2146 throw(AipsError(message));
2147 }
2148 return iter->second;
2149}
2150
2151CountedPtr< Scantable > STMath::gainElevation( const CountedPtr< Scantable >& in,
2152 const vector< float > & coeff,
2153 const std::string & filename,
2154 const std::string& method)
2155{
2156 // Get elevation data from Scantable and convert to degrees
2157 CountedPtr< Scantable > out = getScantable(in, false);
2158 Table& tab = out->table();
2159 ROScalarColumn<Float> elev(tab, "ELEVATION");
2160 Vector<Float> x = elev.getColumn();
2161 x *= Float(180 / C::pi); // Degrees
2162
2163 Vector<Float> coeffs(coeff);
2164 const uInt nc = coeffs.nelements();
2165 if ( filename.length() > 0 && nc > 0 ) {
2166 throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
2167 }
2168
2169 // Correct
2170 if ( nc > 0 || filename.length() == 0 ) {
2171 // Find instrument
2172 Bool throwit = True;
2173 Instrument inst =
2174 STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
2175 throwit);
2176
2177 // Set polynomial
2178 Polynomial<Float>* ppoly = 0;
2179 Vector<Float> coeff;
2180 String msg;
2181 if ( nc > 0 ) {
2182 ppoly = new Polynomial<Float>(nc-1);
2183 coeff = coeffs;
2184 msg = String("user");
2185 } else {
2186 STAttr sdAttr;
2187 coeff = sdAttr.gainElevationPoly(inst);
2188 ppoly = new Polynomial<Float>(coeff.nelements()-1);
2189 msg = String("built in");
2190 }
2191
2192 if ( coeff.nelements() > 0 ) {
2193 ppoly->setCoefficients(coeff);
2194 } else {
2195 delete ppoly;
2196 throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
2197 }
2198 ostringstream oss;
2199 oss << "Making polynomial correction with " << msg << " coefficients:" << endl;
2200 oss << " " << coeff;
2201 pushLog(String(oss));
2202 const uInt nrow = tab.nrow();
2203 Vector<Float> factor(nrow);
2204 for ( uInt i=0; i < nrow; ++i ) {
2205 factor[i] = 1.0 / (*ppoly)(x[i]);
2206 }
2207 delete ppoly;
2208 scaleByVector(tab, factor, true);
2209
2210 } else {
2211 // Read and correct
2212 pushLog("Making correction from ascii Table");
2213 scaleFromAsciiTable(tab, filename, method, x, true);
2214 }
2215 return out;
2216}
2217
2218void STMath::scaleFromAsciiTable(Table& in, const std::string& filename,
2219 const std::string& method,
2220 const Vector<Float>& xout, bool dotsys)
2221{
2222
2223// Read gain-elevation ascii file data into a Table.
2224
2225 String formatString;
2226 Table tbl = readAsciiTable(formatString, Table::Memory, filename, "", "", False);
2227 scaleFromTable(in, tbl, method, xout, dotsys);
2228}
2229
2230void STMath::scaleFromTable(Table& in,
2231 const Table& table,
2232 const std::string& method,
2233 const Vector<Float>& xout, bool dotsys)
2234{
2235
2236 ROScalarColumn<Float> geElCol(table, "ELEVATION");
2237 ROScalarColumn<Float> geFacCol(table, "FACTOR");
2238 Vector<Float> xin = geElCol.getColumn();
2239 Vector<Float> yin = geFacCol.getColumn();
2240 Vector<Bool> maskin(xin.nelements(),True);
2241
2242 // Interpolate (and extrapolate) with desired method
2243
2244 InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
2245
2246 Vector<Float> yout;
2247 Vector<Bool> maskout;
2248 InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout,
2249 xin, yin, maskin, interp,
2250 True, True);
2251
2252 scaleByVector(in, Float(1.0)/yout, dotsys);
2253}
2254
2255void STMath::scaleByVector( Table& in,
2256 const Vector< Float >& factor,
2257 bool dotsys )
2258{
2259 uInt nrow = in.nrow();
2260 if ( factor.nelements() != nrow ) {
2261 throw(AipsError("factors.nelements() != table.nelements()"));
2262 }
2263 ArrayColumn<Float> specCol(in, "SPECTRA");
2264 ArrayColumn<uChar> flagCol(in, "FLAGTRA");
2265 ArrayColumn<Float> tsysCol(in, "TSYS");
2266 for (uInt i=0; i < nrow; ++i) {
2267 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
2268 ma *= factor[i];
2269 specCol.put(i, ma.getArray());
2270 flagCol.put(i, flagsFromMA(ma));
2271 if ( dotsys ) {
2272 Vector<Float> tsys = tsysCol(i);
2273 tsys *= factor[i];
2274 tsysCol.put(i,tsys);
2275 }
2276 }
2277}
2278
2279CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in,
2280 float d, float etaap,
2281 float jyperk )
2282{
2283 CountedPtr< Scantable > out = getScantable(in, false);
2284 Table& tab = in->table();
2285 Unit fluxUnit(tab.keywordSet().asString("FluxUnit"));
2286 Unit K(String("K"));
2287 Unit JY(String("Jy"));
2288
2289 bool tokelvin = true;
2290 Double cfac = 1.0;
2291
2292 if ( fluxUnit == JY ) {
2293 pushLog("Converting to K");
2294 Quantum<Double> t(1.0,fluxUnit);
2295 Quantum<Double> t2 = t.get(JY);
2296 cfac = (t2 / t).getValue(); // value to Jy
2297
2298 tokelvin = true;
2299 out->setFluxUnit("K");
2300 } else if ( fluxUnit == K ) {
2301 pushLog("Converting to Jy");
2302 Quantum<Double> t(1.0,fluxUnit);
2303 Quantum<Double> t2 = t.get(K);
2304 cfac = (t2 / t).getValue(); // value to K
2305
2306 tokelvin = false;
2307 out->setFluxUnit("Jy");
2308 } else {
2309 throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
2310 }
2311 // Make sure input values are converted to either Jy or K first...
2312 Float factor = cfac;
2313
2314 // Select method
2315 if (jyperk > 0.0) {
2316 factor *= jyperk;
2317 if ( tokelvin ) factor = 1.0 / jyperk;
2318 ostringstream oss;
2319 oss << "Jy/K = " << jyperk;
2320 pushLog(String(oss));
2321 Vector<Float> factors(tab.nrow(), factor);
2322 scaleByVector(tab,factors, false);
2323 } else if ( etaap > 0.0) {
2324 if (d < 0) {
2325 Instrument inst =
2326 STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
2327 True);
2328 STAttr sda;
2329 d = sda.diameter(inst);
2330 }
2331 jyperk = STAttr::findJyPerK(etaap, d);
2332 ostringstream oss;
2333 oss << "Jy/K = " << jyperk;
2334 pushLog(String(oss));
2335 factor *= jyperk;
2336 if ( tokelvin ) {
2337 factor = 1.0 / factor;
2338 }
2339 Vector<Float> factors(tab.nrow(), factor);
2340 scaleByVector(tab, factors, False);
2341 } else {
2342
2343 // OK now we must deal with automatic look up of values.
2344 // We must also deal with the fact that the factors need
2345 // to be computed per IF and may be different and may
2346 // change per integration.
2347
2348 pushLog("Looking up conversion factors");
2349 convertBrightnessUnits(out, tokelvin, cfac);
2350 }
2351
2352 return out;
2353}
2354
2355void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in,
2356 bool tokelvin, float cfac )
2357{
2358 Table& table = in->table();
2359 Instrument inst =
2360 STAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True);
2361 TableIterator iter(table, "FREQ_ID");
2362 STFrequencies stfreqs = in->frequencies();
2363 STAttr sdAtt;
2364 while (!iter.pastEnd()) {
2365 Table tab = iter.table();
2366 ArrayColumn<Float> specCol(tab, "SPECTRA");
2367 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2368 ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID");
2369 MEpoch::ROScalarColumn timeCol(tab, "TIME");
2370
2371 uInt freqid; freqidCol.get(0, freqid);
2372 Vector<Float> tmpspec; specCol.get(0, tmpspec);
2373 // STAttr.JyPerK has a Vector interface... change sometime.
2374 Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements()));
2375 for ( uInt i=0; i<tab.nrow(); ++i) {
2376 Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0];
2377 Float factor = cfac * jyperk;
2378 if ( tokelvin ) factor = Float(1.0) / factor;
2379 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
2380 ma *= factor;
2381 specCol.put(i, ma.getArray());
2382 flagCol.put(i, flagsFromMA(ma));
2383 }
2384 ++iter;
2385 }
2386}
2387
2388CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
2389 const std::vector<float>& tau )
2390{
2391 CountedPtr< Scantable > out = getScantable(in, false);
2392
2393 Table outtab = out->table();
2394
2395 const Int ntau = uInt(tau.size());
2396 std::vector<float>::const_iterator tauit = tau.begin();
2397 AlwaysAssert((ntau == 1 || ntau == in->nif() || ntau == in->nif() * in->npol()),
2398 AipsError);
2399 TableIterator iiter(outtab, "IFNO");
2400 while ( !iiter.pastEnd() ) {
2401 Table itab = iiter.table();
2402 TableIterator piter(itab, "POLNO");
2403 while ( !piter.pastEnd() ) {
2404 Table tab = piter.table();
2405 ROScalarColumn<Float> elev(tab, "ELEVATION");
2406 ArrayColumn<Float> specCol(tab, "SPECTRA");
2407 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2408 ArrayColumn<Float> tsysCol(tab, "TSYS");
2409 for ( uInt i=0; i<tab.nrow(); ++i) {
2410 Float zdist = Float(C::pi_2) - elev(i);
2411 Float factor = exp(*tauit/cos(zdist));
2412 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
2413 ma *= factor;
2414 specCol.put(i, ma.getArray());
2415 flagCol.put(i, flagsFromMA(ma));
2416 Vector<Float> tsys;
2417 tsysCol.get(i, tsys);
2418 tsys *= factor;
2419 tsysCol.put(i, tsys);
2420 }
2421 if (ntau == in->nif()*in->npol() ) {
2422 tauit++;
2423 }
2424 piter++;
2425 }
2426 if (ntau >= in->nif() ) {
2427 tauit++;
2428 }
2429 iiter++;
2430 }
2431 return out;
2432}
2433
2434CountedPtr< Scantable > STMath::smoothOther( const CountedPtr< Scantable >& in,
2435 const std::string& kernel,
2436 float width, int order)
2437{
2438 CountedPtr< Scantable > out = getScantable(in, false);
2439 Table table = out->table();
2440
2441 TableIterator iter(table, "IFNO");
2442 while (!iter.pastEnd()) {
2443 Table tab = iter.table();
2444 ArrayColumn<Float> specCol(tab, "SPECTRA");
2445 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2446 Vector<Float> spec;
2447 Vector<uChar> flag;
2448 for (uInt i = 0; i < tab.nrow(); ++i) {
2449 specCol.get(i, spec);
2450 flagCol.get(i, flag);
2451 Vector<Bool> mask(flag.nelements());
2452 convertArray(mask, flag);
2453 Vector<Float> specout;
2454 Vector<Bool> maskout;
2455 if (kernel == "hanning") {
2456 mathutil::hanning(specout, maskout, spec, !mask);
2457 convertArray(flag, !maskout);
2458 } else if (kernel == "rmedian") {
2459 mathutil::runningMedian(specout, maskout, spec , mask, width);
2460 convertArray(flag, maskout);
2461 } else if (kernel == "poly") {
2462 mathutil::polyfit(specout, maskout, spec, !mask, width, order);
2463 convertArray(flag, !maskout);
2464 }
2465
2466 for (uInt j = 0; j < flag.nelements(); ++j) {
2467 uChar userFlag = 1 << 7;
2468 if (maskout[j]==True) userFlag = 0 << 7;
2469 flag(j) = userFlag;
2470 }
2471
2472 flagCol.put(i, flag);
2473 specCol.put(i, specout);
2474 }
2475 ++iter;
2476 }
2477 return out;
2478}
2479
2480CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
2481 const std::string& kernel, float width,
2482 int order)
2483{
2484 if (kernel == "rmedian" || kernel == "hanning" || kernel == "poly") {
2485 return smoothOther(in, kernel, width, order);
2486 }
2487 CountedPtr< Scantable > out = getScantable(in, false);
2488 Table& table = out->table();
2489 VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
2490 // same IFNO should have same no of channels
2491 // this saves overhead
2492 TableIterator iter(table, "IFNO");
2493 while (!iter.pastEnd()) {
2494 Table tab = iter.table();
2495 ArrayColumn<Float> specCol(tab, "SPECTRA");
2496 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2497 Vector<Float> tmpspec; specCol.get(0, tmpspec);
2498 uInt nchan = tmpspec.nelements();
2499 Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
2500 Convolver<Float> conv(kvec, IPosition(1,nchan));
2501 Vector<Float> spec;
2502 Vector<uChar> flag;
2503 for ( uInt i=0; i<tab.nrow(); ++i) {
2504 specCol.get(i, spec);
2505 flagCol.get(i, flag);
2506 Vector<Bool> mask(flag.nelements());
2507 convertArray(mask, flag);
2508 Vector<Float> specout;
2509 mathutil::replaceMaskByZero(specout, mask);
2510 conv.linearConv(specout, spec);
2511 specCol.put(i, specout);
2512 }
2513 ++iter;
2514 }
2515 return out;
2516}
2517
2518CountedPtr< Scantable >
2519 STMath::merge( const std::vector< CountedPtr < Scantable > >& in )
2520{
2521 if ( in.size() < 2 ) {
2522 throw(AipsError("Need at least two scantables to perform a merge."));
2523 }
2524 std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
2525 bool insitu = insitu_;
2526 setInsitu(false);
2527 CountedPtr< Scantable > out = getScantable(*it, false);
2528 setInsitu(insitu);
2529 Table& tout = out->table();
2530 ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
2531 ScalarColumn<uInt> scannocol(tout,"SCANNO"), focusidcol(tout,"FOCUS_ID");
2532 // Renumber SCANNO to be 0-based
2533 Vector<uInt> scannos = scannocol.getColumn();
2534 uInt offset = min(scannos);
2535 scannos -= offset;
2536 scannocol.putColumn(scannos);
2537 uInt newscanno = max(scannos)+1;
2538 ++it;
2539 while ( it != in.end() ){
2540 if ( ! (*it)->conformant(*out) ) {
2541 // non conformant.
2542 //pushLog(String("Warning: Can't merge scantables as header info differs."));
2543 LogIO os( LogOrigin( "STMath", "merge()", WHERE ) ) ;
2544 os << LogIO::SEVERE << "Can't merge scantables as header informations (any one of AntennaName, Equinox, and FluxUnit) differ." << LogIO::EXCEPTION ;
2545 }
2546 out->appendToHistoryTable((*it)->history());
2547 const Table& tab = (*it)->table();
2548
2549 Block<String> cols(3);
2550 cols[0] = String("FREQ_ID");
2551 cols[1] = String("MOLECULE_ID");
2552 cols[2] = String("FOCUS_ID");
2553
2554 TableIterator scanit(tab, "SCANNO");
2555 while (!scanit.pastEnd()) {
2556 ScalarColumn<uInt> thescannocol(scanit.table(),"SCANNO");
2557 Vector<uInt> thescannos(thescannocol.nrow(),newscanno);
2558 thescannocol.putColumn(thescannos);
2559 TableIterator subit(scanit.table(), cols);
2560 while ( !subit.pastEnd() ) {
2561 uInt nrow = tout.nrow();
2562 Table thetab = subit.table();
2563 ROTableRow row(thetab);
2564 Vector<uInt> thecolvals(thetab.nrow());
2565 ScalarColumn<uInt> thefreqidcol(thetab,"FREQ_ID");
2566 ScalarColumn<uInt> themolidcol(thetab, "MOLECULE_ID");
2567 ScalarColumn<uInt> thefocusidcol(thetab,"FOCUS_ID");
2568 // The selected subset of table should have
2569 // the equal FREQ_ID, MOLECULE_ID, and FOCUS_ID values.
2570 const TableRecord& rec = row.get(0);
2571 // Set the proper FREQ_ID
2572 Double rv,rp,inc;
2573 (*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
2574 uInt id;
2575 id = out->frequencies().addEntry(rp, rv, inc);
2576 thecolvals = id;
2577 thefreqidcol.putColumn(thecolvals);
2578 // Set the proper MOLECULE_ID
2579 Vector<String> name,fname;Vector<Double> rf;
2580 (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
2581 id = out->molecules().addEntry(rf, name, fname);
2582 thecolvals = id;
2583 themolidcol.putColumn(thecolvals);
2584 // Set the proper FOCUS_ID
2585 Float fpa,frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
2586 (*it)->focus().getEntry(fpa, fax, ftan, frot, fhand, fmount,fuser,
2587 fxy, fxyp, rec.asuInt("FOCUS_ID"));
2588 id = out->focus().addEntry(fpa, fax, ftan, frot, fhand, fmount,fuser,
2589 fxy, fxyp);
2590 thecolvals = id;
2591 thefocusidcol.putColumn(thecolvals);
2592
2593 tout.addRow(thetab.nrow());
2594 TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
2595
2596 ++subit;
2597 }
2598 ++newscanno;
2599 ++scanit;
2600 }
2601 ++it;
2602 }
2603 return out;
2604}
2605
2606CountedPtr< Scantable >
2607 STMath::invertPhase( const CountedPtr < Scantable >& in )
2608{
2609 return applyToPol(in, &STPol::invertPhase, Float(0.0));
2610}
2611
2612CountedPtr< Scantable >
2613 STMath::rotateXYPhase( const CountedPtr < Scantable >& in, float phase )
2614{
2615 return applyToPol(in, &STPol::rotatePhase, Float(phase));
2616}
2617
2618CountedPtr< Scantable >
2619 STMath::rotateLinPolPhase( const CountedPtr < Scantable >& in, float phase )
2620{
2621 return applyToPol(in, &STPol::rotateLinPolPhase, Float(phase));
2622}
2623
2624CountedPtr< Scantable > STMath::applyToPol( const CountedPtr<Scantable>& in,
2625 STPol::polOperation fptr,
2626 Float phase )
2627{
2628 CountedPtr< Scantable > out = getScantable(in, false);
2629 Table& tout = out->table();
2630 Block<String> cols(4);
2631 cols[0] = String("SCANNO");
2632 cols[1] = String("BEAMNO");
2633 cols[2] = String("IFNO");
2634 cols[3] = String("CYCLENO");
2635 TableIterator iter(tout, cols);
2636 CountedPtr<STPol> stpol = STPol::getPolClass(out->factories_,
2637 out->getPolType() );
2638 while (!iter.pastEnd()) {
2639 Table t = iter.table();
2640 ArrayColumn<Float> speccol(t, "SPECTRA");
2641 ScalarColumn<uInt> focidcol(t, "FOCUS_ID");
2642 Matrix<Float> pols(speccol.getColumn());
2643 try {
2644 stpol->setSpectra(pols);
2645 Float fang,fhand;
2646 fang = in->focusTable_.getTotalAngle(focidcol(0));
2647 fhand = in->focusTable_.getFeedHand(focidcol(0));
2648 stpol->setPhaseCorrections(fang, fhand);
2649 // use a member function pointer in STPol. This only works on
2650 // the STPol pointer itself, not the Counted Pointer so
2651 // derefernce it.
2652 (&(*(stpol))->*fptr)(phase);
2653 speccol.putColumn(stpol->getSpectra());
2654 } catch (AipsError& e) {
2655 //delete stpol;stpol=0;
2656 throw(e);
2657 }
2658 ++iter;
2659 }
2660 //delete stpol;stpol=0;
2661 return out;
2662}
2663
2664CountedPtr< Scantable >
2665 STMath::swapPolarisations( const CountedPtr< Scantable > & in )
2666{
2667 CountedPtr< Scantable > out = getScantable(in, false);
2668 Table& tout = out->table();
2669 Table t0 = tout(tout.col("POLNO") == 0);
2670 Table t1 = tout(tout.col("POLNO") == 1);
2671 if ( t0.nrow() != t1.nrow() )
2672 throw(AipsError("Inconsistent number of polarisations"));
2673 ArrayColumn<Float> speccol0(t0, "SPECTRA");
2674 ArrayColumn<uChar> flagcol0(t0, "FLAGTRA");
2675 ArrayColumn<Float> speccol1(t1, "SPECTRA");
2676 ArrayColumn<uChar> flagcol1(t1, "FLAGTRA");
2677 Matrix<Float> s0 = speccol0.getColumn();
2678 Matrix<uChar> f0 = flagcol0.getColumn();
2679 speccol0.putColumn(speccol1.getColumn());
2680 flagcol0.putColumn(flagcol1.getColumn());
2681 speccol1.putColumn(s0);
2682 flagcol1.putColumn(f0);
2683 return out;
2684}
2685
2686CountedPtr< Scantable >
2687 STMath::averagePolarisations( const CountedPtr< Scantable > & in,
2688 const std::vector<bool>& mask,
2689 const std::string& weight )
2690{
2691 if (in->npol() < 2 )
2692 throw(AipsError("averagePolarisations can only be applied to two or more"
2693 "polarisations"));
2694 bool insitu = insitu_;
2695 setInsitu(false);
2696 CountedPtr< Scantable > pols = getScantable(in, true);
2697 setInsitu(insitu);
2698 Table& tout = pols->table();
2699 std::string taql = "SELECT FROM $1 WHERE POLNO IN [0,1]";
2700 Table tab = tableCommand(taql, in->table());
2701 if (tab.nrow() == 0 )
2702 throw(AipsError("Could not find any rows with POLNO==0 and POLNO==1"));
2703 TableCopy::copyRows(tout, tab);
2704 TableVector<uInt> vec(tout, "POLNO");
2705 vec = 0;
2706 pols->table_.rwKeywordSet().define("nPol", Int(1));
2707 pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
2708 //pols->table_.rwKeywordSet().define("POLTYPE", in->getPolType());
2709 std::vector<CountedPtr<Scantable> > vpols;
2710 vpols.push_back(pols);
2711 CountedPtr< Scantable > out = average(vpols, mask, weight, "SCAN");
2712 return out;
2713}
2714
2715CountedPtr< Scantable >
2716 STMath::averageBeams( const CountedPtr< Scantable > & in,
2717 const std::vector<bool>& mask,
2718 const std::string& weight )
2719{
2720 bool insitu = insitu_;
2721 setInsitu(false);
2722 CountedPtr< Scantable > beams = getScantable(in, false);
2723 setInsitu(insitu);
2724 Table& tout = beams->table();
2725 // give all rows the same BEAMNO
2726 TableVector<uInt> vec(tout, "BEAMNO");
2727 vec = 0;
2728 beams->table_.rwKeywordSet().define("nBeam", Int(1));
2729 std::vector<CountedPtr<Scantable> > vbeams;
2730 vbeams.push_back(beams);
2731 CountedPtr< Scantable > out = average(vbeams, mask, weight, "SCAN");
2732 return out;
2733}
2734
2735
2736CountedPtr< Scantable >
2737 asap::STMath::frequencyAlign( const CountedPtr< Scantable > & in,
2738 const std::string & refTime,
2739 const std::string & method)
2740{
2741 // clone as this is not working insitu
2742 bool insitu = insitu_;
2743 setInsitu(false);
2744 CountedPtr< Scantable > out = getScantable(in, false);
2745 setInsitu(insitu);
2746 Table& tout = out->table();
2747 // Get reference Epoch to time of first row or given String
2748 Unit DAY(String("d"));
2749 MEpoch::Ref epochRef(in->getTimeReference());
2750 MEpoch refEpoch;
2751 if (refTime.length()>0) {
2752 Quantum<Double> qt;
2753 if (MVTime::read(qt,refTime)) {
2754 MVEpoch mv(qt);
2755 refEpoch = MEpoch(mv, epochRef);
2756 } else {
2757 throw(AipsError("Invalid format for Epoch string"));
2758 }
2759 } else {
2760 refEpoch = in->timeCol_(0);
2761 }
2762 MPosition refPos = in->getAntennaPosition();
2763
2764 InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
2765 /*
2766 // Comment from MV.
2767 // the following code has been commented out because different FREQ_IDs have to be aligned together even
2768 // if the frame doesn't change. So far, lack of this check didn't cause any problems.
2769 // test if user frame is different to base frame
2770 if ( in->frequencies().getFrameString(true)
2771 == in->frequencies().getFrameString(false) ) {
2772 throw(AipsError("Can't convert as no output frame has been set"
2773 " (use set_freqframe) or it is aligned already."));
2774 }
2775 */
2776 MFrequency::Types system = in->frequencies().getFrame();
2777 MVTime mvt(refEpoch.getValue());
2778 String epochout = mvt.string(MVTime::YMD) + String(" (") + refEpoch.getRefString() + String(")");
2779 ostringstream oss;
2780 oss << "Aligned at reference Epoch " << epochout
2781 << " in frame " << MFrequency::showType(system);
2782 pushLog(String(oss));
2783 // set up the iterator
2784 Block<String> cols(4);
2785 // select by constant direction
2786 cols[0] = String("SRCNAME");
2787 cols[1] = String("BEAMNO");
2788 // select by IF ( no of channels varies over this )
2789 cols[2] = String("IFNO");
2790 // select by restfrequency
2791 cols[3] = String("MOLECULE_ID");
2792 TableIterator iter(tout, cols);
2793 while ( !iter.pastEnd() ) {
2794 Table t = iter.table();
2795 MDirection::ROScalarColumn dirCol(t, "DIRECTION");
2796 TableIterator fiter(t, "FREQ_ID");
2797 // determine nchan from the first row. This should work as
2798 // we are iterating over BEAMNO and IFNO // we should have constant direction
2799
2800 ROArrayColumn<Float> sCol(t, "SPECTRA");
2801 const MDirection direction = dirCol(0);
2802 const uInt nchan = sCol(0).nelements();
2803
2804 // skip operations if there is nothing to align
2805 if (fiter.pastEnd()) {
2806 continue;
2807 }
2808
2809 Table ftab = fiter.table();
2810 // align all frequency ids with respect to the first encountered id
2811 ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
2812 // get the SpectralCoordinate for the freqid, which we are iterating over
2813 SpectralCoordinate sC = in->frequencies().getSpectralCoordinate(freqidCol(0));
2814 FrequencyAligner<Float> fa( sC, nchan, refEpoch,
2815 direction, refPos, system );
2816 // realign the SpectralCoordinate and put into the output Scantable
2817 Vector<String> units(1);
2818 units = String("Hz");
2819 Bool linear=True;
2820 SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
2821 sc2.setWorldAxisUnits(units);
2822 const uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
2823 sc2.referenceValue()[0],
2824 sc2.increment()[0]);
2825 while ( !fiter.pastEnd() ) {
2826 ftab = fiter.table();
2827 // spectral coordinate for the current FREQ_ID
2828 ScalarColumn<uInt> freqidCol2(ftab, "FREQ_ID");
2829 sC = in->frequencies().getSpectralCoordinate(freqidCol2(0));
2830 // create the "global" abcissa for alignment with same FREQ_ID
2831 Vector<Double> abc(nchan);
2832 for (uInt i=0; i<nchan; i++) {
2833 Double w;
2834 sC.toWorld(w,Double(i));
2835 abc[i] = w;
2836 }
2837 TableVector<uInt> tvec(ftab, "FREQ_ID");
2838 // assign new frequency id to all rows
2839 tvec = id;
2840 // cache abcissa for same time stamps, so iterate over those
2841 TableIterator timeiter(ftab, "TIME");
2842 while ( !timeiter.pastEnd() ) {
2843 Table tab = timeiter.table();
2844 ArrayColumn<Float> specCol(tab, "SPECTRA");
2845 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2846 MEpoch::ROScalarColumn timeCol(tab, "TIME");
2847 // use align abcissa cache after the first row
2848 // these rows should be just be POLNO
2849 bool first = true;
2850 for (int i=0; i<int(tab.nrow()); ++i) {
2851 // input values
2852 Vector<uChar> flag = flagCol(i);
2853 Vector<Bool> mask(flag.shape());
2854 Vector<Float> specOut, spec;
2855 spec = specCol(i);
2856 Vector<Bool> maskOut;Vector<uChar> flagOut;
2857 convertArray(mask, flag);
2858 // alignment
2859 Bool ok = fa.align(specOut, maskOut, abc, spec,
2860 mask, timeCol(i), !first,
2861 interp, False);
2862 (void) ok; // unused stop compiler nagging
2863 // back into scantable
2864 flagOut.resize(maskOut.nelements());
2865 convertArray(flagOut, maskOut);
2866 flagCol.put(i, flagOut);
2867 specCol.put(i, specOut);
2868 // start abcissa caching
2869 first = false;
2870 }
2871 // next timestamp
2872 ++timeiter;
2873 }
2874 // next FREQ_ID
2875 ++fiter;
2876 }
2877 // next aligner
2878 ++iter;
2879 }
2880 // set this afterwards to ensure we are doing insitu correctly.
2881 out->frequencies().setFrame(system, true);
2882 return out;
2883}
2884
2885CountedPtr<Scantable>
2886 asap::STMath::convertPolarisation( const CountedPtr<Scantable>& in,
2887 const std::string & newtype )
2888{
2889 if (in->npol() != 2 && in->npol() != 4)
2890 throw(AipsError("Can only convert two or four polarisations."));
2891 if ( in->getPolType() == newtype )
2892 throw(AipsError("No need to convert."));
2893 if ( ! in->selector_.empty() )
2894 throw(AipsError("Can only convert whole scantable. Unset the selection."));
2895 bool insitu = insitu_;
2896 setInsitu(false);
2897 CountedPtr< Scantable > out = getScantable(in, true);
2898 setInsitu(insitu);
2899 Table& tout = out->table();
2900 tout.rwKeywordSet().define("POLTYPE", String(newtype));
2901
2902 Block<String> cols(4);
2903 cols[0] = "SCANNO";
2904 cols[1] = "CYCLENO";
2905 cols[2] = "BEAMNO";
2906 cols[3] = "IFNO";
2907 TableIterator it(in->originalTable_, cols);
2908 String basetype = in->getPolType();
2909 STPol* stpol = STPol::getPolClass(in->factories_, basetype);
2910 try {
2911 while ( !it.pastEnd() ) {
2912 Table tab = it.table();
2913 uInt row = tab.rowNumbers()[0];
2914 stpol->setSpectra(in->getPolMatrix(row));
2915 Float fang,fhand;
2916 fang = in->focusTable_.getTotalAngle(in->mfocusidCol_(row));
2917 fhand = in->focusTable_.getFeedHand(in->mfocusidCol_(row));
2918 stpol->setPhaseCorrections(fang, fhand);
2919 Int npolout = 0;
2920 for (uInt i=0; i<tab.nrow(); ++i) {
2921 Vector<Float> outvec = stpol->getSpectrum(i, newtype);
2922 if ( outvec.nelements() > 0 ) {
2923 tout.addRow();
2924 TableCopy::copyRows(tout, tab, tout.nrow()-1, 0, 1);
2925 ArrayColumn<Float> sCol(tout,"SPECTRA");
2926 ScalarColumn<uInt> pCol(tout,"POLNO");
2927 sCol.put(tout.nrow()-1 ,outvec);
2928 pCol.put(tout.nrow()-1 ,uInt(npolout));
2929 npolout++;
2930 }
2931 }
2932 tout.rwKeywordSet().define("nPol", npolout);
2933 ++it;
2934 }
2935 } catch (AipsError& e) {
2936 delete stpol;
2937 throw(e);
2938 }
2939 delete stpol;
2940 return out;
2941}
2942
2943CountedPtr< Scantable >
2944 asap::STMath::mxExtract( const CountedPtr< Scantable > & in,
2945 const std::string & scantype )
2946{
2947 bool insitu = insitu_;
2948 setInsitu(false);
2949 CountedPtr< Scantable > out = getScantable(in, true);
2950 setInsitu(insitu);
2951 Table& tout = out->table();
2952 std::string taql = "SELECT FROM $1 WHERE BEAMNO != REFBEAMNO";
2953 if (scantype == "on") {
2954 taql = "SELECT FROM $1 WHERE BEAMNO == REFBEAMNO";
2955 }
2956 Table tab = tableCommand(taql, in->table());
2957 TableCopy::copyRows(tout, tab);
2958 if (scantype == "on") {
2959 // re-index SCANNO to 0
2960 TableVector<uInt> vec(tout, "SCANNO");
2961 vec = 0;
2962 }
2963 return out;
2964}
2965
2966std::vector<float>
2967 asap::STMath::fft( const casa::CountedPtr< Scantable > & in,
2968 const std::vector<int>& whichrow,
2969 bool getRealImag )
2970{
2971 std::vector<float> res;
2972 Table tab = in->table();
2973 std::vector<bool> mask;
2974
2975 if (whichrow.size() < 1) { // for all rows (by default)
2976 int nrow = int(tab.nrow());
2977 for (int i = 0; i < nrow; ++i) {
2978 res = in->execFFT(i, mask, getRealImag);
2979 }
2980 } else { // for specified rows
2981 for (uInt i = 0; i < whichrow.size(); ++i) {
2982 res = in->execFFT(i, mask, getRealImag);
2983 }
2984 }
2985
2986 return res;
2987}
2988
2989
2990CountedPtr<Scantable>
2991 asap::STMath::lagFlag( const CountedPtr<Scantable>& in,
2992 double start, double end,
2993 const std::string& mode )
2994{
2995 CountedPtr<Scantable> out = getScantable(in, false);
2996 Table& tout = out->table();
2997 TableIterator iter(tout, "FREQ_ID");
2998 FFTServer<Float,Complex> ffts;
2999
3000 while ( !iter.pastEnd() ) {
3001 Table tab = iter.table();
3002 Double rp,rv,inc;
3003 ROTableRow row(tab);
3004 const TableRecord& rec = row.get(0);
3005 uInt freqid = rec.asuInt("FREQ_ID");
3006 out->frequencies().getEntry(rp, rv, inc, freqid);
3007 ArrayColumn<Float> specCol(tab, "SPECTRA");
3008 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
3009
3010 for (int i=0; i<int(tab.nrow()); ++i) {
3011 Vector<Float> spec = specCol(i);
3012 Vector<uChar> flag = flagCol(i);
3013 std::vector<bool> mask;
3014 for (uInt j = 0; j < flag.nelements(); ++j) {
3015 mask.push_back(!(flag[j]>0));
3016 }
3017 mathutil::doZeroOrderInterpolation(spec, mask);
3018
3019 Vector<Complex> lags;
3020 ffts.fft0(lags, spec);
3021
3022 Int lag0(start+0.5);
3023 Int lag1(end+0.5);
3024 if (mode == "frequency") {
3025 lag0 = Int(spec.nelements()*abs(inc)/(start)+0.5);
3026 lag1 = Int(spec.nelements()*abs(inc)/(end)+0.5);
3027 }
3028 Int lstart = max(0, lag0);
3029 Int lend = min(Int(lags.nelements()-1), lag1);
3030 if (lstart == lend) {
3031 lags[lstart] = Complex(0.0);
3032 } else {
3033 if (lstart > lend) {
3034 Int tmp = lend;
3035 lend = lstart;
3036 lstart = tmp;
3037 }
3038 for (int j=lstart; j <=lend ;++j) {
3039 lags[j] = Complex(0.0);
3040 }
3041 }
3042
3043 ffts.fft0(spec, lags);
3044
3045 specCol.put(i, spec);
3046 }
3047 ++iter;
3048 }
3049 return out;
3050}
3051
3052// Averaging spectra with different channel/resolution
3053CountedPtr<Scantable>
3054STMath::new_average( const std::vector<CountedPtr<Scantable> >& in,
3055 const bool& compel,
3056 const std::vector<bool>& mask,
3057 const std::string& weight,
3058 const std::string& avmode )
3059 throw ( casa::AipsError )
3060{
3061 LogIO os( LogOrigin( "STMath", "new_average()", WHERE ) ) ;
3062 if ( avmode == "SCAN" && in.size() != 1 )
3063 throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
3064 "Use merge first."));
3065
3066 // 2012/02/17 TN
3067 // Since STGrid is implemented, average doesn't consider direction
3068 // when accumulating
3069 // check if OTF observation
3070// String obstype = in[0]->getHeader().obstype ;
3071// Double tol = 0.0 ;
3072// if ( obstype.find( "OTF" ) != String::npos ) {
3073// tol = TOL_OTF ;
3074// }
3075// else {
3076// tol = TOL_POINT ;
3077// }
3078
3079 CountedPtr<Scantable> out ; // processed result
3080 if ( compel ) {
3081 std::vector< CountedPtr<Scantable> > newin ; // input for average process
3082 uInt insize = in.size() ; // number of input scantables
3083
3084 // TEST: do normal average in each table before IF grouping
3085 os << "Do preliminary averaging" << LogIO::POST ;
3086 vector< CountedPtr<Scantable> > tmpin( insize ) ;
3087 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3088 vector< CountedPtr<Scantable> > v( 1, in[itable] ) ;
3089 tmpin[itable] = average( v, mask, weight, avmode ) ;
3090 }
3091
3092 // warning
3093 os << "Average spectra with different spectral resolution" << LogIO::POST ;
3094
3095 // temporarily set coordinfo
3096 vector<string> oldinfo( insize ) ;
3097 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3098 vector<string> coordinfo = in[itable]->getCoordInfo() ;
3099 oldinfo[itable] = coordinfo[0] ;
3100 coordinfo[0] = "Hz" ;
3101 tmpin[itable]->setCoordInfo( coordinfo ) ;
3102 }
3103
3104 // columns
3105 ScalarColumn<uInt> freqIDCol ;
3106 ScalarColumn<uInt> ifnoCol ;
3107 ScalarColumn<uInt> scannoCol ;
3108
3109
3110 // check IF frequency coverage
3111 // freqid: list of FREQ_ID, which is used, in each table
3112 // iffreq: list of minimum and maximum frequency for each FREQ_ID in
3113 // each table
3114 // freqid[insize][numIF]
3115 // freqid: [[id00, id01, ...],
3116 // [id10, id11, ...],
3117 // ...
3118 // [idn0, idn1, ...]]
3119 // iffreq[insize][numIF*2]
3120 // iffreq: [[min_id00, max_id00, min_id01, max_id01, ...],
3121 // [min_id10, max_id10, min_id11, max_id11, ...],
3122 // ...
3123 // [min_idn0, max_idn0, min_idn1, max_idn1, ...]]
3124 //os << "Check IF settings in each table" << LogIO::POST ;
3125 vector< vector<uInt> > freqid( insize );
3126 vector< vector<double> > iffreq( insize ) ;
3127 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3128 uInt rows = tmpin[itable]->nrow() ;
3129 uInt freqnrows = tmpin[itable]->frequencies().table().nrow() ;
3130 for ( uInt irow = 0 ; irow < rows ; irow++ ) {
3131 if ( freqid[itable].size() == freqnrows ) {
3132 break ;
3133 }
3134 else {
3135 freqIDCol.attach( tmpin[itable]->table(), "FREQ_ID" ) ;
3136 ifnoCol.attach( tmpin[itable]->table(), "IFNO" ) ;
3137 uInt id = freqIDCol( irow ) ;
3138 if ( freqid[itable].size() == 0 || count( freqid[itable].begin(), freqid[itable].end(), id ) == 0 ) {
3139 //os << "itable = " << itable << ": IF " << id << " is included in the list" << LogIO::POST ;
3140 vector<double> abcissa = tmpin[itable]->getAbcissa( irow ) ;
3141 freqid[itable].push_back( id ) ;
3142 iffreq[itable].push_back( abcissa[0] - 0.5 * ( abcissa[1] - abcissa[0] ) ) ;
3143 iffreq[itable].push_back( abcissa[abcissa.size()-1] + 0.5 * ( abcissa[1] - abcissa[0] ) ) ;
3144 }
3145 }
3146 }
3147 }
3148
3149 // debug
3150 //os << "IF settings summary:" << endl ;
3151 //for ( uInt i = 0 ; i < freqid.size() ; i++ ) {
3152 //os << " Table" << i << endl ;
3153 //for ( uInt j = 0 ; j < freqid[i].size() ; j++ ) {
3154 //os << " id = " << freqid[i][j] << " (min,max) = (" << iffreq[i][2*j] << "," << iffreq[i][2*j+1] << ")" << endl ;
3155 //}
3156 //}
3157 //os << endl ;
3158 //os.post() ;
3159
3160 // IF grouping based on their frequency coverage
3161 // ifgrp: list of table index and FREQ_ID for all members in each IF group
3162 // ifgfreq: list of minimum and maximum frequency in each IF group
3163 // ifgrp[numgrp][nummember*2]
3164 // ifgrp: [[table00, freqrow00, table01, freqrow01, ...],
3165 // [table10, freqrow10, table11, freqrow11, ...],
3166 // ...
3167 // [tablen0, freqrown0, tablen1, freqrown1, ...]]
3168 // ifgfreq[numgrp*2]
3169 // ifgfreq: [min0_grp0, max0_grp0, min1_grp1, max1_grp1, ...]
3170 //os << "IF grouping based on their frequency coverage" << LogIO::POST ;
3171 vector< vector<uInt> > ifgrp ;
3172 vector<double> ifgfreq ;
3173
3174 // parameter for IF grouping
3175 // groupmode = OR retrieve all region
3176 // AND only retrieve overlaped region
3177 //string groupmode = "AND" ;
3178 string groupmode = "OR" ;
3179 uInt sizecr = 0 ;
3180 if ( groupmode == "AND" )
3181 sizecr = 2 ;
3182 else if ( groupmode == "OR" )
3183 sizecr = 0 ;
3184
3185 vector<double> sortedfreq ;
3186 for ( uInt i = 0 ; i < iffreq.size() ; i++ ) {
3187 for ( uInt j = 0 ; j < iffreq[i].size() ; j++ ) {
3188 if ( count( sortedfreq.begin(), sortedfreq.end(), iffreq[i][j] ) == 0 )
3189 sortedfreq.push_back( iffreq[i][j] ) ;
3190 }
3191 }
3192 sort( sortedfreq.begin(), sortedfreq.end() ) ;
3193 for ( vector<double>::iterator i = sortedfreq.begin() ; i != sortedfreq.end()-1 ; i++ ) {
3194 ifgfreq.push_back( *i ) ;
3195 ifgfreq.push_back( *(i+1) ) ;
3196 }
3197 ifgrp.resize( ifgfreq.size()/2 ) ;
3198 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3199 for ( uInt iif = 0 ; iif < freqid[itable].size() ; iif++ ) {
3200 double range0 = iffreq[itable][2*iif] ;
3201 double range1 = iffreq[itable][2*iif+1] ;
3202 for ( uInt j = 0 ; j < ifgrp.size() ; j++ ) {
3203 double fmin = max( range0, ifgfreq[2*j] ) ;
3204 double fmax = min( range1, ifgfreq[2*j+1] ) ;
3205 if ( fmin < fmax ) {
3206 ifgrp[j].push_back( itable ) ;
3207 ifgrp[j].push_back( freqid[itable][iif] ) ;
3208 }
3209 }
3210 }
3211 }
3212 vector< vector<uInt> >::iterator fiter = ifgrp.begin() ;
3213 vector<double>::iterator giter = ifgfreq.begin() ;
3214 while( fiter != ifgrp.end() ) {
3215 if ( fiter->size() <= sizecr ) {
3216 fiter = ifgrp.erase( fiter ) ;
3217 giter = ifgfreq.erase( giter ) ;
3218 giter = ifgfreq.erase( giter ) ;
3219 }
3220 else {
3221 fiter++ ;
3222 advance( giter, 2 ) ;
3223 }
3224 }
3225
3226 // Grouping continuous IF groups (without frequency gap)
3227 // freqgrp: list of IF group indexes in each frequency group
3228 // freqrange: list of minimum and maximum frequency in each frequency group
3229 // freqgrp[numgrp][nummember]
3230 // freqgrp: [[ifgrp00, ifgrp01, ifgrp02, ...],
3231 // [ifgrp10, ifgrp11, ifgrp12, ...],
3232 // ...
3233 // [ifgrpn0, ifgrpn1, ifgrpn2, ...]]
3234 // freqrange[numgrp*2]
3235 // freqrange: [min_grp0, max_grp0, min_grp1, max_grp1, ...]
3236 vector< vector<uInt> > freqgrp ;
3237 double freqrange = 0.0 ;
3238 uInt grpnum = 0 ;
3239 for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
3240 // Assumed that ifgfreq was sorted
3241 if ( grpnum != 0 && freqrange == ifgfreq[2*i] ) {
3242 freqgrp[grpnum-1].push_back( i ) ;
3243 }
3244 else {
3245 vector<uInt> grp0( 1, i ) ;
3246 freqgrp.push_back( grp0 ) ;
3247 grpnum++ ;
3248 }
3249 freqrange = ifgfreq[2*i+1] ;
3250 }
3251
3252
3253 // print IF groups
3254 ostringstream oss ;
3255 oss << "IF Group summary: " << endl ;
3256 oss << " GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
3257 for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
3258 oss << " GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
3259 for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
3260 oss << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ;
3261 }
3262 oss << endl ;
3263 }
3264 oss << endl ;
3265 os << oss.str() << LogIO::POST ;
3266
3267 // print frequency group
3268 oss.str("") ;
3269 oss << "Frequency Group summary: " << endl ;
3270 oss << " GROUP_ID [FREQ_MIN, FREQ_MAX]: IF_GROUP_ID" << endl ;
3271 for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
3272 oss << " GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*freqgrp[i][0]] << "," << ifgfreq[2*freqgrp[i][freqgrp[i].size()-1]+1] << "]: " ;
3273 for ( uInt j = 0 ; j < freqgrp[i].size() ; j++ ) {
3274 oss << freqgrp[i][j] << " " ;
3275 }
3276 oss << endl ;
3277 }
3278 oss << endl ;
3279 os << oss.str() << LogIO::POST ;
3280
3281 // membership check
3282 // groups: list of IF group indexes whose frequency range overlaps with
3283 // that of each table and IF
3284 // groups[numtable][numIF][nummembership]
3285 // groups: [[[grp, grp,...], [grp, grp,...],...],
3286 // [[grp, grp,...], [grp, grp,...],...],
3287 // ...
3288 // [[grp, grp,...], [grp, grp,...],...]]
3289 vector< vector< vector<uInt> > > groups( insize ) ;
3290 for ( uInt i = 0 ; i < insize ; i++ ) {
3291 groups[i].resize( freqid[i].size() ) ;
3292 }
3293 for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
3294 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
3295 uInt tableid = ifgrp[igrp][2*imem] ;
3296 vector<uInt>::iterator iter = find( freqid[tableid].begin(), freqid[tableid].end(), ifgrp[igrp][2*imem+1] ) ;
3297 if ( iter != freqid[tableid].end() ) {
3298 uInt rowid = distance( freqid[tableid].begin(), iter ) ;
3299 groups[tableid][rowid].push_back( igrp ) ;
3300 }
3301 }
3302 }
3303
3304 // print membership
3305 //oss.str("") ;
3306 //for ( uInt i = 0 ; i < insize ; i++ ) {
3307 //oss << "Table " << i << endl ;
3308 //for ( uInt j = 0 ; j < groups[i].size() ; j++ ) {
3309 //oss << " FREQ_ID " << setw( 2 ) << freqid[i][j] << ": " ;
3310 //for ( uInt k = 0 ; k < groups[i][j].size() ; k++ ) {
3311 //oss << setw( 2 ) << groups[i][j][k] << " " ;
3312 //}
3313 //oss << endl ;
3314 //}
3315 //}
3316 //os << oss.str() << LogIO::POST ;
3317
3318 // set back coordinfo
3319 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3320 vector<string> coordinfo = tmpin[itable]->getCoordInfo() ;
3321 coordinfo[0] = oldinfo[itable] ;
3322 tmpin[itable]->setCoordInfo( coordinfo ) ;
3323 }
3324
3325 // Create additional table if needed
3326 bool oldInsitu = insitu_ ;
3327 setInsitu( false ) ;
3328 vector< vector<uInt> > addrow( insize ) ;
3329 vector<uInt> addtable( insize, 0 ) ;
3330 vector<uInt> newtableids( insize ) ;
3331 vector<uInt> newifids( insize, 0 ) ;
3332 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3333 //os << "Table " << itable << ": " ;
3334 for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
3335 addrow[itable].push_back( groups[itable][ifrow].size()-1 ) ;
3336 //os << addrow[itable][ifrow] << " " ;
3337 }
3338 addtable[itable] = *max_element( addrow[itable].begin(), addrow[itable].end() ) ;
3339 //os << "(" << addtable[itable] << ")" << LogIO::POST ;
3340 }
3341 newin.resize( insize ) ;
3342 copy( tmpin.begin(), tmpin.end(), newin.begin() ) ;
3343 for ( uInt i = 0 ; i < insize ; i++ ) {
3344 newtableids[i] = i ;
3345 }
3346 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3347 for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
3348 CountedPtr<Scantable> add = getScantable( newin[itable], false ) ;
3349 vector<int> freqidlist ;
3350 for ( uInt i = 0 ; i < groups[itable].size() ; i++ ) {
3351 if ( groups[itable][i].size() > iadd + 1 ) {
3352 freqidlist.push_back( freqid[itable][i] ) ;
3353 }
3354 }
3355 stringstream taqlstream ;
3356 taqlstream << "SELECT FROM $1 WHERE FREQ_ID IN [" ;
3357 for ( uInt i = 0 ; i < freqidlist.size() ; i++ ) {
3358 taqlstream << freqidlist[i] ;
3359 if ( i < freqidlist.size() - 1 )
3360 taqlstream << "," ;
3361 else
3362 taqlstream << "]" ;
3363 }
3364 string taql = taqlstream.str() ;
3365 //os << "taql = " << taql << LogIO::POST ;
3366 STSelector selector = STSelector() ;
3367 selector.setTaQL( taql ) ;
3368 add->setSelection( selector ) ;
3369 newin.push_back( add ) ;
3370 newtableids.push_back( itable ) ;
3371 newifids.push_back( iadd + 1 ) ;
3372 }
3373 }
3374
3375 // udpate ifgrp
3376 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3377 for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
3378 for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
3379 if ( groups[itable][ifrow].size() > iadd + 1 ) {
3380 uInt igrp = groups[itable][ifrow][iadd+1] ;
3381 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
3382 if ( ifgrp[igrp][2*imem] == newtableids[iadd+insize] && ifgrp[igrp][2*imem+1] == freqid[newtableids[iadd+insize]][ifrow] ) {
3383 ifgrp[igrp][2*imem] = insize + iadd ;
3384 }
3385 }
3386 }
3387 }
3388 }
3389 }
3390
3391 // print IF groups again for debug
3392 //oss.str( "" ) ;
3393 //oss << "IF Group summary: " << endl ;
3394 //oss << " GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
3395 //for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
3396 //oss << " GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
3397 //for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
3398 //oss << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ;
3399 //}
3400 //oss << endl ;
3401 //}
3402 //oss << endl ;
3403 //os << oss.str() << LogIO::POST ;
3404
3405 // reset SCANNO and IFNO/FREQ_ID: IF is reset by the result of sortation
3406 os << "All scan number is set to 0" << LogIO::POST ;
3407 //os << "All IF number is set to IF group index" << LogIO::POST ;
3408 insize = newin.size() ;
3409 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3410 uInt rows = newin[itable]->nrow() ;
3411 Table &tmpt = newin[itable]->table() ;
3412 freqIDCol.attach( tmpt, "FREQ_ID" ) ;
3413 scannoCol.attach( tmpt, "SCANNO" ) ;
3414 ifnoCol.attach( tmpt, "IFNO" ) ;
3415 for ( uInt irow=0 ; irow < rows ; irow++ ) {
3416 scannoCol.put( irow, 0 ) ;
3417 uInt freqID = freqIDCol( irow ) ;
3418 vector<uInt>::iterator iter = find( freqid[newtableids[itable]].begin(), freqid[newtableids[itable]].end(), freqID ) ;
3419 if ( iter != freqid[newtableids[itable]].end() ) {
3420 uInt index = distance( freqid[newtableids[itable]].begin(), iter ) ;
3421 ifnoCol.put( irow, groups[newtableids[itable]][index][newifids[itable]] ) ;
3422 }
3423 else {
3424 throw(AipsError("IF grouping was wrong in additional tables.")) ;
3425 }
3426 }
3427 }
3428 oldinfo.resize( insize ) ;
3429 setInsitu( oldInsitu ) ;
3430
3431 // temporarily set coordinfo
3432 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3433 vector<string> coordinfo = newin[itable]->getCoordInfo() ;
3434 oldinfo[itable] = coordinfo[0] ;
3435 coordinfo[0] = "Hz" ;
3436 newin[itable]->setCoordInfo( coordinfo ) ;
3437 }
3438
3439 // save column values in the vector
3440 vector< vector<uInt> > freqTableIdVec( insize ) ;
3441 vector< vector<uInt> > freqIdVec( insize ) ;
3442 vector< vector<uInt> > ifNoVec( insize ) ;
3443 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3444 ScalarColumn<uInt> freqIDs ;
3445 freqIDs.attach( newin[itable]->frequencies().table(), "ID" ) ;
3446 ifnoCol.attach( newin[itable]->table(), "IFNO" ) ;
3447 freqIDCol.attach( newin[itable]->table(), "FREQ_ID" ) ;
3448 for ( uInt irow = 0 ; irow < newin[itable]->frequencies().table().nrow() ; irow++ ) {
3449 freqTableIdVec[itable].push_back( freqIDs( irow ) ) ;
3450 }
3451 for ( uInt irow = 0 ; irow < newin[itable]->table().nrow() ; irow++ ) {
3452 freqIdVec[itable].push_back( freqIDCol( irow ) ) ;
3453 ifNoVec[itable].push_back( ifnoCol( irow ) ) ;
3454 }
3455 }
3456
3457 // reset spectra and flagtra: pick up common part of frequency coverage
3458 //os << "Pick common frequency range and align resolution" << LogIO::POST ;
3459 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3460 uInt rows = newin[itable]->nrow() ;
3461 int nminchan = -1 ;
3462 int nmaxchan = -1 ;
3463 vector<uInt> freqIdUpdate ;
3464 for ( uInt irow = 0 ; irow < rows ; irow++ ) {
3465 uInt ifno = ifNoVec[itable][irow] ; // IFNO is reset by group index
3466 double minfreq = ifgfreq[2*ifno] ;
3467 double maxfreq = ifgfreq[2*ifno+1] ;
3468 //os << "frequency range: [" << minfreq << "," << maxfreq << "]" << LogIO::POST ;
3469 vector<double> abcissa = newin[itable]->getAbcissa( irow ) ;
3470 int nchan = abcissa.size() ;
3471 double resol = abcissa[1] - abcissa[0] ;
3472 //os << "abcissa range : [" << abcissa[0] << "," << abcissa[nchan-1] << "]" << LogIO::POST ;
3473 if ( minfreq <= abcissa[0] )
3474 nminchan = 0 ;
3475 else {
3476 //double cfreq = ( minfreq - abcissa[0] ) / resol ;
3477 double cfreq = ( minfreq - abcissa[0] + 0.5 * resol ) / resol ;
3478 nminchan = int(cfreq) + ( ( cfreq - int(cfreq) <= 0.5 ) ? 0 : 1 ) ;
3479 }
3480 if ( maxfreq >= abcissa[abcissa.size()-1] )
3481 nmaxchan = abcissa.size() - 1 ;
3482 else {
3483 //double cfreq = ( abcissa[abcissa.size()-1] - maxfreq ) / resol ;
3484 double cfreq = ( abcissa[abcissa.size()-1] - maxfreq + 0.5 * resol ) / resol ;
3485 nmaxchan = abcissa.size() - 1 - int(cfreq) - ( ( cfreq - int(cfreq) >= 0.5 ) ? 1 : 0 ) ;
3486 }
3487 //os << "channel range (" << irow << "): [" << nminchan << "," << nmaxchan << "]" << LogIO::POST ;
3488 if ( nmaxchan > nminchan ) {
3489 newin[itable]->reshapeSpectrum( nminchan, nmaxchan, irow ) ;
3490 int newchan = nmaxchan - nminchan + 1 ;
3491 if ( count( freqIdUpdate.begin(), freqIdUpdate.end(), freqIdVec[itable][irow] ) == 0 && newchan < nchan )
3492 freqIdUpdate.push_back( freqIdVec[itable][irow] ) ;
3493 }
3494 else {
3495 throw(AipsError("Failed to pick up common part of frequency range.")) ;
3496 }
3497 }
3498 for ( uInt i = 0 ; i < freqIdUpdate.size() ; i++ ) {
3499 uInt freqId = freqIdUpdate[i] ;
3500 Double refpix ;
3501 Double refval ;
3502 Double increment ;
3503
3504 // update row
3505 newin[itable]->frequencies().getEntry( refpix, refval, increment, freqId ) ;
3506 refval = refval - ( refpix - nminchan ) * increment ;
3507 refpix = 0 ;
3508 newin[itable]->frequencies().setEntry( refpix, refval, increment, freqId ) ;
3509 }
3510 }
3511
3512
3513 // reset spectra and flagtra: align spectral resolution
3514 //os << "Align spectral resolution" << LogIO::POST ;
3515 // gmaxdnu: the coarsest frequency resolution in the frequency group
3516 // gmemid: member index that have a resolution equal to gmaxdnu
3517 // gmaxdnu[numfreqgrp]
3518 // gmaxdnu: [dnu0, dnu1, ...]
3519 // gmemid[numfreqgrp]
3520 // gmemid: [id0, id1, ...]
3521 vector<double> gmaxdnu( freqgrp.size(), 0.0 ) ;
3522 vector<uInt> gmemid( freqgrp.size(), 0 ) ;
3523 for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
3524 double maxdnu = 0.0 ; // maximum (coarsest) frequency resolution
3525 int minchan = INT_MAX ; // minimum channel number
3526 Double refpixref = -1 ; // reference of 'reference pixel'
3527 Double refvalref = -1 ; // reference of 'reference frequency'
3528 Double refinc = -1 ; // reference frequency resolution
3529 uInt refreqid ;
3530 uInt reftable = INT_MAX;
3531 // process only if group member > 1
3532 if ( ifgrp[igrp].size() > 2 ) {
3533 // find minchan and maxdnu in each group
3534 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
3535 uInt tableid = ifgrp[igrp][2*imem] ;
3536 uInt rowid = ifgrp[igrp][2*imem+1] ;
3537 vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
3538 if ( iter != freqIdVec[tableid].end() ) {
3539 uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
3540 vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
3541 int nchan = abcissa.size() ;
3542 double dnu = abcissa[1] - abcissa[0] ;
3543 //os << "GROUP " << igrp << " (" << tableid << "," << rowid << "): nchan = " << nchan << " (minchan = " << minchan << ")" << LogIO::POST ;
3544 if ( nchan < minchan ) {
3545 minchan = nchan ;
3546 maxdnu = dnu ;
3547 newin[tableid]->frequencies().getEntry( refpixref, refvalref, refinc, rowid ) ;
3548 refreqid = rowid ;
3549 reftable = tableid ;
3550 }
3551 }
3552 }
3553 // regrid spectra in each group
3554 os << "GROUP " << igrp << endl ;
3555 os << " Channel number is adjusted to " << minchan << endl ;
3556 os << " Corresponding frequency resolution is " << maxdnu << "Hz" << LogIO::POST ;
3557 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
3558 uInt tableid = ifgrp[igrp][2*imem] ;
3559 uInt rowid = ifgrp[igrp][2*imem+1] ;
3560 freqIDCol.attach( newin[tableid]->table(), "FREQ_ID" ) ;
3561 //os << "tableid = " << tableid << " rowid = " << rowid << ": " << LogIO::POST ;
3562 //os << " regridChannel applied to " ;
3563 //if ( tableid != reftable )
3564 refreqid = newin[tableid]->frequencies().addEntry( refpixref, refvalref, refinc ) ;
3565 for ( uInt irow = 0 ; irow < newin[tableid]->table().nrow() ; irow++ ) {
3566 uInt tfreqid = freqIdVec[tableid][irow] ;
3567 if ( tfreqid == rowid ) {
3568 //os << irow << " " ;
3569 newin[tableid]->regridChannel( minchan, maxdnu, irow ) ;
3570 freqIDCol.put( irow, refreqid ) ;
3571 freqIdVec[tableid][irow] = refreqid ;
3572 }
3573 }
3574 //os << LogIO::POST ;
3575 }
3576 }
3577 else {
3578 uInt tableid = ifgrp[igrp][0] ;
3579 uInt rowid = ifgrp[igrp][1] ;
3580 vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
3581 if ( iter != freqIdVec[tableid].end() ) {
3582 uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
3583 vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
3584 minchan = abcissa.size() ;
3585 maxdnu = abcissa[1] - abcissa[0] ;
3586 }
3587 }
3588 for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
3589 if ( count( freqgrp[i].begin(), freqgrp[i].end(), igrp ) > 0 ) {
3590 if ( maxdnu > gmaxdnu[i] ) {
3591 gmaxdnu[i] = maxdnu ;
3592 gmemid[i] = igrp ;
3593 }
3594 break ;
3595 }
3596 }
3597 }
3598
3599 // set back coordinfo
3600 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
3601 vector<string> coordinfo = newin[itable]->getCoordInfo() ;
3602 coordinfo[0] = oldinfo[itable] ;
3603 newin[itable]->setCoordInfo( coordinfo ) ;
3604 }
3605
3606 // accumulate all rows into the first table
3607 // NOTE: assumed in.size() = 1
3608 vector< CountedPtr<Scantable> > tmp( 1 ) ;
3609 if ( newin.size() == 1 )
3610 tmp[0] = newin[0] ;
3611 else
3612 tmp[0] = merge( newin ) ;
3613
3614 //return tmp[0] ;
3615
3616 // average
3617 CountedPtr<Scantable> tmpout = average( tmp, mask, weight, avmode ) ;
3618
3619 //return tmpout ;
3620
3621 // combine frequency group
3622 os << "Combine spectra based on frequency grouping" << LogIO::POST ;
3623 os << "IFNO is renumbered as frequency group ID (see above)" << LogIO::POST ;
3624 vector<string> coordinfo = tmpout->getCoordInfo() ;
3625 oldinfo[0] = coordinfo[0] ;
3626 coordinfo[0] = "Hz" ;
3627 tmpout->setCoordInfo( coordinfo ) ;
3628 // create proformas of output table
3629 stringstream taqlstream ;
3630 taqlstream << "SELECT FROM $1 WHERE IFNO IN [" ;
3631 for ( uInt i = 0 ; i < gmemid.size() ; i++ ) {
3632 taqlstream << gmemid[i] ;
3633 if ( i < gmemid.size() - 1 )
3634 taqlstream << "," ;
3635 else
3636 taqlstream << "]" ;
3637 }
3638 string taql = taqlstream.str() ;
3639 //os << "taql = " << taql << LogIO::POST ;
3640 STSelector selector = STSelector() ;
3641 selector.setTaQL( taql ) ;
3642 oldInsitu = insitu_ ;
3643 setInsitu( false ) ;
3644 out = getScantable( tmpout, false ) ;
3645 setInsitu( oldInsitu ) ;
3646 out->setSelection( selector ) ;
3647 // regrid rows
3648 ifnoCol.attach( tmpout->table(), "IFNO" ) ;
3649 for ( uInt irow = 0 ; irow < tmpout->table().nrow() ; irow++ ) {
3650 uInt ifno = ifnoCol( irow ) ;
3651 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3652 if ( count( freqgrp[igrp].begin(), freqgrp[igrp].end(), ifno ) > 0 ) {
3653 vector<double> abcissa = tmpout->getAbcissa( irow ) ;
3654 double bw = ( abcissa[1] - abcissa[0] ) * abcissa.size() ;
3655 int nchan = (int)( bw / gmaxdnu[igrp] ) ;
3656 tmpout->regridChannel( nchan, gmaxdnu[igrp], irow ) ;
3657 break ;
3658 }
3659 }
3660 }
3661 // combine spectra
3662 ArrayColumn<Float> specColOut ;
3663 specColOut.attach( out->table(), "SPECTRA" ) ;
3664 ArrayColumn<uChar> flagColOut ;
3665 flagColOut.attach( out->table(), "FLAGTRA" ) ;
3666 ScalarColumn<uInt> ifnoColOut ;
3667 ifnoColOut.attach( out->table(), "IFNO" ) ;
3668 ScalarColumn<uInt> polnoColOut ;
3669 polnoColOut.attach( out->table(), "POLNO" ) ;
3670 ScalarColumn<uInt> freqidColOut ;
3671 freqidColOut.attach( out->table(), "FREQ_ID" ) ;
3672// MDirection::ScalarColumn dirColOut ;
3673// dirColOut.attach( out->table(), "DIRECTION" ) ;
3674 Table &tab = tmpout->table() ;
3675 Block<String> cols(1);
3676 cols[0] = String("POLNO") ;
3677 TableIterator iter( tab, cols ) ;
3678 bool done = false ;
3679 vector< vector<uInt> > sizes( freqgrp.size() ) ;
3680 while( !iter.pastEnd() ) {
3681 vector< vector<Float> > specout( freqgrp.size() ) ;
3682 vector< vector<uChar> > flagout( freqgrp.size() ) ;
3683 ArrayColumn<Float> specCols ;
3684 specCols.attach( iter.table(), "SPECTRA" ) ;
3685 ArrayColumn<uChar> flagCols ;
3686 flagCols.attach( iter.table(), "FLAGTRA" ) ;
3687 ifnoCol.attach( iter.table(), "IFNO" ) ;
3688 ScalarColumn<uInt> polnos ;
3689 polnos.attach( iter.table(), "POLNO" ) ;
3690// MDirection::ScalarColumn dircol ;
3691// dircol.attach( iter.table(), "DIRECTION" ) ;
3692 uInt polno = polnos( 0 ) ;
3693 //os << "POLNO iteration: " << polno << LogIO::POST ;
3694// for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3695// sizes[igrp].resize( freqgrp[igrp].size() ) ;
3696// for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
3697// for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
3698// uInt ifno = ifnoCol( irow ) ;
3699// if ( ifno == freqgrp[igrp][imem] ) {
3700// Vector<Float> spec = specCols( irow ) ;
3701// Vector<uChar> flag = flagCols( irow ) ;
3702// vector<Float> svec ;
3703// spec.tovector( svec ) ;
3704// vector<uChar> fvec ;
3705// flag.tovector( fvec ) ;
3706// //os << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << LogIO::POST ;
3707// specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ;
3708// flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ;
3709// //os << "specout[" << igrp << "].size() = " << specout[igrp].size() << LogIO::POST ;
3710// sizes[igrp][imem] = spec.nelements() ;
3711// }
3712// }
3713// }
3714// for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
3715// uInt ifout = ifnoColOut( irow ) ;
3716// uInt polout = polnoColOut( irow ) ;
3717// if ( ifout == gmemid[igrp] && polout == polno ) {
3718// // set SPECTRA and FRAGTRA
3719// Vector<Float> newspec( specout[igrp] ) ;
3720// Vector<uChar> newflag( flagout[igrp] ) ;
3721// specColOut.put( irow, newspec ) ;
3722// flagColOut.put( irow, newflag ) ;
3723// // IFNO renumbering
3724// ifnoColOut.put( irow, igrp ) ;
3725// }
3726// }
3727// }
3728 // get a list of number of channels for each frequency group member
3729 if ( !done ) {
3730 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3731 sizes[igrp].resize( freqgrp[igrp].size() ) ;
3732 for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
3733 for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
3734 uInt ifno = ifnoCol( irow ) ;
3735 if ( ifno == freqgrp[igrp][imem] ) {
3736 Vector<Float> spec = specCols( irow ) ;
3737 sizes[igrp][imem] = spec.nelements() ;
3738 break ;
3739 }
3740 }
3741 }
3742 }
3743 done = true ;
3744 }
3745 // combine spectra
3746 for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
3747 uInt polout = polnoColOut( irow ) ;
3748 if ( polout == polno ) {
3749 uInt ifout = ifnoColOut( irow ) ;
3750// Vector<Double> direction = dirColOut(irow).getAngle(Unit(String("rad"))).getValue() ;
3751 uInt igrp ;
3752 for ( uInt jgrp = 0 ; jgrp < freqgrp.size() ; jgrp++ ) {
3753 if ( ifout == gmemid[jgrp] ) {
3754 igrp = jgrp ;
3755 break ;
3756 }
3757 }
3758 for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
3759 for ( uInt jrow = 0 ; jrow < iter.table().nrow() ; jrow++ ) {
3760 uInt ifno = ifnoCol( jrow ) ;
3761 // 2012/02/17 TN
3762 // Since STGrid is implemented, average doesn't consider direction
3763 // when accumulating
3764// Vector<Double> tdir = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
3765// //if ( ifno == freqgrp[igrp][imem] && allTrue( tdir == direction ) ) {
3766// Double dx = tdir[0] - direction[0] ;
3767// Double dy = tdir[1] - direction[1] ;
3768// Double dd = sqrt( dx * dx + dy * dy ) ;
3769 //if ( ifno == freqgrp[igrp][imem] && allNearAbs( tdir, direction, tol ) ) {
3770// if ( ifno == freqgrp[igrp][imem] && dd <= tol ) {
3771 if ( ifno == freqgrp[igrp][imem] ) {
3772 Vector<Float> spec = specCols( jrow ) ;
3773 Vector<uChar> flag = flagCols( jrow ) ;
3774 vector<Float> svec ;
3775 spec.tovector( svec ) ;
3776 vector<uChar> fvec ;
3777 flag.tovector( fvec ) ;
3778 //os << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << LogIO::POST ;
3779 specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ;
3780 flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ;
3781 //os << "specout[" << igrp << "].size() = " << specout[igrp].size() << LogIO::POST ;
3782 }
3783 }
3784 }
3785 // set SPECTRA and FRAGTRA
3786 Vector<Float> newspec( specout[igrp] ) ;
3787 Vector<uChar> newflag( flagout[igrp] ) ;
3788 specColOut.put( irow, newspec ) ;
3789 flagColOut.put( irow, newflag ) ;
3790 // IFNO renumbering
3791 ifnoColOut.put( irow, igrp ) ;
3792 }
3793 }
3794 iter++ ;
3795 }
3796 // update FREQUENCIES subtable
3797 vector<bool> updated( freqgrp.size(), false ) ;
3798 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3799 uInt index = 0 ;
3800 uInt pixShift = 0 ;
3801 while ( freqgrp[igrp][index] != gmemid[igrp] ) {
3802 pixShift += sizes[igrp][index++] ;
3803 }
3804 for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
3805 if ( ifnoColOut( irow ) == gmemid[igrp] && !updated[igrp] ) {
3806 uInt freqidOut = freqidColOut( irow ) ;
3807 //os << "freqgrp " << igrp << " freqidOut = " << freqidOut << LogIO::POST ;
3808 double refpix ;
3809 double refval ;
3810 double increm ;
3811 out->frequencies().getEntry( refpix, refval, increm, freqidOut ) ;
3812 refpix += pixShift ;
3813 out->frequencies().setEntry( refpix, refval, increm, freqidOut ) ;
3814 updated[igrp] = true ;
3815 }
3816 }
3817 }
3818
3819 //out = tmpout ;
3820
3821 coordinfo = tmpout->getCoordInfo() ;
3822 coordinfo[0] = oldinfo[0] ;
3823 tmpout->setCoordInfo( coordinfo ) ;
3824 }
3825 else {
3826 // simple average
3827 out = average( in, mask, weight, avmode ) ;
3828 }
3829
3830 return out;
3831}
3832
3833CountedPtr<Scantable> STMath::cwcal( const CountedPtr<Scantable>& s,
3834 const String calmode,
3835 const String antname )
3836{
3837 // frequency switch
3838 if ( calmode == "fs" ) {
3839 return cwcalfs( s, antname ) ;
3840 }
3841 else {
3842 vector<bool> masks = s->getMask( 0 ) ;
3843 vector<int> types ;
3844
3845 // sky scan
3846 STSelector sel = STSelector() ;
3847 types.push_back( SrcType::SKY ) ;
3848 sel.setTypes( types ) ;
3849 s->setSelection( sel ) ;
3850 vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ;
3851 CountedPtr<Scantable> asky = average( tmp, masks, "TINT", "SCAN" ) ;
3852 s->unsetSelection() ;
3853 sel.reset() ;
3854 types.clear() ;
3855
3856 // hot scan
3857 types.push_back( SrcType::HOT ) ;
3858 sel.setTypes( types ) ;
3859 s->setSelection( sel ) ;
3860 tmp.clear() ;
3861 tmp.push_back( getScantable( s, false ) ) ;
3862 CountedPtr<Scantable> ahot = average( tmp, masks, "TINT", "SCAN" ) ;
3863 s->unsetSelection() ;
3864 sel.reset() ;
3865 types.clear() ;
3866
3867 // cold scan
3868 CountedPtr<Scantable> acold ;
3869// types.push_back( SrcType::COLD ) ;
3870// sel.setTypes( types ) ;
3871// s->setSelection( sel ) ;
3872// tmp.clear() ;
3873// tmp.push_back( getScantable( s, false ) ) ;
3874// CountedPtr<Scantable> acold = average( tmp, masks, "TINT", "SCNAN" ) ;
3875// s->unsetSelection() ;
3876// sel.reset() ;
3877// types.clear() ;
3878
3879 // off scan
3880 types.push_back( SrcType::PSOFF ) ;
3881 sel.setTypes( types ) ;
3882 s->setSelection( sel ) ;
3883 tmp.clear() ;
3884 tmp.push_back( getScantable( s, false ) ) ;
3885 CountedPtr<Scantable> aoff = average( tmp, masks, "TINT", "SCAN" ) ;
3886 s->unsetSelection() ;
3887 sel.reset() ;
3888 types.clear() ;
3889
3890 // on scan
3891 bool insitu = insitu_ ;
3892 insitu_ = false ;
3893 CountedPtr<Scantable> out = getScantable( s, true ) ;
3894 insitu_ = insitu ;
3895 types.push_back( SrcType::PSON ) ;
3896 sel.setTypes( types ) ;
3897 s->setSelection( sel ) ;
3898 TableCopy::copyRows( out->table(), s->table() ) ;
3899 s->unsetSelection() ;
3900 sel.reset() ;
3901 types.clear() ;
3902
3903 // process each on scan
3904 ArrayColumn<Float> tsysCol ;
3905 tsysCol.attach( out->table(), "TSYS" ) ;
3906 for ( int i = 0 ; i < out->nrow() ; i++ ) {
3907 vector<float> sp = getCalibratedSpectra( out, aoff, asky, ahot, acold, i, antname ) ;
3908 out->setSpectrum( sp, i ) ;
3909 string reftime = out->getTime( i ) ;
3910 vector<int> ii( 1, out->getIF( i ) ) ;
3911 vector<int> ib( 1, out->getBeam( i ) ) ;
3912 vector<int> ip( 1, out->getPol( i ) ) ;
3913 sel.setIFs( ii ) ;
3914 sel.setBeams( ib ) ;
3915 sel.setPolarizations( ip ) ;
3916 asky->setSelection( sel ) ;
3917 vector<float> sptsys = getTsysFromTime( reftime, asky, "linear" ) ;
3918 const Vector<Float> Vtsys( sptsys ) ;
3919 tsysCol.put( i, Vtsys ) ;
3920 asky->unsetSelection() ;
3921 sel.reset() ;
3922 }
3923
3924 // flux unit
3925 out->setFluxUnit( "K" ) ;
3926
3927 return out ;
3928 }
3929}
3930
3931CountedPtr<Scantable> STMath::almacal( const CountedPtr<Scantable>& s,
3932 const String calmode )
3933{
3934 // frequency switch
3935 if ( calmode == "fs" ) {
3936 return almacalfs( s ) ;
3937 }
3938 else {
3939// double t0, t1 ;
3940 // t0 = mathutil::gettimeofday_sec() ;
3941 vector<bool> masks = s->getMask( 0 ) ;
3942
3943 // off scan
3944 STSelector sel = STSelector() ;
3945 vector<int> types( 1, SrcType::PSOFF ) ;
3946 //types.push_back( SrcType::PSOFF ) ;
3947 sel.setTypes( types ) ;
3948 s->setSelection( sel ) ;
3949 // TODO 2010/01/08 TN
3950 // Grouping by time should be needed before averaging.
3951 // Each group must have own unique SCANNO (should be renumbered).
3952 // See PIPELINE/SDCalibration.py
3953 CountedPtr<Scantable> soff = getScantable( s, false ) ;
3954 Table ttab = soff->table() ;
3955// String tmpname = File::newUniqueName( "./", "temp" ).baseName() ;
3956// Table ttab = s->table().copyToMemoryTable( tmpname, False ) ;
3957 ROScalarColumn<Double> *timeCol = new ROScalarColumn<Double>( ttab, "TIME" ) ;
3958 uInt nrow = timeCol->nrow() ;
3959 Vector<Double> timeSep = timeCol->getColumn() ;
3960 delete timeCol ;
3961 for ( uInt i = nrow-2 ; i > 0 ; i-- ) {
3962 timeSep[i] -= timeSep[i-1] ;
3963 }
3964 Vector<Double> interval = soff->integrCol_.getColumn() ;
3965 interval /= 86400.0 ;
3966 Block<uInt> gaplist( 100 ) ;
3967 uInt glidx = 0 ;
3968 uInt glsize = gaplist.size() ;
3969 for ( uInt i = 0 ; i < nrow - 1 ; i++ ) {
3970 double gap = 2.0 * timeSep[i+1] / ( interval[i] + interval[i+1] ) ;
3971 //cout << "gap[" << i << "]=" << setw(5) << gap << endl ;
3972 if ( gap > 1.1 ) {
3973// cout << "detected gap between " << i << " and " << i+1 << endl ;
3974 gaplist[glidx] = i ;
3975 glidx++ ;
3976 }
3977 if ( glidx >= glsize ) {
3978// cout << "resize gaplist" << endl ;
3979 glsize += 100 ;
3980 gaplist.resize( glsize ) ;
3981 }
3982 }
3983 gaplist[glidx] = nrow - 1 ;
3984 glidx++ ;
3985// cout << "gaplist = " << Vector<uInt>(IPosition(1,glidx),gaplist.storage()) << endl ;
3986 uInt newid = 0 ;
3987 Vector<uInt> newscanno( nrow, 0 ) ;
3988 uInt *p = newscanno.data() ;
3989 IPosition pos( 1 ) ;
3990 for ( uInt i = 1 ; i < glidx ; i++ ) {
3991 pos[0] = gaplist[i] - gaplist[i-1] ;
3992 Vector<uInt> scnslice( pos, p+gaplist[i-1]+1, SHARE ) ;
3993 scnslice = i ;
3994 }
3995 soff->scanCol_.putColumn( newscanno ) ;
3996// cout << "new scancol = " << soff->scanCol_.getColumn() << endl ;
3997// double t2 = mathutil::gettimeofday_sec() ;
3998 vector< CountedPtr<Scantable> > tmp( 1, soff ) ;
3999 CountedPtr<Scantable> aoff = average( tmp, masks, "TINT", "SCAN" ) ;
4000 // double t3 = mathutil::gettimeofday_sec() ;
4001 //cout << "aoff.nrow = " << aoff->nrow() << endl ;
4002 s->unsetSelection() ;
4003 sel.reset() ;
4004 // t1 = mathutil::gettimeofday_sec() ;
4005 // cout << "elapsed time for off averaging: " << t1-t0 << " sec" << endl ;
4006 // cout << " elapsed time for average(): " << t3-t2 << " sec" << endl ;
4007
4008 // on scan
4009// t0 = mathutil::gettimeofday_sec() ;
4010 bool insitu = insitu_ ;
4011 insitu_ = false ;
4012 CountedPtr<Scantable> out = getScantable( s, true ) ;
4013 insitu_ = insitu ;
4014 types[0] = SrcType::PSON ;
4015 sel.setTypes( types ) ;
4016 s->setSelection( sel ) ;
4017 //TableCopy::copyRows( out->table(), s->table() ) ;
4018 out->table().addRow( s->nrow() ) ;
4019 copyRows( out->table(), s->table(), 0, 0, s->nrow(), False ) ;
4020 //s->unsetSelection() ;
4021 sel.reset() ;
4022// t1 = mathutil::gettimeofday_sec() ;
4023// cout << "elapsed time for preparing output table: " << t1-t0 << " sec" << endl ;
4024
4025 // process each on scan
4026 // t0 = mathutil::gettimeofday_sec() ;
4027// for ( int i = 0 ; i < out->nrow() ; i++ ) {
4028// vector<float> sp = getCalibratedSpectra( out, aoff, i ) ;
4029// out->setSpectrum( sp, i ) ;
4030// }
4031
4032 // using STIdxIterAcc
4033 vector<string> cols( 3 ) ;
4034 cols[0] = "BEAMNO" ;
4035 cols[1] = "POLNO" ;
4036 cols[2] = "IFNO" ;
4037 STIdxIter *iter = new STIdxIterAcc( out, cols ) ;
4038 while ( !iter->pastEnd() ) {
4039 Vector<uInt> ids = iter->current() ;
4040 stringstream ss ;
4041 ss << "SELECT FROM $1 WHERE "
4042 << "BEAMNO==" << ids[0] << "&&"
4043 << "POLNO==" << ids[1] << "&&"
4044 << "IFNO==" << ids[2] ;
4045 //cout << "TaQL string: " << ss.str() << endl ;
4046 sel.setTaQL( ss.str() ) ;
4047 aoff->setSelection( sel ) ;
4048 Vector<uInt> rows = iter->getRows( SHARE ) ;
4049 // out should be an exact copy of s except that SPECTRA column is empty
4050 calibrateALMA( out, s, aoff, rows ) ;
4051 aoff->unsetSelection() ;
4052 sel.reset() ;
4053 iter->next() ;
4054 }
4055 delete iter ;
4056 s->unsetSelection() ;
4057
4058 // t1 = mathutil::gettimeofday_sec() ;
4059 // cout << "elapsed time for calibration: " << t1-t0 << " sec" << endl ;
4060
4061 // flux unit
4062 out->setFluxUnit( "K" ) ;
4063
4064 return out ;
4065 }
4066}
4067
4068CountedPtr<Scantable> STMath::cwcalfs( const CountedPtr<Scantable>& s,
4069 const String antname )
4070{
4071 vector<int> types ;
4072
4073 // APEX calibration mode
4074 int apexcalmode = 1 ;
4075
4076 if ( antname.find( "APEX" ) != string::npos ) {
4077 // check if off scan exists or not
4078 STSelector sel = STSelector() ;
4079 //sel.setName( offstr1 ) ;
4080 types.push_back( SrcType::FLOOFF ) ;
4081 sel.setTypes( types ) ;
4082 try {
4083 s->setSelection( sel ) ;
4084 }
4085 catch ( AipsError &e ) {
4086 apexcalmode = 0 ;
4087 }
4088 sel.reset() ;
4089 }
4090 s->unsetSelection() ;
4091 types.clear() ;
4092
4093 vector<bool> masks = s->getMask( 0 ) ;
4094 CountedPtr<Scantable> ssig, sref ;
4095 CountedPtr<Scantable> out ;
4096
4097 if ( antname.find( "APEX" ) != string::npos ) {
4098 // APEX calibration
4099 // sky scan
4100 STSelector sel = STSelector() ;
4101 types.push_back( SrcType::FLOSKY ) ;
4102 sel.setTypes( types ) ;
4103 s->setSelection( sel ) ;
4104 vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ;
4105 CountedPtr<Scantable> askylo = average( tmp, masks, "TINT", "SCAN" ) ;
4106 s->unsetSelection() ;
4107 sel.reset() ;
4108 types.clear() ;
4109 types.push_back( SrcType::FHISKY ) ;
4110 sel.setTypes( types ) ;
4111 s->setSelection( sel ) ;
4112 tmp.clear() ;
4113 tmp.push_back( getScantable( s, false ) ) ;
4114 CountedPtr<Scantable> askyhi = average( tmp, masks, "TINT", "SCAN" ) ;
4115 s->unsetSelection() ;
4116 sel.reset() ;
4117 types.clear() ;
4118
4119 // hot scan
4120 types.push_back( SrcType::FLOHOT ) ;
4121 sel.setTypes( types ) ;
4122 s->setSelection( sel ) ;
4123 tmp.clear() ;
4124 tmp.push_back( getScantable( s, false ) ) ;
4125 CountedPtr<Scantable> ahotlo = average( tmp, masks, "TINT", "SCAN" ) ;
4126 s->unsetSelection() ;
4127 sel.reset() ;
4128 types.clear() ;
4129 types.push_back( SrcType::FHIHOT ) ;
4130 sel.setTypes( types ) ;
4131 s->setSelection( sel ) ;
4132 tmp.clear() ;
4133 tmp.push_back( getScantable( s, false ) ) ;
4134 CountedPtr<Scantable> ahothi = average( tmp, masks, "TINT", "SCAN" ) ;
4135 s->unsetSelection() ;
4136 sel.reset() ;
4137 types.clear() ;
4138
4139 // cold scan
4140 CountedPtr<Scantable> acoldlo, acoldhi ;
4141// types.push_back( SrcType::FLOCOLD ) ;
4142// sel.setTypes( types ) ;
4143// s->setSelection( sel ) ;
4144// tmp.clear() ;
4145// tmp.push_back( getScantable( s, false ) ) ;
4146// CountedPtr<Scantable> acoldlo = average( tmp, masks, "TINT", "SCAN" ) ;
4147// s->unsetSelection() ;
4148// sel.reset() ;
4149// types.clear() ;
4150// types.push_back( SrcType::FHICOLD ) ;
4151// sel.setTypes( types ) ;
4152// s->setSelection( sel ) ;
4153// tmp.clear() ;
4154// tmp.push_back( getScantable( s, false ) ) ;
4155// CountedPtr<Scantable> acoldhi = average( tmp, masks, "TINT", "SCAN" ) ;
4156// s->unsetSelection() ;
4157// sel.reset() ;
4158// types.clear() ;
4159
4160 // ref scan
4161 bool insitu = insitu_ ;
4162 insitu_ = false ;
4163 sref = getScantable( s, true ) ;
4164 insitu_ = insitu ;
4165 types.push_back( SrcType::FSLO ) ;
4166 sel.setTypes( types ) ;
4167 s->setSelection( sel ) ;
4168 TableCopy::copyRows( sref->table(), s->table() ) ;
4169 s->unsetSelection() ;
4170 sel.reset() ;
4171 types.clear() ;
4172
4173 // sig scan
4174 insitu_ = false ;
4175 ssig = getScantable( s, true ) ;
4176 insitu_ = insitu ;
4177 types.push_back( SrcType::FSHI ) ;
4178 sel.setTypes( types ) ;
4179 s->setSelection( sel ) ;
4180 TableCopy::copyRows( ssig->table(), s->table() ) ;
4181 s->unsetSelection() ;
4182 sel.reset() ;
4183 types.clear() ;
4184
4185 if ( apexcalmode == 0 ) {
4186 // APEX fs data without off scan
4187 // process each sig and ref scan
4188 ArrayColumn<Float> tsysCollo ;
4189 tsysCollo.attach( ssig->table(), "TSYS" ) ;
4190 ArrayColumn<Float> tsysColhi ;
4191 tsysColhi.attach( sref->table(), "TSYS" ) ;
4192 for ( int i = 0 ; i < ssig->nrow() ; i++ ) {
4193 vector< CountedPtr<Scantable> > sky( 2 ) ;
4194 sky[0] = askylo ;
4195 sky[1] = askyhi ;
4196 vector< CountedPtr<Scantable> > hot( 2 ) ;
4197 hot[0] = ahotlo ;
4198 hot[1] = ahothi ;
4199 vector< CountedPtr<Scantable> > cold( 2 ) ;
4200 //cold[0] = acoldlo ;
4201 //cold[1] = acoldhi ;
4202 vector<float> sp = getFSCalibratedSpectra( ssig, sref, sky, hot, cold, i ) ;
4203 ssig->setSpectrum( sp, i ) ;
4204 string reftime = ssig->getTime( i ) ;
4205 vector<int> ii( 1, ssig->getIF( i ) ) ;
4206 vector<int> ib( 1, ssig->getBeam( i ) ) ;
4207 vector<int> ip( 1, ssig->getPol( i ) ) ;
4208 sel.setIFs( ii ) ;
4209 sel.setBeams( ib ) ;
4210 sel.setPolarizations( ip ) ;
4211 askylo->setSelection( sel ) ;
4212 vector<float> sptsys = getTsysFromTime( reftime, askylo, "linear" ) ;
4213 const Vector<Float> Vtsyslo( sptsys ) ;
4214 tsysCollo.put( i, Vtsyslo ) ;
4215 askylo->unsetSelection() ;
4216 sel.reset() ;
4217 sky[0] = askyhi ;
4218 sky[1] = askylo ;
4219 hot[0] = ahothi ;
4220 hot[1] = ahotlo ;
4221 cold[0] = acoldhi ;
4222 cold[1] = acoldlo ;
4223 sp = getFSCalibratedSpectra( sref, ssig, sky, hot, cold, i ) ;
4224 sref->setSpectrum( sp, i ) ;
4225 reftime = sref->getTime( i ) ;
4226 ii[0] = sref->getIF( i ) ;
4227 ib[0] = sref->getBeam( i ) ;
4228 ip[0] = sref->getPol( i ) ;
4229 sel.setIFs( ii ) ;
4230 sel.setBeams( ib ) ;
4231 sel.setPolarizations( ip ) ;
4232 askyhi->setSelection( sel ) ;
4233 sptsys = getTsysFromTime( reftime, askyhi, "linear" ) ;
4234 const Vector<Float> Vtsyshi( sptsys ) ;
4235 tsysColhi.put( i, Vtsyshi ) ;
4236 askyhi->unsetSelection() ;
4237 sel.reset() ;
4238 }
4239 }
4240 else if ( apexcalmode == 1 ) {
4241 // APEX fs data with off scan
4242 // off scan
4243 types.push_back( SrcType::FLOOFF ) ;
4244 sel.setTypes( types ) ;
4245 s->setSelection( sel ) ;
4246 tmp.clear() ;
4247 tmp.push_back( getScantable( s, false ) ) ;
4248 CountedPtr<Scantable> aofflo = average( tmp, masks, "TINT", "SCAN" ) ;
4249 s->unsetSelection() ;
4250 sel.reset() ;
4251 types.clear() ;
4252 types.push_back( SrcType::FHIOFF ) ;
4253 sel.setTypes( types ) ;
4254 s->setSelection( sel ) ;
4255 tmp.clear() ;
4256 tmp.push_back( getScantable( s, false ) ) ;
4257 CountedPtr<Scantable> aoffhi = average( tmp, masks, "TINT", "SCAN" ) ;
4258 s->unsetSelection() ;
4259 sel.reset() ;
4260 types.clear() ;
4261
4262 // process each sig and ref scan
4263 ArrayColumn<Float> tsysCollo ;
4264 tsysCollo.attach( ssig->table(), "TSYS" ) ;
4265 ArrayColumn<Float> tsysColhi ;
4266 tsysColhi.attach( sref->table(), "TSYS" ) ;
4267 for ( int i = 0 ; i < ssig->nrow() ; i++ ) {
4268 vector<float> sp = getCalibratedSpectra( ssig, aofflo, askylo, ahotlo, acoldlo, i, antname ) ;
4269 ssig->setSpectrum( sp, i ) ;
4270 sp = getCalibratedSpectra( sref, aoffhi, askyhi, ahothi, acoldhi, i, antname ) ;
4271 string reftime = ssig->getTime( i ) ;
4272 vector<int> ii( 1, ssig->getIF( i ) ) ;
4273 vector<int> ib( 1, ssig->getBeam( i ) ) ;
4274 vector<int> ip( 1, ssig->getPol( i ) ) ;
4275 sel.setIFs( ii ) ;
4276 sel.setBeams( ib ) ;
4277 sel.setPolarizations( ip ) ;
4278 askylo->setSelection( sel ) ;
4279 vector<float> sptsys = getTsysFromTime( reftime, askylo, "linear" ) ;
4280 const Vector<Float> Vtsyslo( sptsys ) ;
4281 tsysCollo.put( i, Vtsyslo ) ;
4282 askylo->unsetSelection() ;
4283 sel.reset() ;
4284 sref->setSpectrum( sp, i ) ;
4285 reftime = sref->getTime( i ) ;
4286 ii[0] = sref->getIF( i ) ;
4287 ib[0] = sref->getBeam( i ) ;
4288 ip[0] = sref->getPol( i ) ;
4289 sel.setIFs( ii ) ;
4290 sel.setBeams( ib ) ;
4291 sel.setPolarizations( ip ) ;
4292 askyhi->setSelection( sel ) ;
4293 sptsys = getTsysFromTime( reftime, askyhi, "linear" ) ;
4294 const Vector<Float> Vtsyshi( sptsys ) ;
4295 tsysColhi.put( i, Vtsyshi ) ;
4296 askyhi->unsetSelection() ;
4297 sel.reset() ;
4298 }
4299 }
4300 }
4301 else {
4302 // non-APEX fs data
4303 // sky scan
4304 STSelector sel = STSelector() ;
4305 types.push_back( SrcType::SKY ) ;
4306 sel.setTypes( types ) ;
4307 s->setSelection( sel ) ;
4308 vector< CountedPtr<Scantable> > tmp( 1, getScantable( s, false ) ) ;
4309 CountedPtr<Scantable> asky = average( tmp, masks, "TINT", "SCAN" ) ;
4310 s->unsetSelection() ;
4311 sel.reset() ;
4312 types.clear() ;
4313
4314 // hot scan
4315 types.push_back( SrcType::HOT ) ;
4316 sel.setTypes( types ) ;
4317 s->setSelection( sel ) ;
4318 tmp.clear() ;
4319 tmp.push_back( getScantable( s, false ) ) ;
4320 CountedPtr<Scantable> ahot = average( tmp, masks, "TINT", "SCAN" ) ;
4321 s->unsetSelection() ;
4322 sel.reset() ;
4323 types.clear() ;
4324
4325 // cold scan
4326 CountedPtr<Scantable> acold ;
4327// types.push_back( SrcType::COLD ) ;
4328// sel.setTypes( types ) ;
4329// s->setSelection( sel ) ;
4330// tmp.clear() ;
4331// tmp.push_back( getScantable( s, false ) ) ;
4332// CountedPtr<Scantable> acold = average( tmp, masks, "TINT", "SCAN" ) ;
4333// s->unsetSelection() ;
4334// sel.reset() ;
4335// types.clear() ;
4336
4337 // ref scan
4338 bool insitu = insitu_ ;
4339 insitu_ = false ;
4340 sref = getScantable( s, true ) ;
4341 insitu_ = insitu ;
4342 types.push_back( SrcType::FSOFF ) ;
4343 sel.setTypes( types ) ;
4344 s->setSelection( sel ) ;
4345 TableCopy::copyRows( sref->table(), s->table() ) ;
4346 s->unsetSelection() ;
4347 sel.reset() ;
4348 types.clear() ;
4349
4350 // sig scan
4351 insitu_ = false ;
4352 ssig = getScantable( s, true ) ;
4353 insitu_ = insitu ;
4354 types.push_back( SrcType::FSON ) ;
4355 sel.setTypes( types ) ;
4356 s->setSelection( sel ) ;
4357 TableCopy::copyRows( ssig->table(), s->table() ) ;
4358 s->unsetSelection() ;
4359 sel.reset() ;
4360 types.clear() ;
4361
4362 // process each sig and ref scan
4363 ArrayColumn<Float> tsysColsig ;
4364 tsysColsig.attach( ssig->table(), "TSYS" ) ;
4365 ArrayColumn<Float> tsysColref ;
4366 tsysColref.attach( ssig->table(), "TSYS" ) ;
4367 for ( int i = 0 ; i < ssig->nrow() ; i++ ) {
4368 vector<float> sp = getFSCalibratedSpectra( ssig, sref, asky, ahot, acold, i ) ;
4369 ssig->setSpectrum( sp, i ) ;
4370 string reftime = ssig->getTime( i ) ;
4371 vector<int> ii( 1, ssig->getIF( i ) ) ;
4372 vector<int> ib( 1, ssig->getBeam( i ) ) ;
4373 vector<int> ip( 1, ssig->getPol( i ) ) ;
4374 sel.setIFs( ii ) ;
4375 sel.setBeams( ib ) ;
4376 sel.setPolarizations( ip ) ;
4377 asky->setSelection( sel ) ;
4378 vector<float> sptsys = getTsysFromTime( reftime, asky, "linear" ) ;
4379 const Vector<Float> Vtsys( sptsys ) ;
4380 tsysColsig.put( i, Vtsys ) ;
4381 asky->unsetSelection() ;
4382 sel.reset() ;
4383 sp = getFSCalibratedSpectra( sref, ssig, asky, ahot, acold, i ) ;
4384 sref->setSpectrum( sp, i ) ;
4385 tsysColref.put( i, Vtsys ) ;
4386 }
4387 }
4388
4389 // do folding if necessary
4390 Table sigtab = ssig->table() ;
4391 Table reftab = sref->table() ;
4392 ScalarColumn<uInt> sigifnoCol ;
4393 ScalarColumn<uInt> refifnoCol ;
4394 ScalarColumn<uInt> sigfidCol ;
4395 ScalarColumn<uInt> reffidCol ;
4396 Int nchan = (Int)ssig->nchan() ;
4397 sigifnoCol.attach( sigtab, "IFNO" ) ;
4398 refifnoCol.attach( reftab, "IFNO" ) ;
4399 sigfidCol.attach( sigtab, "FREQ_ID" ) ;
4400 reffidCol.attach( reftab, "FREQ_ID" ) ;
4401 Vector<uInt> sfids( sigfidCol.getColumn() ) ;
4402 Vector<uInt> rfids( reffidCol.getColumn() ) ;
4403 vector<uInt> sfids_unique ;
4404 vector<uInt> rfids_unique ;
4405 vector<uInt> sifno_unique ;
4406 vector<uInt> rifno_unique ;
4407 for ( uInt i = 0 ; i < sfids.nelements() ; i++ ) {
4408 if ( count( sfids_unique.begin(), sfids_unique.end(), sfids[i] ) == 0 ) {
4409 sfids_unique.push_back( sfids[i] ) ;
4410 sifno_unique.push_back( ssig->getIF( i ) ) ;
4411 }
4412 if ( count( rfids_unique.begin(), rfids_unique.end(), rfids[i] ) == 0 ) {
4413 rfids_unique.push_back( rfids[i] ) ;
4414 rifno_unique.push_back( sref->getIF( i ) ) ;
4415 }
4416 }
4417 double refpix_sig, refval_sig, increment_sig ;
4418 double refpix_ref, refval_ref, increment_ref ;
4419 vector< CountedPtr<Scantable> > tmp( sfids_unique.size() ) ;
4420 for ( uInt i = 0 ; i < sfids_unique.size() ; i++ ) {
4421 ssig->frequencies().getEntry( refpix_sig, refval_sig, increment_sig, sfids_unique[i] ) ;
4422 sref->frequencies().getEntry( refpix_ref, refval_ref, increment_ref, rfids_unique[i] ) ;
4423 if ( refpix_sig == refpix_ref ) {
4424 double foffset = refval_ref - refval_sig ;
4425 int choffset = static_cast<int>(foffset/increment_sig) ;
4426 double doffset = foffset / increment_sig ;
4427 if ( abs(choffset) >= nchan ) {
4428 LogIO os( LogOrigin( "STMath", "cwcalfs", WHERE ) ) ;
4429 os << "FREQ_ID=[" << sfids_unique[i] << "," << rfids_unique[i] << "]: out-band frequency switching, no folding" << LogIO::POST ;
4430 os << "Just return signal data" << LogIO::POST ;
4431 //std::vector< CountedPtr<Scantable> > tabs ;
4432 //tabs.push_back( ssig ) ;
4433 //tabs.push_back( sref ) ;
4434 //out = merge( tabs ) ;
4435 tmp[i] = ssig ;
4436 }
4437 else {
4438 STSelector sel = STSelector() ;
4439 vector<int> v( 1, sifno_unique[i] ) ;
4440 sel.setIFs( v ) ;
4441 ssig->setSelection( sel ) ;
4442 sel.reset() ;
4443 v[0] = rifno_unique[i] ;
4444 sel.setIFs( v ) ;
4445 sref->setSelection( sel ) ;
4446 sel.reset() ;
4447 if ( antname.find( "APEX" ) != string::npos ) {
4448 tmp[i] = dofold( ssig, sref, 0.5*doffset, -0.5*doffset ) ;
4449 //tmp[i] = dofold( ssig, sref, doffset ) ;
4450 }
4451 else {
4452 tmp[i] = dofold( ssig, sref, doffset ) ;
4453 }
4454 ssig->unsetSelection() ;
4455 sref->unsetSelection() ;
4456 }
4457 }
4458 }
4459
4460 if ( tmp.size() > 1 ) {
4461 out = merge( tmp ) ;
4462 }
4463 else {
4464 out = tmp[0] ;
4465 }
4466
4467 // flux unit
4468 out->setFluxUnit( "K" ) ;
4469
4470 return out ;
4471}
4472
4473CountedPtr<Scantable> STMath::almacalfs( const CountedPtr<Scantable>& s )
4474{
4475 (void) s; //currently unused
4476 CountedPtr<Scantable> out ;
4477
4478 return out ;
4479}
4480
4481vector<float> STMath::getSpectrumFromTime( string reftime,
4482 CountedPtr<Scantable>& s,
4483 string mode )
4484{
4485 LogIO os( LogOrigin( "STMath", "getSpectrumFromTime", WHERE ) ) ;
4486 vector<float> sp ;
4487
4488 if ( s->nrow() == 0 ) {
4489 os << LogIO::SEVERE << "No spectra in the input scantable. Return empty spectrum." << LogIO::POST ;
4490 return sp ;
4491 }
4492 else if ( s->nrow() == 1 ) {
4493 //os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ;
4494 return s->getSpectrum( 0 ) ;
4495 }
4496 else {
4497 vector<int> idx = getRowIdFromTime( reftime, s ) ;
4498 if ( mode == "before" ) {
4499 int id = -1 ;
4500 if ( idx[0] != -1 ) {
4501 id = idx[0] ;
4502 }
4503 else if ( idx[1] != -1 ) {
4504 os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
4505 id = idx[1] ;
4506 }
4507 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4508 sp = s->getSpectrum( id ) ;
4509 }
4510 else if ( mode == "after" ) {
4511 int id = -1 ;
4512 if ( idx[1] != -1 ) {
4513 id = idx[1] ;
4514 }
4515 else if ( idx[0] != -1 ) {
4516 os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
4517 id = idx[1] ;
4518 }
4519 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4520 sp = s->getSpectrum( id ) ;
4521 }
4522 else if ( mode == "nearest" ) {
4523 int id = -1 ;
4524 if ( idx[0] == -1 ) {
4525 id = idx[1] ;
4526 }
4527 else if ( idx[1] == -1 ) {
4528 id = idx[0] ;
4529 }
4530 else if ( idx[0] == idx[1] ) {
4531 id = idx[0] ;
4532 }
4533 else {
4534 //double t0 = getMJD( s->getTime( idx[0] ) ) ;
4535 //double t1 = getMJD( s->getTime( idx[1] ) ) ;
4536 double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
4537 double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
4538 double tref = getMJD( reftime ) ;
4539 if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
4540 id = idx[1] ;
4541 }
4542 else {
4543 id = idx[0] ;
4544 }
4545 }
4546 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4547 sp = s->getSpectrum( id ) ;
4548 }
4549 else if ( mode == "linear" ) {
4550 if ( idx[0] == -1 ) {
4551 // use after
4552 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
4553 int id = idx[1] ;
4554 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4555 sp = s->getSpectrum( id ) ;
4556 }
4557 else if ( idx[1] == -1 ) {
4558 // use before
4559 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
4560 int id = idx[0] ;
4561 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4562 sp = s->getSpectrum( id ) ;
4563 }
4564 else if ( idx[0] == idx[1] ) {
4565 // use before
4566 //os << "No need to interporate." << LogIO::POST ;
4567 int id = idx[0] ;
4568 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4569 sp = s->getSpectrum( id ) ;
4570 }
4571 else {
4572 // do interpolation
4573 //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
4574 //double t0 = getMJD( s->getTime( idx[0] ) ) ;
4575 //double t1 = getMJD( s->getTime( idx[1] ) ) ;
4576 double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
4577 double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
4578 double tref = getMJD( reftime ) ;
4579 vector<float> sp0 = s->getSpectrum( idx[0] ) ;
4580 vector<float> sp1 = s->getSpectrum( idx[1] ) ;
4581 for ( unsigned int i = 0 ; i < sp0.size() ; i++ ) {
4582 float v = ( sp1[i] - sp0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + sp0[i] ;
4583 sp.push_back( v ) ;
4584 }
4585 }
4586 }
4587 else {
4588 os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
4589 }
4590 return sp ;
4591 }
4592}
4593
4594Vector<Float> STMath::getSpectrumFromTime( double reftime,
4595 Vector<Double> &timeVec,
4596 const CountedPtr<Scantable>& s,
4597 string mode )
4598{
4599 LogIO os( LogOrigin( "STMath", "getSpectrumFromTime", WHERE ) ) ;
4600 Vector<Float> sp ;
4601
4602 if ( s->nrow() == 0 ) {
4603 os << LogIO::SEVERE << "No spectra in the input scantable. Return empty spectrum." << LogIO::POST ;
4604 return sp ;
4605 }
4606 else if ( s->nrow() == 1 ) {
4607 //os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ;
4608 s->specCol_.get( 0, sp ) ;
4609 return sp ;
4610 }
4611 else {
4612 vector<int> idx = getRowIdFromTime( reftime, timeVec ) ;
4613 if ( mode == "before" ) {
4614 int id = -1 ;
4615 if ( idx[0] != -1 ) {
4616 id = idx[0] ;
4617 }
4618 else if ( idx[1] != -1 ) {
4619 os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
4620 id = idx[1] ;
4621 }
4622 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4623 s->specCol_.get( id, sp ) ;
4624 //sp = s->getSpectrum( id ) ;
4625 }
4626 else if ( mode == "after" ) {
4627 int id = -1 ;
4628 if ( idx[1] != -1 ) {
4629 id = idx[1] ;
4630 }
4631 else if ( idx[0] != -1 ) {
4632 os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
4633 id = idx[1] ;
4634 }
4635 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4636 s->specCol_.get( id, sp ) ;
4637 //sp = s->getSpectrum( id ) ;
4638 }
4639 else if ( mode == "nearest" ) {
4640 int id = -1 ;
4641 if ( idx[0] == -1 ) {
4642 id = idx[1] ;
4643 }
4644 else if ( idx[1] == -1 ) {
4645 id = idx[0] ;
4646 }
4647 else if ( idx[0] == idx[1] ) {
4648 id = idx[0] ;
4649 }
4650 else {
4651 double t0 = timeVec[idx[0]] ;
4652 double t1 = timeVec[idx[1]] ;
4653 double tref = reftime ;
4654 if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
4655 id = idx[1] ;
4656 }
4657 else {
4658 id = idx[0] ;
4659 }
4660 }
4661 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4662 s->specCol_.get( id, sp ) ;
4663 //sp = s->getSpectrum( id ) ;
4664 }
4665 else if ( mode == "linear" ) {
4666 if ( idx[0] == -1 ) {
4667 // use after
4668 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
4669 int id = idx[1] ;
4670 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4671 s->specCol_.get( id, sp ) ;
4672 //sp = s->getSpectrum( id ) ;
4673 }
4674 else if ( idx[1] == -1 ) {
4675 // use before
4676 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
4677 int id = idx[0] ;
4678 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4679 s->specCol_.get( id, sp ) ;
4680 //sp = s->getSpectrum( id ) ;
4681 }
4682 else if ( idx[0] == idx[1] ) {
4683 // use before
4684 //os << "No need to interporate." << LogIO::POST ;
4685 int id = idx[0] ;
4686 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4687 s->specCol_.get( id, sp ) ;
4688 //sp = s->getSpectrum( id ) ;
4689 }
4690 else {
4691 // do interpolation
4692 //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
4693 double t0 = timeVec[idx[0]] ;
4694 double t1 = timeVec[idx[1]] ;
4695 double tref = reftime ;
4696 s->specCol_.get( idx[0], sp ) ;
4697 //sp = s->getSpectrum( idx[0] ) ;
4698 //vector<float> sp1 = s->getSpectrum( idx[1] ) ;
4699 Vector<Float> sp1 = s->specCol_( idx[1] ) ;
4700 double tfactor = ( tref - t0 ) / ( t1 - t0 ) ;
4701 for ( unsigned int i = 0 ; i < sp.size() ; i++ ) {
4702 sp[i] = ( sp1[i] - sp[i] ) * tfactor + sp[i] ;
4703 }
4704 }
4705 }
4706 else {
4707 os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
4708 }
4709 return sp ;
4710 }
4711}
4712
4713Vector<Float> STMath::getSpectrumFromTime( double reftime,
4714 const Vector<Double> &timeVec,
4715 const Matrix<Float>& spectra,
4716 string mode )
4717{
4718 LogIO os( LogOrigin( "STMath", "getSpectrumFromTime", WHERE ) ) ;
4719 Vector<Float> sp ;
4720 uInt ncol = spectra.ncolumn() ;
4721
4722 if ( ncol == 0 ) {
4723 os << LogIO::SEVERE << "No spectra in the input scantable. Return empty spectrum." << LogIO::POST ;
4724 return sp ;
4725 }
4726 else if ( ncol == 1 ) {
4727 //os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ;
4728 //s->specCol_.get( 0, sp ) ;
4729 sp.reference( spectra.column( 0 ) ) ;
4730 return sp ;
4731 }
4732 else {
4733 vector<int> idx = getRowIdFromTime( reftime, timeVec ) ;
4734 if ( mode == "before" ) {
4735 int id = -1 ;
4736 if ( idx[0] != -1 ) {
4737 id = idx[0] ;
4738 }
4739 else if ( idx[1] != -1 ) {
4740 os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
4741 id = idx[1] ;
4742 }
4743 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4744 //s->specCol_.get( id, sp ) ;
4745 //sp = s->getSpectrum( id ) ;
4746 sp.reference( spectra.column( id ) ) ;
4747 }
4748 else if ( mode == "after" ) {
4749 int id = -1 ;
4750 if ( idx[1] != -1 ) {
4751 id = idx[1] ;
4752 }
4753 else if ( idx[0] != -1 ) {
4754 os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
4755 id = idx[1] ;
4756 }
4757 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4758 //s->specCol_.get( id, sp ) ;
4759 //sp = s->getSpectrum( id ) ;
4760 sp.reference( spectra.column( id ) ) ;
4761 }
4762 else if ( mode == "nearest" ) {
4763 int id = -1 ;
4764 if ( idx[0] == -1 ) {
4765 id = idx[1] ;
4766 }
4767 else if ( idx[1] == -1 ) {
4768 id = idx[0] ;
4769 }
4770 else if ( idx[0] == idx[1] ) {
4771 id = idx[0] ;
4772 }
4773 else {
4774 double t0 = timeVec[idx[0]] ;
4775 double t1 = timeVec[idx[1]] ;
4776 double tref = reftime ;
4777 if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
4778 id = idx[1] ;
4779 }
4780 else {
4781 id = idx[0] ;
4782 }
4783 }
4784 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4785 //s->specCol_.get( id, sp ) ;
4786 //sp = s->getSpectrum( id ) ;
4787 sp.reference( spectra.column( id ) ) ;
4788 }
4789 else if ( mode == "linear" ) {
4790 if ( idx[0] == -1 ) {
4791 // use after
4792 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
4793 int id = idx[1] ;
4794 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4795 //s->specCol_.get( id, sp ) ;
4796 //sp = s->getSpectrum( id ) ;
4797 sp.reference( spectra.column( id ) ) ;
4798 }
4799 else if ( idx[1] == -1 ) {
4800 // use before
4801 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
4802 int id = idx[0] ;
4803 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4804 //s->specCol_.get( id, sp ) ;
4805 //sp = s->getSpectrum( id ) ;
4806 sp.reference( spectra.column( id ) ) ;
4807 }
4808 else if ( idx[0] == idx[1] ) {
4809 // use before
4810 //os << "No need to interporate." << LogIO::POST ;
4811 int id = idx[0] ;
4812 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4813 //s->specCol_.get( id, sp ) ;
4814 //sp = s->getSpectrum( id ) ;
4815 sp.reference( spectra.column( id ) ) ;
4816 }
4817 else {
4818 // do interpolation
4819 //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
4820 double t0 = timeVec[idx[0]] ;
4821 double t1 = timeVec[idx[1]] ;
4822 double tref = reftime ;
4823 //s->specCol_.get( idx[0], sp ) ;
4824 //sp = s->getSpectrum( idx[0] ) ;
4825 //vector<float> sp1 = s->getSpectrum( idx[1] ) ;
4826 sp = spectra.column( idx[0] ).copy() ;
4827 //Vector<Float> sp1 = s->specCol_( idx[1] ) ;
4828 Vector<Float> sp1( spectra.column( idx[1] ) ) ;
4829 double tfactor = ( tref - t0 ) / ( t1 - t0 ) ;
4830 for ( unsigned int i = 0 ; i < sp.size() ; i++ ) {
4831 sp[i] = ( sp1[i] - sp[i] ) * tfactor + sp[i] ;
4832 }
4833 }
4834 }
4835 else {
4836 os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
4837 }
4838 return sp ;
4839 }
4840}
4841
4842void STMath::getSpectrumFromTime2( Vector<Float> &sp,
4843 double reftime,
4844 Vector<Double> &timeVec,
4845 const CountedPtr<Scantable>& s,
4846 string mode )
4847{
4848 LogIO os( LogOrigin( "STMath", "getSpectrumFromTime", WHERE ) ) ;
4849 //Vector<Float> sp ;
4850
4851 if ( s->nrow() == 0 ) {
4852 os << LogIO::SEVERE << "No spectra in the input scantable. Return empty spectrum." << LogIO::POST ;
4853 //return sp ;
4854 }
4855 else if ( s->nrow() == 1 ) {
4856 //os << "use row " << 0 << " (scanno = " << s->getScan( 0 ) << ")" << LogIO::POST ;
4857 //return s->getSpectrum( 0 ) ;
4858 s->specCol_.get( 0, sp ) ;
4859 //return sp ;
4860 }
4861 else {
4862 vector<int> idx = getRowIdFromTime( reftime, timeVec ) ;
4863 if ( mode == "before" ) {
4864 int id = -1 ;
4865 if ( idx[0] != -1 ) {
4866 id = idx[0] ;
4867 }
4868 else if ( idx[1] != -1 ) {
4869 os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
4870 id = idx[1] ;
4871 }
4872 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4873 //sp = s->getSpectrum( id ) ;
4874 s->specCol_.get( id, sp ) ;
4875 }
4876 else if ( mode == "after" ) {
4877 int id = -1 ;
4878 if ( idx[1] != -1 ) {
4879 id = idx[1] ;
4880 }
4881 else if ( idx[0] != -1 ) {
4882 os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
4883 id = idx[1] ;
4884 }
4885 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4886 //sp = s->getSpectrum( id ) ;
4887 s->specCol_.get( id, sp ) ;
4888 }
4889 else if ( mode == "nearest" ) {
4890 int id = -1 ;
4891 if ( idx[0] == -1 ) {
4892 id = idx[1] ;
4893 }
4894 else if ( idx[1] == -1 ) {
4895 id = idx[0] ;
4896 }
4897 else if ( idx[0] == idx[1] ) {
4898 id = idx[0] ;
4899 }
4900 else {
4901 double t0 = timeVec[idx[0]] ;
4902 double t1 = timeVec[idx[1]] ;
4903 double tref = reftime ;
4904 if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
4905 id = idx[1] ;
4906 }
4907 else {
4908 id = idx[0] ;
4909 }
4910 }
4911 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4912 //sp = s->getSpectrum( id ) ;
4913 s->specCol_.get( id, sp ) ;
4914 }
4915 else if ( mode == "linear" ) {
4916 if ( idx[0] == -1 ) {
4917 // use after
4918 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
4919 int id = idx[1] ;
4920 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4921 //sp = s->getSpectrum( id ) ;
4922 s->specCol_.get( id, sp ) ;
4923 }
4924 else if ( idx[1] == -1 ) {
4925 // use before
4926 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
4927 int id = idx[0] ;
4928 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4929 //sp = s->getSpectrum( id ) ;
4930 s->specCol_.get( id, sp ) ;
4931 }
4932 else if ( idx[0] == idx[1] ) {
4933 // use before
4934 //os << "No need to interporate." << LogIO::POST ;
4935 int id = idx[0] ;
4936 //os << "use row " << id << " (scanno = " << s->getScan( id ) << ")" << LogIO::POST ;
4937 //sp = s->getSpectrum( id ) ;
4938 s->specCol_.get( id, sp ) ;
4939 }
4940 else {
4941 // do interpolation
4942 //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
4943 //double t0 = timeVec[idx[0]] ;
4944 //double t1 = timeVec[idx[1]] ;
4945 //double tref = reftime ;
4946 //sp = s->getSpectrum( idx[0] ) ;
4947 s->specCol_.get( idx[0], sp ) ;
4948 //vector<float> sp1 = s->getSpectrum( idx[1] ) ;
4949 Vector<Float> sp1 = s->specCol_( idx[1] ) ;
4950 //double tfactor = ( tref - t0 ) / ( t1 - t0 ) ;
4951 double tfactor = ( reftime - timeVec[idx[0]] ) / ( timeVec[idx[1]] - timeVec[idx[0]] ) ;
4952 for ( unsigned int i = 0 ; i < sp.size() ; i++ ) {
4953 sp[i] = ( sp1[i] - sp[i] ) * tfactor + sp[i] ;
4954 }
4955 }
4956 }
4957 else {
4958 os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
4959 }
4960 //return sp ;
4961 }
4962}
4963
4964double STMath::getMJD( string strtime )
4965{
4966 if ( strtime.find("/") == string::npos ) {
4967 // MJD time string
4968 return atof( strtime.c_str() ) ;
4969 }
4970 else {
4971 // string in YYYY/MM/DD/HH:MM:SS format
4972 uInt year = atoi( strtime.substr( 0, 4 ).c_str() ) ;
4973 uInt month = atoi( strtime.substr( 5, 2 ).c_str() ) ;
4974 uInt day = atoi( strtime.substr( 8, 2 ).c_str() ) ;
4975 uInt hour = atoi( strtime.substr( 11, 2 ).c_str() ) ;
4976 uInt minute = atoi( strtime.substr( 14, 2 ).c_str() ) ;
4977 uInt sec = atoi( strtime.substr( 17, 2 ).c_str() ) ;
4978 Time t( year, month, day, hour, minute, sec ) ;
4979 return t.modifiedJulianDay() ;
4980 }
4981}
4982
4983vector<int> STMath::getRowIdFromTime( string reftime, CountedPtr<Scantable> &s )
4984{
4985 double reft = getMJD( reftime ) ;
4986 double dtmin = 1.0e100 ;
4987 double dtmax = -1.0e100 ;
4988 vector<double> dt ;
4989 int just_before = -1 ;
4990 int just_after = -1 ;
4991 for ( int i = 0 ; i < s->nrow() ; i++ ) {
4992 dt.push_back( getMJD( s->getTime( i ) ) - reft ) ;
4993 }
4994 for ( unsigned int i = 0 ; i < dt.size() ; i++ ) {
4995 if ( dt[i] > 0.0 ) {
4996 // after reftime
4997 if ( dt[i] < dtmin ) {
4998 just_after = i ;
4999 dtmin = dt[i] ;
5000 }
5001 }
5002 else if ( dt[i] < 0.0 ) {
5003 // before reftime
5004 if ( dt[i] > dtmax ) {
5005 just_before = i ;
5006 dtmax = dt[i] ;
5007 }
5008 }
5009 else {
5010 // just a reftime
5011 just_before = i ;
5012 just_after = i ;
5013 dtmax = 0 ;
5014 dtmin = 0 ;
5015 break ;
5016 }
5017 }
5018
5019 vector<int> v ;
5020 v.push_back( just_before ) ;
5021 v.push_back( just_after ) ;
5022
5023 return v ;
5024}
5025
5026vector<int> STMath::getRowIdFromTime( double reftime, const Vector<Double> &t )
5027{
5028// double reft = reftime ;
5029 double dtmin = 1.0e100 ;
5030 double dtmax = -1.0e100 ;
5031// vector<double> dt ;
5032 int just_before = -1 ;
5033 int just_after = -1 ;
5034 //cout << setprecision(24) << reft << endl ;
5035// ROScalarColumn<Double> timeCol( s->table(), "TIME" ) ;
5036// for ( int i = 0 ; i < s->nrow() ; i++ ) {
5037// cout << setprecision(24) << timeCol(i) << endl ;
5038// //dt.push_back( getMJD( s->getTime( i ) ) - reft ) ;
5039// dt.push_back( timeCol(i) - reft ) ;
5040// }
5041 Vector<Double> dt = t - reftime ;
5042 for ( unsigned int i = 0 ; i < dt.size() ; i++ ) {
5043 if ( dt[i] > 0.0 ) {
5044 // after reftime
5045 if ( dt[i] < dtmin ) {
5046 just_after = i ;
5047 dtmin = dt[i] ;
5048 }
5049 }
5050 else if ( dt[i] < 0.0 ) {
5051 // before reftime
5052 if ( dt[i] > dtmax ) {
5053 just_before = i ;
5054 dtmax = dt[i] ;
5055 }
5056 }
5057 else {
5058 // just a reftime
5059 just_before = i ;
5060 just_after = i ;
5061 dtmax = 0 ;
5062 dtmin = 0 ;
5063 break ;
5064 }
5065 }
5066
5067 vector<int> v(2) ;
5068 v[0] = just_before ;
5069 v[1] = just_after ;
5070
5071 return v ;
5072}
5073
5074vector<float> STMath::getTcalFromTime( string reftime,
5075 CountedPtr<Scantable>& s,
5076 string mode )
5077{
5078 LogIO os( LogOrigin( "STMath", "getTcalFromTime", WHERE ) ) ;
5079 vector<float> tcal ;
5080 STTcal tcalTable = s->tcal() ;
5081 String time ;
5082 Vector<Float> tcalval ;
5083 if ( s->nrow() == 0 ) {
5084 os << LogIO::SEVERE << "No row in the input scantable. Return empty tcal." << LogIO::POST ;
5085 return tcal ;
5086 }
5087 else if ( s->nrow() == 1 ) {
5088 uInt tcalid = s->getTcalId( 0 ) ;
5089 //os << "use row " << 0 << " (tcalid = " << tcalid << ")" << LogIO::POST ;
5090 tcalTable.getEntry( time, tcalval, tcalid ) ;
5091 tcalval.tovector( tcal ) ;
5092 return tcal ;
5093 }
5094 else {
5095 vector<int> idx = getRowIdFromTime( reftime, s ) ;
5096 if ( mode == "before" ) {
5097 int id = -1 ;
5098 if ( idx[0] != -1 ) {
5099 id = idx[0] ;
5100 }
5101 else if ( idx[1] != -1 ) {
5102 os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
5103 id = idx[1] ;
5104 }
5105 uInt tcalid = s->getTcalId( id ) ;
5106 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
5107 tcalTable.getEntry( time, tcalval, tcalid ) ;
5108 tcalval.tovector( tcal ) ;
5109 }
5110 else if ( mode == "after" ) {
5111 int id = -1 ;
5112 if ( idx[1] != -1 ) {
5113 id = idx[1] ;
5114 }
5115 else if ( idx[0] != -1 ) {
5116 os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
5117 id = idx[1] ;
5118 }
5119 uInt tcalid = s->getTcalId( id ) ;
5120 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
5121 tcalTable.getEntry( time, tcalval, tcalid ) ;
5122 tcalval.tovector( tcal ) ;
5123 }
5124 else if ( mode == "nearest" ) {
5125 int id = -1 ;
5126 if ( idx[0] == -1 ) {
5127 id = idx[1] ;
5128 }
5129 else if ( idx[1] == -1 ) {
5130 id = idx[0] ;
5131 }
5132 else if ( idx[0] == idx[1] ) {
5133 id = idx[0] ;
5134 }
5135 else {
5136 //double t0 = getMJD( s->getTime( idx[0] ) ) ;
5137 //double t1 = getMJD( s->getTime( idx[1] ) ) ;
5138 double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
5139 double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
5140 double tref = getMJD( reftime ) ;
5141 if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
5142 id = idx[1] ;
5143 }
5144 else {
5145 id = idx[0] ;
5146 }
5147 }
5148 uInt tcalid = s->getTcalId( id ) ;
5149 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
5150 tcalTable.getEntry( time, tcalval, tcalid ) ;
5151 tcalval.tovector( tcal ) ;
5152 }
5153 else if ( mode == "linear" ) {
5154 if ( idx[0] == -1 ) {
5155 // use after
5156 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
5157 int id = idx[1] ;
5158 uInt tcalid = s->getTcalId( id ) ;
5159 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
5160 tcalTable.getEntry( time, tcalval, tcalid ) ;
5161 tcalval.tovector( tcal ) ;
5162 }
5163 else if ( idx[1] == -1 ) {
5164 // use before
5165 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
5166 int id = idx[0] ;
5167 uInt tcalid = s->getTcalId( id ) ;
5168 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
5169 tcalTable.getEntry( time, tcalval, tcalid ) ;
5170 tcalval.tovector( tcal ) ;
5171 }
5172 else if ( idx[0] == idx[1] ) {
5173 // use before
5174 //os << "No need to interporate." << LogIO::POST ;
5175 int id = idx[0] ;
5176 uInt tcalid = s->getTcalId( id ) ;
5177 //os << "use row " << id << " (tcalid = " << tcalid << ")" << LogIO::POST ;
5178 tcalTable.getEntry( time, tcalval, tcalid ) ;
5179 tcalval.tovector( tcal ) ;
5180 }
5181 else {
5182 // do interpolation
5183 //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
5184 //double t0 = getMJD( s->getTime( idx[0] ) ) ;
5185 //double t1 = getMJD( s->getTime( idx[1] ) ) ;
5186 double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
5187 double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
5188 double tref = getMJD( reftime ) ;
5189 vector<float> tcal0 ;
5190 vector<float> tcal1 ;
5191 uInt tcalid0 = s->getTcalId( idx[0] ) ;
5192 uInt tcalid1 = s->getTcalId( idx[1] ) ;
5193 tcalTable.getEntry( time, tcalval, tcalid0 ) ;
5194 tcalval.tovector( tcal0 ) ;
5195 tcalTable.getEntry( time, tcalval, tcalid1 ) ;
5196 tcalval.tovector( tcal1 ) ;
5197 for ( unsigned int i = 0 ; i < tcal0.size() ; i++ ) {
5198 float v = ( tcal1[i] - tcal0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + tcal0[i] ;
5199 tcal.push_back( v ) ;
5200 }
5201 }
5202 }
5203 else {
5204 os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
5205 }
5206 return tcal ;
5207 }
5208}
5209
5210vector<float> STMath::getTsysFromTime( string reftime,
5211 CountedPtr<Scantable>& s,
5212 string mode )
5213{
5214 LogIO os( LogOrigin( "STMath", "getTsysFromTime", WHERE ) ) ;
5215 ArrayColumn<Float> tsysCol ;
5216 tsysCol.attach( s->table(), "TSYS" ) ;
5217 vector<float> tsys ;
5218 String time ;
5219 Vector<Float> tsysval ;
5220 if ( s->nrow() == 0 ) {
5221 os << LogIO::SEVERE << "No row in the input scantable. Return empty tsys." << LogIO::POST ;
5222 return tsys ;
5223 }
5224 else if ( s->nrow() == 1 ) {
5225 //os << "use row " << 0 << LogIO::POST ;
5226 tsysval = tsysCol( 0 ) ;
5227 tsysval.tovector( tsys ) ;
5228 return tsys ;
5229 }
5230 else {
5231 vector<int> idx = getRowIdFromTime( reftime, s ) ;
5232 if ( mode == "before" ) {
5233 int id = -1 ;
5234 if ( idx[0] != -1 ) {
5235 id = idx[0] ;
5236 }
5237 else if ( idx[1] != -1 ) {
5238 os << LogIO::WARN << "Failed to find a scan before reftime. return a spectrum just after the reftime." << LogIO::POST ;
5239 id = idx[1] ;
5240 }
5241 //os << "use row " << id << LogIO::POST ;
5242 tsysval = tsysCol( id ) ;
5243 tsysval.tovector( tsys ) ;
5244 }
5245 else if ( mode == "after" ) {
5246 int id = -1 ;
5247 if ( idx[1] != -1 ) {
5248 id = idx[1] ;
5249 }
5250 else if ( idx[0] != -1 ) {
5251 os << LogIO::WARN << "Failed to find a scan after reftime. return a spectrum just before the reftime." << LogIO::POST ;
5252 id = idx[1] ;
5253 }
5254 //os << "use row " << id << LogIO::POST ;
5255 tsysval = tsysCol( id ) ;
5256 tsysval.tovector( tsys ) ;
5257 }
5258 else if ( mode == "nearest" ) {
5259 int id = -1 ;
5260 if ( idx[0] == -1 ) {
5261 id = idx[1] ;
5262 }
5263 else if ( idx[1] == -1 ) {
5264 id = idx[0] ;
5265 }
5266 else if ( idx[0] == idx[1] ) {
5267 id = idx[0] ;
5268 }
5269 else {
5270 //double t0 = getMJD( s->getTime( idx[0] ) ) ;
5271 //double t1 = getMJD( s->getTime( idx[1] ) ) ;
5272 double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
5273 double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
5274 double tref = getMJD( reftime ) ;
5275 if ( abs( t0 - tref ) > abs( t1 - tref ) ) {
5276 id = idx[1] ;
5277 }
5278 else {
5279 id = idx[0] ;
5280 }
5281 }
5282 //os << "use row " << id << LogIO::POST ;
5283 tsysval = tsysCol( id ) ;
5284 tsysval.tovector( tsys ) ;
5285 }
5286 else if ( mode == "linear" ) {
5287 if ( idx[0] == -1 ) {
5288 // use after
5289 os << LogIO::WARN << "Failed to interpolate. return a spectrum just after the reftime." << LogIO::POST ;
5290 int id = idx[1] ;
5291 //os << "use row " << id << LogIO::POST ;
5292 tsysval = tsysCol( id ) ;
5293 tsysval.tovector( tsys ) ;
5294 }
5295 else if ( idx[1] == -1 ) {
5296 // use before
5297 os << LogIO::WARN << "Failed to interpolate. return a spectrum just before the reftime." << LogIO::POST ;
5298 int id = idx[0] ;
5299 //os << "use row " << id << LogIO::POST ;
5300 tsysval = tsysCol( id ) ;
5301 tsysval.tovector( tsys ) ;
5302 }
5303 else if ( idx[0] == idx[1] ) {
5304 // use before
5305 //os << "No need to interporate." << LogIO::POST ;
5306 int id = idx[0] ;
5307 //os << "use row " << id << LogIO::POST ;
5308 tsysval = tsysCol( id ) ;
5309 tsysval.tovector( tsys ) ;
5310 }
5311 else {
5312 // do interpolation
5313 //os << "interpolate between " << idx[0] << " and " << idx[1] << " (scanno: " << s->getScan( idx[0] ) << ", " << s->getScan( idx[1] ) << ")" << LogIO::POST ;
5314 //double t0 = getMJD( s->getTime( idx[0] ) ) ;
5315 //double t1 = getMJD( s->getTime( idx[1] ) ) ;
5316 double t0 = s->getEpoch( idx[0] ).get( Unit( "d" ) ).getValue() ;
5317 double t1 = s->getEpoch( idx[1] ).get( Unit( "d" ) ).getValue() ;
5318 double tref = getMJD( reftime ) ;
5319 vector<float> tsys0 ;
5320 vector<float> tsys1 ;
5321 tsysval = tsysCol( idx[0] ) ;
5322 tsysval.tovector( tsys0 ) ;
5323 tsysval = tsysCol( idx[1] ) ;
5324 tsysval.tovector( tsys1 ) ;
5325 for ( unsigned int i = 0 ; i < tsys0.size() ; i++ ) {
5326 float v = ( tsys1[i] - tsys0[i] ) / ( t1 - t0 ) * ( tref - t0 ) + tsys0[i] ;
5327 tsys.push_back( v ) ;
5328 }
5329 }
5330 }
5331 else {
5332 os << LogIO::SEVERE << "Unknown mode" << LogIO::POST ;
5333 }
5334 return tsys ;
5335 }
5336}
5337
5338vector<float> STMath::getCalibratedSpectra( CountedPtr<Scantable>& on,
5339 CountedPtr<Scantable>& off,
5340 CountedPtr<Scantable>& sky,
5341 CountedPtr<Scantable>& hot,
5342 CountedPtr<Scantable>& cold,
5343 int index,
5344 string antname )
5345{
5346 (void) cold; //currently unused
5347 string reftime = on->getTime( index ) ;
5348 vector<int> ii( 1, on->getIF( index ) ) ;
5349 vector<int> ib( 1, on->getBeam( index ) ) ;
5350 vector<int> ip( 1, on->getPol( index ) ) ;
5351 STSelector sel = STSelector() ;
5352 sel.setIFs( ii ) ;
5353 sel.setBeams( ib ) ;
5354 sel.setPolarizations( ip ) ;
5355 sky->setSelection( sel ) ;
5356 hot->setSelection( sel ) ;
5357 //cold->setSelection( sel ) ;
5358 off->setSelection( sel ) ;
5359 vector<float> spsky = getSpectrumFromTime( reftime, sky, "linear" ) ;
5360 vector<float> sphot = getSpectrumFromTime( reftime, hot, "linear" ) ;
5361 //vector<float> spcold = getSpectrumFromTime( reftime, cold, "linear" ) ;
5362 vector<float> spoff = getSpectrumFromTime( reftime, off, "linear" ) ;
5363 vector<float> spec = on->getSpectrum( index ) ;
5364 vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ;
5365 vector<float> sp( tcal.size() ) ;
5366 if ( antname.find( "APEX" ) != string::npos ) {
5367 // using gain array
5368 for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
5369 float v = ( ( spec[j] - spoff[j] ) / spoff[j] )
5370 * ( spsky[j] / ( sphot[j] - spsky[j] ) ) * tcal[j] ;
5371 sp[j] = v ;
5372 }
5373 }
5374 else {
5375 // Chopper-Wheel calibration (Ulich & Haas 1976)
5376 for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
5377 float v = ( spec[j] - spoff[j] ) / ( sphot[j] - spsky[j] ) * tcal[j] ;
5378 sp[j] = v ;
5379 }
5380 }
5381 sel.reset() ;
5382 sky->unsetSelection() ;
5383 hot->unsetSelection() ;
5384 //cold->unsetSelection() ;
5385 off->unsetSelection() ;
5386
5387 return sp ;
5388}
5389
5390vector<float> STMath::getCalibratedSpectra( CountedPtr<Scantable>& on,
5391 CountedPtr<Scantable>& off,
5392 int index )
5393{
5394// string reftime = on->getTime( index ) ;
5395 ROScalarColumn<Double> timeCol( off->table(), "TIME" ) ;
5396 Vector<Double> timeVec = timeCol.getColumn() ;
5397 //ROTableColumn timeCol( on->table(), "TIME" ) ;
5398 timeCol.attach( on->table(), "TIME" ) ;
5399 double reftime = timeCol.asdouble(index) ;
5400 vector<int> ii( 1, on->getIF( index ) ) ;
5401 vector<int> ib( 1, on->getBeam( index ) ) ;
5402 vector<int> ip( 1, on->getPol( index ) ) ;
5403 STSelector sel = STSelector() ;
5404 sel.setIFs( ii ) ;
5405 sel.setBeams( ib ) ;
5406 sel.setPolarizations( ip ) ;
5407 off->setSelection( sel ) ;
5408 Vector<Float> spoff = getSpectrumFromTime( reftime, timeVec, off, "linear" ) ;
5409 vector<float> spec = on->getSpectrum( index ) ;
5410 //vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ;
5411 //vector<float> tsys = on->getTsysVec( index ) ;
5412 ArrayColumn<Float> tsysCol( on->table(), "TSYS" ) ;
5413 Vector<Float> tsys = tsysCol( index ) ;
5414 vector<float> sp( spec.size() ) ;
5415 // ALMA Calibration
5416 //
5417 // Ta* = Tsys * ( ON - OFF ) / OFF
5418 //
5419 // 2010/01/07 Takeshi Nakazato
5420 unsigned int tsyssize = tsys.nelements() ;
5421 unsigned int spsize = sp.size() ;
5422 for ( unsigned int j = 0 ; j < sp.size() ; j++ ) {
5423 float tscale = 0.0 ;
5424 if ( tsyssize == spsize )
5425 tscale = tsys[j] ;
5426 else
5427 tscale = tsys[0] ;
5428 float v = tscale * ( spec[j] - spoff[j] ) / spoff[j] ;
5429 sp[j] = v ;
5430 }
5431 sel.reset() ;
5432 off->unsetSelection() ;
5433
5434 return sp ;
5435}
5436
5437void STMath::calibrateALMA( CountedPtr<Scantable>& out,
5438 const CountedPtr<Scantable>& on,
5439 const CountedPtr<Scantable>& off,
5440 Vector<uInt>& rows )
5441{
5442 // 2012/05/22 TN
5443 // Assume that out has empty SPECTRA column
5444
5445 // if rows is empty, just return
5446 if ( rows.nelements() == 0 )
5447 return ;
5448 ROScalarColumn<Double> timeCol( off->table(), "TIME" ) ;
5449 Vector<Double> timeVec = timeCol.getColumn() ;
5450 timeCol.attach( on->table(), "TIME" ) ;
5451 ROArrayColumn<Float> tsysCol( on->table(), "TSYS" ) ;
5452 unsigned int spsize = on->nchan( on->getIF(rows[0]) ) ;
5453 // I know that the data is contiguous
5454 const uInt *p = rows.data() ;
5455 for ( int irow = 0 ; irow < rows.nelements() ; irow++ ) {
5456 double reftime = timeCol.asdouble(*p) ;
5457 Vector<Float> spoff = getSpectrumFromTime( reftime, timeVec, off, "linear" ) ;
5458 Vector<Float> spec = on->specCol_( *p ) ;
5459 Vector<Float> tsys = tsysCol( *p ) ;
5460 // ALMA Calibration
5461 //
5462 // Ta* = Tsys * ( ON - OFF ) / OFF
5463 //
5464 // 2010/01/07 Takeshi Nakazato
5465 unsigned int tsyssize = tsys.nelements() ;
5466 for ( unsigned int j = 0 ; j < spsize ; j++ ) {
5467 spec[j] = ( spec[j] - spoff[j] ) / spoff[j] ;
5468 if ( tsyssize == spsize )
5469 spec[j] *= tsys[j] ;
5470 else
5471 spec[j] *= tsys[0] ;
5472 }
5473 out->specCol_.put( *p, spec ) ;
5474 p++ ;
5475 }
5476}
5477
5478vector<float> STMath::getFSCalibratedSpectra( CountedPtr<Scantable>& sig,
5479 CountedPtr<Scantable>& ref,
5480 CountedPtr<Scantable>& sky,
5481 CountedPtr<Scantable>& hot,
5482 CountedPtr<Scantable>& cold,
5483 int index )
5484{
5485 (void) cold; //currently unused
5486 string reftime = sig->getTime( index ) ;
5487 vector<int> ii( 1, sig->getIF( index ) ) ;
5488 vector<int> ib( 1, sig->getBeam( index ) ) ;
5489 vector<int> ip( 1, sig->getPol( index ) ) ;
5490 vector<int> ic( 1, sig->getScan( index ) ) ;
5491 STSelector sel = STSelector() ;
5492 sel.setIFs( ii ) ;
5493 sel.setBeams( ib ) ;
5494 sel.setPolarizations( ip ) ;
5495 sky->setSelection( sel ) ;
5496 hot->setSelection( sel ) ;
5497 //cold->setSelection( sel ) ;
5498 vector<float> spsky = getSpectrumFromTime( reftime, sky, "linear" ) ;
5499 vector<float> sphot = getSpectrumFromTime( reftime, hot, "linear" ) ;
5500 //vector<float> spcold = getSpectrumFromTime( reftime, cold, "linear" ) ;
5501 vector<float> spref = ref->getSpectrum( index ) ;
5502 vector<float> spsig = sig->getSpectrum( index ) ;
5503 vector<float> tcal = getTcalFromTime( reftime, sky, "linear" ) ;
5504 vector<float> sp( tcal.size() ) ;
5505 for ( unsigned int j = 0 ; j < tcal.size() ; j++ ) {
5506 float v = tcal[j] * spsky[j] / ( sphot[j] - spsky[j] ) * ( spsig[j] - spref[j] ) / spref[j] ;
5507 sp[j] = v ;
5508 }
5509 sel.reset() ;
5510 sky->unsetSelection() ;
5511 hot->unsetSelection() ;
5512 //cold->unsetSelection() ;
5513
5514 return sp ;
5515}
5516
5517vector<float> STMath::getFSCalibratedSpectra( CountedPtr<Scantable>& sig,
5518 CountedPtr<Scantable>& ref,
5519 vector< CountedPtr<Scantable> >& sky,
5520 vector< CountedPtr<Scantable> >& hot,
5521 vector< CountedPtr<Scantable> >& cold,
5522 int index )
5523{
5524 (void) cold; //currently unused
5525 string reftime = sig->getTime( index ) ;
5526 vector<int> ii( 1, sig->getIF( index ) ) ;
5527 vector<int> ib( 1, sig->getBeam( index ) ) ;
5528 vector<int> ip( 1, sig->getPol( index ) ) ;
5529 vector<int> ic( 1, sig->getScan( index ) ) ;
5530 STSelector sel = STSelector() ;
5531 sel.setIFs( ii ) ;
5532 sel.setBeams( ib ) ;
5533 sel.setPolarizations( ip ) ;
5534 sky[0]->setSelection( sel ) ;
5535 hot[0]->setSelection( sel ) ;
5536 //cold[0]->setSelection( sel ) ;
5537 vector<float> spskys = getSpectrumFromTime( reftime, sky[0], "linear" ) ;
5538 vector<float> sphots = getSpectrumFromTime( reftime, hot[0], "linear" ) ;
5539 //vector<float> spcolds = getSpectrumFromTime( reftime, cold[0], "linear" ) ;
5540 vector<float> tcals = getTcalFromTime( reftime, sky[0], "linear" ) ;
5541 sel.reset() ;
5542 ii[0] = ref->getIF( index ) ;
5543 sel.setIFs( ii ) ;
5544 sel.setBeams( ib ) ;
5545 sel.setPolarizations( ip ) ;
5546 sky[1]->setSelection( sel ) ;
5547 hot[1]->setSelection( sel ) ;
5548 //cold[1]->setSelection( sel ) ;
5549 vector<float> spskyr = getSpectrumFromTime( reftime, sky[1], "linear" ) ;
5550 vector<float> sphotr = getSpectrumFromTime( reftime, hot[1], "linear" ) ;
5551 //vector<float> spcoldr = getSpectrumFromTime( reftime, cold[1], "linear" ) ;
5552 vector<float> tcalr = getTcalFromTime( reftime, sky[1], "linear" ) ;
5553 vector<float> spref = ref->getSpectrum( index ) ;
5554 vector<float> spsig = sig->getSpectrum( index ) ;
5555 vector<float> sp( tcals.size() ) ;
5556 for ( unsigned int j = 0 ; j < tcals.size() ; j++ ) {
5557 float v = tcals[j] * spsig[j] / ( sphots[j] - spskys[j] ) - tcalr[j] * spref[j] / ( sphotr[j] - spskyr[j] ) ;
5558 sp[j] = v ;
5559 }
5560 sel.reset() ;
5561 sky[0]->unsetSelection() ;
5562 hot[0]->unsetSelection() ;
5563 //cold[0]->unsetSelection() ;
5564 sky[1]->unsetSelection() ;
5565 hot[1]->unsetSelection() ;
5566 //cold[1]->unsetSelection() ;
5567
5568 return sp ;
5569}
5570
5571void STMath::copyRows( Table &out,
5572 const Table &in,
5573 uInt startout,
5574 uInt startin,
5575 uInt nrow,
5576 Bool copySpectra,
5577 Bool copyFlagtra,
5578 Bool copyTsys )
5579{
5580 uInt nexclude = 0 ;
5581 Block<String> excludeColsBlock( 3 ) ;
5582 if ( !copySpectra ) {
5583 excludeColsBlock[nexclude] = "SPECTRA" ;
5584 nexclude++ ;
5585 }
5586 if ( !copyFlagtra ) {
5587 excludeColsBlock[nexclude] = "FLAGTRA" ;
5588 nexclude++ ;
5589 }
5590 if ( !copyTsys ) {
5591 excludeColsBlock[nexclude] = "TSYS" ;
5592 nexclude++ ;
5593 }
5594 // if ( nexclude < 3 ) {
5595 // excludeCols.resize( nexclude, True ) ;
5596 // }
5597 Vector<String> excludeCols( IPosition(1,nexclude),
5598 excludeColsBlock.storage(),
5599 SHARE ) ;
5600// cout << "excludeCols=" << excludeCols << endl ;
5601 TableRow rowout( out, excludeCols, True ) ;
5602 ROTableRow rowin( in, excludeCols, True ) ;
5603 uInt rin = startin ;
5604 uInt rout = startout ;
5605 for ( uInt i = 0 ; i < nrow ; i++ ) {
5606 rowin.get( rin ) ;
5607 rowout.putMatchingFields( rout, rowin.record() ) ;
5608 rin++ ;
5609 rout++ ;
5610 }
5611}
Note: See TracBrowser for help on using the repository browser.