source: trunk/src/STMath.cpp@ 2804

Last change on this file since 2804 was 2695, checked in by Takeshi Nakazato, 12 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...

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