source: branches/alma/src/STMath.cpp@ 1633

Last change on this file since 1633 was 1633, checked in by Takeshi Nakazato, 16 years ago

New Development: No

JIRA Issue: Yes CAS-1422

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: Defined calibrate() and almacal() in python/asapmath.py

Defined cwcal() in STMath class

Test Programs: List test programs

Put in Release Notes: Yes

Module(s): task_sdaverage

Description: Describe your changes here...

asapmath.py is changed to be able to calibrate data for APEX.
This modification is tentatively considered a calibration of ALMA
single-dish data. However, it must be updated in the future since
I don't know how raw ALMA data are provided and I have to change
code along with read ALMA data.

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