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

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

New Development: No

JIRA Issue: No

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...

Do not resize Vector everytime when new element is added.
Instead, we can allocate enough memory first.


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