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

Last change on this file since 1614 was 1611, checked in by Takeshi Nakazato, 15 years ago

New Development: No

JIRA Issue: Yes CAS-1423

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Averaging methods are changed to refer DIRECTION all the time.
This change affects the behavior of sdaverage when timeaverage is
set to True.
Tolerance for direction is set to 1.0e-15 rad for OTF data and
is set to 1.0e-5 rad (around 2 arcsec) for non-OTF data.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 107.7 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 "MathUtils.h"
47#include "RowAccumulator.h"
48#include "STAttr.h"
49#include "STMath.h"
50#include "STSelector.h"
51
52using namespace casa;
53
54using namespace asap;
55
56// tolerance for direction comparison (rad)
57#define TOL_OTF 1.0e-15
58#define TOL_POINT 1.0e-5 // 2 arcsec
59
60STMath::STMath(bool insitu) :
61 insitu_(insitu)
62{
63}
64
65
66STMath::~STMath()
67{
68}
69
70CountedPtr<Scantable>
71STMath::average( const std::vector<CountedPtr<Scantable> >& in,
72 const std::vector<bool>& mask,
73 const std::string& weight,
74 const std::string& avmode)
75{
76 if ( avmode == "SCAN" && in.size() != 1 )
77 throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
78 "Use merge first."));
79 WeightType wtype = stringToWeight(weight);
80
81 // check if OTF observation
82 String obstype = in[0]->getHeader().obstype ;
83 //bool otfscan = false ;
84 Double tol = 0.0 ;
85 if ( obstype.find( "OTF" ) != String::npos ) {
86 //cout << "OTF scan" << endl ;
87 //otfscan = true ;
88 tol = TOL_OTF ;
89 }
90 else {
91 tol = TOL_POINT ;
92 }
93
94 // output
95 // clone as this is non insitu
96 bool insitu = insitu_;
97 setInsitu(false);
98 CountedPtr< Scantable > out = getScantable(in[0], true);
99 setInsitu(insitu);
100 std::vector<CountedPtr<Scantable> >::const_iterator stit = in.begin();
101 ++stit;
102 while ( stit != in.end() ) {
103 out->appendToHistoryTable((*stit)->history());
104 ++stit;
105 }
106
107 Table& tout = out->table();
108
109 /// @todo check if all scantables are conformant
110
111 ArrayColumn<Float> specColOut(tout,"SPECTRA");
112 ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
113 ArrayColumn<Float> tsysColOut(tout,"TSYS");
114 ScalarColumn<Double> mjdColOut(tout,"TIME");
115 ScalarColumn<Double> intColOut(tout,"INTERVAL");
116 ScalarColumn<uInt> cycColOut(tout,"CYCLENO");
117 ScalarColumn<uInt> scanColOut(tout,"SCANNO");
118
119 // set up the output table rows. These are based on the structure of the
120 // FIRST scantable in the vector
121 const Table& baset = in[0]->table();
122
123 Block<String> cols(3);
124 cols[0] = String("BEAMNO");
125 cols[1] = String("IFNO");
126 cols[2] = String("POLNO");
127 if ( avmode == "SOURCE" ) {
128 cols.resize(4);
129 cols[3] = String("SRCNAME");
130 }
131 if ( avmode == "SCAN" && in.size() == 1) {
132 //cols.resize(4);
133 //cols[3] = String("SCANNO");
134 cols.resize(5);
135 cols[3] = String("SRCNAME");
136 cols[4] = String("SCANNO");
137 }
138 uInt outrowCount = 0;
139 TableIterator iter(baset, cols);
140// int count = 0 ;
141 while (!iter.pastEnd()) {
142 Table subt = iter.table();
143// // copy the first row of this selection into the new table
144// tout.addRow();
145// TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
146// // re-index to 0
147// if ( avmode != "SCAN" && avmode != "SOURCE" ) {
148// scanColOut.put(outrowCount, uInt(0));
149// }
150// ++outrowCount;
151 MDirection::ScalarColumn dircol ;
152 dircol.attach( subt, "DIRECTION" ) ;
153 Int length = subt.nrow() ;
154 vector< Vector<Double> > dirs ;
155 vector<int> indexes ;
156 for ( Int i = 0 ; i < length ; i++ ) {
157 Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
158 //cout << setw( 3 ) << count++ << ": " ;
159 //cout << "[" << t[0] << "," << t[1] << "]" << endl ;
160 bool adddir = true ;
161 for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
162 //if ( allTrue( t == dirs[j] ) ) {
163 Double dx = t[0] - dirs[j][0] ;
164 Double dy = t[1] - dirs[j][1] ;
165 Double dd = sqrt( dx * dx + dy * dy ) ;
166 //if ( allNearAbs( t, dirs[j], tol ) ) {
167 if ( dd <= tol ) {
168 adddir = false ;
169 break ;
170 }
171 }
172 if ( adddir ) {
173 dirs.push_back( t ) ;
174 indexes.push_back( i ) ;
175 }
176 }
177 uInt rowNum = dirs.size() ;
178 //cout << "dirs.size()=" << dirs.size() << endl ;
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 //cout << "removeRows.size()=" << removeRows.size() << endl ;
236 if ( removeRows.size() != 0 ) {
237 //cout << "[" ;
238 //for ( uInt irow=0 ; irow<removeRows.size()-1 ; irow++ )
239 //cout << removeRows[irow] << "," ;
240 //cout << removeRows[removeRows.size()-1] << "]" << endl ;
241 subt.removeRow( removeRows ) ;
242 }
243
244 if ( nrsubt == removeRows.size() )
245 throw(AipsError("Averaging data is empty.")) ;
246
247 specCol.attach(subt,"SPECTRA");
248 flagCol.attach(subt,"FLAGTRA");
249 tsysCol.attach(subt,"TSYS");
250 intCol.attach(subt,"INTERVAL");
251 mjdCol.attach(subt,"TIME");
252 Vector<Float> spec,tsys;
253 Vector<uChar> flag;
254 Double inter,time;
255 for (uInt k = 0; k < subt.nrow(); ++k ) {
256 flagCol.get(k, flag);
257 Vector<Bool> bflag(flag.shape());
258 convertArray(bflag, flag);
259 /*
260 if ( allEQ(bflag, True) ) {
261 continue;//don't accumulate
262 }
263 */
264 specCol.get(k, spec);
265 tsysCol.get(k, tsys);
266 intCol.get(k, inter);
267 mjdCol.get(k, time);
268 // spectrum has to be added last to enable weighting by the other values
269 acc.add(spec, !bflag, tsys, inter, time);
270 }
271 }
272 const Vector<Bool>& msk = acc.getMask();
273 if ( allEQ(msk, False) ) {
274 uint n = rowstodelete.nelements();
275 rowstodelete.resize(n+1, True);
276 rowstodelete[n] = i;
277 continue;
278 }
279 //write out
280 if (acc.state()) {
281 Vector<uChar> flg(msk.shape());
282 convertArray(flg, !msk);
283 flagColOut.put(i, flg);
284 specColOut.put(i, acc.getSpectrum());
285 tsysColOut.put(i, acc.getTsys());
286 intColOut.put(i, acc.getInterval());
287 mjdColOut.put(i, acc.getTime());
288 // we should only have one cycle now -> reset it to be 0
289 // frequency switched data has different CYCLENO for different IFNO
290 // which requires resetting this value
291 cycColOut.put(i, uInt(0));
292 } else {
293 ostringstream oss;
294 oss << "For output row="<<i<<", all input rows of data are flagged. no averaging" << endl;
295 pushLog(String(oss));
296 }
297 acc.reset();
298 }
299 if (rowstodelete.nelements() > 0) {
300 cout << rowstodelete << endl;
301 tout.removeRow(rowstodelete);
302 if (tout.nrow() == 0) {
303 throw(AipsError("Can't average fully flagged data."));
304 }
305 }
306 return out;
307}
308
309CountedPtr< Scantable >
310 STMath::averageChannel( const CountedPtr < Scantable > & in,
311 const std::string & mode,
312 const std::string& avmode )
313{
314 // check if OTF observation
315 String obstype = in->getHeader().obstype ;
316 //bool otfscan = false ;
317 Double tol = 0.0 ;
318 if ( obstype.find( "OTF" ) != String::npos ) {
319 //cout << "OTF scan" << endl ;
320 //otfscan = true ;
321 tol = TOL_OTF ;
322 }
323 else {
324 tol = TOL_POINT ;
325 }
326
327 // clone as this is non insitu
328 bool insitu = insitu_;
329 setInsitu(false);
330 CountedPtr< Scantable > out = getScantable(in, true);
331 setInsitu(insitu);
332 Table& tout = out->table();
333 ArrayColumn<Float> specColOut(tout,"SPECTRA");
334 ArrayColumn<uChar> flagColOut(tout,"FLAGTRA");
335 ArrayColumn<Float> tsysColOut(tout,"TSYS");
336 ScalarColumn<uInt> scanColOut(tout,"SCANNO");
337 ScalarColumn<Double> intColOut(tout, "INTERVAL");
338 Table tmp = in->table().sort("BEAMNO");
339 Block<String> cols(3);
340 cols[0] = String("BEAMNO");
341 cols[1] = String("IFNO");
342 cols[2] = String("POLNO");
343 if ( avmode == "SCAN") {
344 cols.resize(4);
345 cols[3] = String("SCANNO");
346 }
347 uInt outrowCount = 0;
348 uChar userflag = 1 << 7;
349 TableIterator iter(tmp, cols);
350 while (!iter.pastEnd()) {
351 Table subt = iter.table();
352 ROArrayColumn<Float> specCol, tsysCol;
353 ROArrayColumn<uChar> flagCol;
354 ROScalarColumn<Double> intCol(subt, "INTERVAL");
355 specCol.attach(subt,"SPECTRA");
356 flagCol.attach(subt,"FLAGTRA");
357 tsysCol.attach(subt,"TSYS");
358// tout.addRow();
359// TableCopy::copyRows(tout, subt, outrowCount, 0, 1);
360// if ( avmode != "SCAN") {
361// scanColOut.put(outrowCount, uInt(0));
362// }
363// Vector<Float> tmp;
364// specCol.get(0, tmp);
365// uInt nchan = tmp.nelements();
366// // have to do channel by channel here as MaskedArrMath
367// // doesn't have partialMedians
368// Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
369// Vector<Float> outspec(nchan);
370// Vector<uChar> outflag(nchan,0);
371// Vector<Float> outtsys(1);/// @fixme when tsys is channel based
372// for (uInt i=0; i<nchan; ++i) {
373// Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
374// MaskedArray<Float> ma = maskedArray(specs,flags);
375// outspec[i] = median(ma);
376// if ( allEQ(ma.getMask(), False) )
377// outflag[i] = userflag;// flag data
378// }
379// outtsys[0] = median(tsysCol.getColumn());
380// specColOut.put(outrowCount, outspec);
381// flagColOut.put(outrowCount, outflag);
382// tsysColOut.put(outrowCount, outtsys);
383// Double intsum = sum(intCol.getColumn());
384// intColOut.put(outrowCount, intsum);
385// ++outrowCount;
386// ++iter;
387 MDirection::ScalarColumn dircol ;
388 dircol.attach( subt, "DIRECTION" ) ;
389 Int length = subt.nrow() ;
390 vector< Vector<Double> > dirs ;
391 vector<int> indexes ;
392 for ( Int i = 0 ; i < length ; i++ ) {
393 Vector<Double> t = dircol(i).getAngle(Unit(String("rad"))).getValue() ;
394 bool adddir = true ;
395 for ( uInt j = 0 ; j < dirs.size() ; j++ ) {
396 //if ( allTrue( t == dirs[j] ) ) {
397 Double dx = t[0] - dirs[j][0] ;
398 Double dy = t[1] - dirs[j][1] ;
399 Double dd = sqrt( dx * dx + dy * dy ) ;
400 //if ( allNearAbs( t, dirs[j], tol ) ) {
401 if ( dd <= tol ) {
402 adddir = false ;
403 break ;
404 }
405 }
406 if ( adddir ) {
407 dirs.push_back( t ) ;
408 indexes.push_back( i ) ;
409 }
410 }
411 uInt rowNum = dirs.size() ;
412 tout.addRow( rowNum );
413 for ( uInt i = 0 ; i < rowNum ; i++ ) {
414 TableCopy::copyRows(tout, subt, outrowCount+i, indexes[i], 1) ;
415 if ( avmode != "SCAN") {
416 //scanColOut.put(outrowCount+i, uInt(0));
417 }
418 }
419 MDirection::ScalarColumn dircolOut ;
420 dircolOut.attach( tout, "DIRECTION" ) ;
421 for ( uInt irow = 0 ; irow < rowNum ; irow++ ) {
422 Vector<Double> t = dircolOut(outrowCount+irow).getAngle(Unit(String("rad"))).getValue() ;
423 Vector<Float> tmp;
424 specCol.get(0, tmp);
425 uInt nchan = tmp.nelements();
426 // have to do channel by channel here as MaskedArrMath
427 // doesn't have partialMedians
428 Vector<uChar> flags = flagCol.getColumn(Slicer(Slice(0)));
429 // mask spectra for different DIRECTION
430 //cout << "irow=" << outrowCount+irow << ": flagged [" ;
431 for ( uInt jrow = 0 ; jrow < subt.nrow() ; jrow++ ) {
432 Vector<Double> direction = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
433 //if ( t[0] != direction[0] || t[1] != direction[1] ) {
434 Double dx = t[0] - direction[0] ;
435 Double dy = t[1] - direction[1] ;
436 Double dd = sqrt( dx * dx + dy * dy ) ;
437 //if ( !allNearAbs( t, direction, tol ) ) {
438 if ( dd > tol ) {
439 //cout << jrow << " " ;
440 flags[jrow] = userflag ;
441 }
442 }
443 //cout << "]" << endl ;
444 //cout << "flags=" << flags << endl ;
445 Vector<Float> outspec(nchan);
446 Vector<uChar> outflag(nchan,0);
447 Vector<Float> outtsys(1);/// @fixme when tsys is channel based
448 for (uInt i=0; i<nchan; ++i) {
449 Vector<Float> specs = specCol.getColumn(Slicer(Slice(i)));
450 MaskedArray<Float> ma = maskedArray(specs,flags);
451 outspec[i] = median(ma);
452 if ( allEQ(ma.getMask(), False) )
453 outflag[i] = userflag;// flag data
454 }
455 outtsys[0] = median(tsysCol.getColumn());
456 specColOut.put(outrowCount+irow, outspec);
457 flagColOut.put(outrowCount+irow, outflag);
458 tsysColOut.put(outrowCount+irow, outtsys);
459 Vector<Double> integ = intCol.getColumn() ;
460 MaskedArray<Double> mi = maskedArray( integ, flags ) ;
461 Double intsum = sum(mi);
462 intColOut.put(outrowCount+irow, intsum);
463 }
464 outrowCount += rowNum ;
465 ++iter;
466 }
467 return out;
468}
469
470CountedPtr< Scantable > STMath::getScantable(const CountedPtr< Scantable >& in,
471 bool droprows)
472{
473 if (insitu_) {
474 return in;
475 }
476 else {
477 // clone
478 return CountedPtr<Scantable>(new Scantable(*in, Bool(droprows)));
479 }
480}
481
482CountedPtr< Scantable > STMath::unaryOperate( const CountedPtr< Scantable >& in,
483 float val,
484 const std::string& mode,
485 bool tsys )
486{
487 CountedPtr< Scantable > out = getScantable(in, false);
488 Table& tab = out->table();
489 ArrayColumn<Float> specCol(tab,"SPECTRA");
490 ArrayColumn<Float> tsysCol(tab,"TSYS");
491 for (uInt i=0; i<tab.nrow(); ++i) {
492 Vector<Float> spec;
493 Vector<Float> ts;
494 specCol.get(i, spec);
495 tsysCol.get(i, ts);
496 if (mode == "MUL" || mode == "DIV") {
497 if (mode == "DIV") val = 1.0/val;
498 spec *= val;
499 specCol.put(i, spec);
500 if ( tsys ) {
501 ts *= val;
502 tsysCol.put(i, ts);
503 }
504 } else if ( mode == "ADD" || mode == "SUB") {
505 if (mode == "SUB") val *= -1.0;
506 spec += val;
507 specCol.put(i, spec);
508 if ( tsys ) {
509 ts += val;
510 tsysCol.put(i, ts);
511 }
512 }
513 }
514 return out;
515}
516
517CountedPtr<Scantable> STMath::binaryOperate(const CountedPtr<Scantable>& left,
518 const CountedPtr<Scantable>& right,
519 const std::string& mode)
520{
521 bool insitu = insitu_;
522 if ( ! left->conformant(*right) ) {
523 throw(AipsError("'left' and 'right' scantables are not conformant."));
524 }
525 setInsitu(false);
526 CountedPtr< Scantable > out = getScantable(left, false);
527 setInsitu(insitu);
528 Table& tout = out->table();
529 Block<String> coln(5);
530 coln[0] = "SCANNO"; coln[1] = "CYCLENO"; coln[2] = "BEAMNO";
531 coln[3] = "IFNO"; coln[4] = "POLNO";
532 Table tmpl = tout.sort(coln);
533 Table tmpr = right->table().sort(coln);
534 ArrayColumn<Float> lspecCol(tmpl,"SPECTRA");
535 ROArrayColumn<Float> rspecCol(tmpr,"SPECTRA");
536 ArrayColumn<uChar> lflagCol(tmpl,"FLAGTRA");
537 ROArrayColumn<uChar> rflagCol(tmpr,"FLAGTRA");
538
539 for (uInt i=0; i<tout.nrow(); ++i) {
540 Vector<Float> lspecvec, rspecvec;
541 Vector<uChar> lflagvec, rflagvec;
542 lspecvec = lspecCol(i); rspecvec = rspecCol(i);
543 lflagvec = lflagCol(i); rflagvec = rflagCol(i);
544 MaskedArray<Float> mleft = maskedArray(lspecvec, lflagvec);
545 MaskedArray<Float> mright = maskedArray(rspecvec, rflagvec);
546 if (mode == "ADD") {
547 mleft += mright;
548 } else if ( mode == "SUB") {
549 mleft -= mright;
550 } else if ( mode == "MUL") {
551 mleft *= mright;
552 } else if ( mode == "DIV") {
553 mleft /= mright;
554 } else {
555 throw(AipsError("Illegal binary operator"));
556 }
557 lspecCol.put(i, mleft.getArray());
558 }
559 return out;
560}
561
562
563
564MaskedArray<Float> STMath::maskedArray( const Vector<Float>& s,
565 const Vector<uChar>& f)
566{
567 Vector<Bool> mask;
568 mask.resize(f.shape());
569 convertArray(mask, f);
570 return MaskedArray<Float>(s,!mask);
571}
572
573MaskedArray<Double> STMath::maskedArray( const Vector<Double>& s,
574 const Vector<uChar>& f)
575{
576 Vector<Bool> mask;
577 mask.resize(f.shape());
578 convertArray(mask, f);
579 return MaskedArray<Double>(s,!mask);
580}
581
582Vector<uChar> STMath::flagsFromMA(const MaskedArray<Float>& ma)
583{
584 const Vector<Bool>& m = ma.getMask();
585 Vector<uChar> flags(m.shape());
586 convertArray(flags, !m);
587 return flags;
588}
589
590CountedPtr< Scantable > STMath::autoQuotient( const CountedPtr< Scantable >& in,
591 const std::string & mode,
592 bool preserve )
593{
594 /// @todo make other modes available
595 /// modes should be "nearest", "pair"
596 // make this operation non insitu
597 const Table& tin = in->table();
598 Table ons = tin(tin.col("SRCTYPE") == Int(0));
599 Table offs = tin(tin.col("SRCTYPE") == Int(1));
600 if ( offs.nrow() == 0 )
601 throw(AipsError("No 'off' scans present."));
602 // put all "on" scans into output table
603
604 bool insitu = insitu_;
605 setInsitu(false);
606 CountedPtr< Scantable > out = getScantable(in, true);
607 setInsitu(insitu);
608 Table& tout = out->table();
609
610 TableCopy::copyRows(tout, ons);
611 TableRow row(tout);
612 ROScalarColumn<Double> offtimeCol(offs, "TIME");
613 ArrayColumn<Float> outspecCol(tout, "SPECTRA");
614 ROArrayColumn<Float> outtsysCol(tout, "TSYS");
615 ArrayColumn<uChar> outflagCol(tout, "FLAGTRA");
616 for (uInt i=0; i < tout.nrow(); ++i) {
617 const TableRecord& rec = row.get(i);
618 Double ontime = rec.asDouble("TIME");
619 Table presel = offs(offs.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
620 && offs.col("IFNO") == Int(rec.asuInt("IFNO"))
621 && offs.col("POLNO") == Int(rec.asuInt("POLNO")) );
622 ROScalarColumn<Double> offtimeCol(presel, "TIME");
623
624 Double mindeltat = min(abs(offtimeCol.getColumn() - ontime));
625 // Timestamp may vary within a cycle ???!!!
626 // increase this by 0.01 sec in case of rounding errors...
627 // There might be a better way to do this.
628 // fix to this fix. TIME is MJD, so 1.0d not 1.0s
629 mindeltat += 0.01/24./60./60.;
630 Table sel = presel( abs(presel.col("TIME")-ontime) <= mindeltat);
631
632 if ( sel.nrow() < 1 ) {
633 throw(AipsError("No closest in time found... This could be a rounding "
634 "issue. Try quotient instead."));
635 }
636 TableRow offrow(sel);
637 const TableRecord& offrec = offrow.get(0);//should only be one row
638 RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
639 RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
640 RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
641 /// @fixme this assumes tsys is a scalar not vector
642 Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
643 Vector<Float> specon, tsyson;
644 outtsysCol.get(i, tsyson);
645 outspecCol.get(i, specon);
646 Vector<uChar> flagon;
647 outflagCol.get(i, flagon);
648 MaskedArray<Float> mon = maskedArray(specon, flagon);
649 MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
650 MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
651 if (preserve) {
652 quot -= tsysoffscalar;
653 } else {
654 quot -= tsyson[0];
655 }
656 outspecCol.put(i, quot.getArray());
657 outflagCol.put(i, flagsFromMA(quot));
658 }
659 // renumber scanno
660 TableIterator it(tout, "SCANNO");
661 uInt i = 0;
662 while ( !it.pastEnd() ) {
663 Table t = it.table();
664 TableVector<uInt> vec(t, "SCANNO");
665 vec = i;
666 ++i;
667 ++it;
668 }
669 return out;
670}
671
672
673CountedPtr< Scantable > STMath::quotient( const CountedPtr< Scantable > & on,
674 const CountedPtr< Scantable > & off,
675 bool preserve )
676{
677 bool insitu = insitu_;
678 if ( ! on->conformant(*off) ) {
679 throw(AipsError("'on' and 'off' scantables are not conformant."));
680 }
681 setInsitu(false);
682 CountedPtr< Scantable > out = getScantable(on, false);
683 setInsitu(insitu);
684 Table& tout = out->table();
685 const Table& toff = off->table();
686 TableIterator sit(tout, "SCANNO");
687 TableIterator s2it(toff, "SCANNO");
688 while ( !sit.pastEnd() ) {
689 Table ton = sit.table();
690 TableRow row(ton);
691 Table t = s2it.table();
692 ArrayColumn<Float> outspecCol(ton, "SPECTRA");
693 ROArrayColumn<Float> outtsysCol(ton, "TSYS");
694 ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
695 for (uInt i=0; i < ton.nrow(); ++i) {
696 const TableRecord& rec = row.get(i);
697 Table offsel = t( t.col("BEAMNO") == Int(rec.asuInt("BEAMNO"))
698 && t.col("IFNO") == Int(rec.asuInt("IFNO"))
699 && t.col("POLNO") == Int(rec.asuInt("POLNO")) );
700 if ( offsel.nrow() == 0 )
701 throw AipsError("STMath::quotient: no matching off");
702 TableRow offrow(offsel);
703 const TableRecord& offrec = offrow.get(0);//should be ncycles - take first
704 RORecordFieldPtr< Array<Float> > specoff(offrec, "SPECTRA");
705 RORecordFieldPtr< Array<Float> > tsysoff(offrec, "TSYS");
706 RORecordFieldPtr< Array<uChar> > flagoff(offrec, "FLAGTRA");
707 Float tsysoffscalar = (*tsysoff)(IPosition(1,0));
708 Vector<Float> specon, tsyson;
709 outtsysCol.get(i, tsyson);
710 outspecCol.get(i, specon);
711 Vector<uChar> flagon;
712 outflagCol.get(i, flagon);
713 MaskedArray<Float> mon = maskedArray(specon, flagon);
714 MaskedArray<Float> moff = maskedArray(*specoff, *flagoff);
715 MaskedArray<Float> quot = (tsysoffscalar * mon / moff);
716 if (preserve) {
717 quot -= tsysoffscalar;
718 } else {
719 quot -= tsyson[0];
720 }
721 outspecCol.put(i, quot.getArray());
722 outflagCol.put(i, flagsFromMA(quot));
723 }
724 ++sit;
725 ++s2it;
726 // take the first off for each on scan which doesn't have a
727 // matching off scan
728 // non <= noff: matching pairs, non > noff matching pairs then first off
729 if ( s2it.pastEnd() ) s2it.reset();
730 }
731 return out;
732}
733
734// dototalpower (migration of GBTIDL procedure dototalpower.pro)
735// calibrate the CAL on-off pair. It calculate Tsys and average CAL on-off subintegrations
736// do it for each cycles in a specific scan.
737CountedPtr< Scantable > STMath::dototalpower( const CountedPtr< Scantable >& calon,
738 const CountedPtr< Scantable >& caloff, Float tcal )
739{
740if ( ! calon->conformant(*caloff) ) {
741 throw(AipsError("'CAL on' and 'CAL off' scantables are not conformant."));
742 }
743 setInsitu(false);
744 CountedPtr< Scantable > out = getScantable(caloff, false);
745 Table& tout = out->table();
746 const Table& tcon = calon->table();
747 Vector<Float> tcalout;
748 Vector<Float> tcalout2; //debug
749
750 if ( tout.nrow() != tcon.nrow() ) {
751 throw(AipsError("Mismatch in number of rows to form cal on - off pair."));
752 }
753 // iteration by scanno or cycle no.
754 TableIterator sit(tout, "SCANNO");
755 TableIterator s2it(tcon, "SCANNO");
756 while ( !sit.pastEnd() ) {
757 Table toff = sit.table();
758 TableRow row(toff);
759 Table t = s2it.table();
760 ScalarColumn<Double> outintCol(toff, "INTERVAL");
761 ArrayColumn<Float> outspecCol(toff, "SPECTRA");
762 ArrayColumn<Float> outtsysCol(toff, "TSYS");
763 ArrayColumn<uChar> outflagCol(toff, "FLAGTRA");
764 ROScalarColumn<uInt> outtcalIdCol(toff, "TCAL_ID");
765 ROScalarColumn<uInt> outpolCol(toff, "POLNO");
766 ROScalarColumn<Double> onintCol(t, "INTERVAL");
767 ROArrayColumn<Float> onspecCol(t, "SPECTRA");
768 ROArrayColumn<Float> ontsysCol(t, "TSYS");
769 ROArrayColumn<uChar> onflagCol(t, "FLAGTRA");
770 //ROScalarColumn<uInt> ontcalIdCol(t, "TCAL_ID");
771
772 for (uInt i=0; i < toff.nrow(); ++i) {
773 //skip these checks -> assumes the data order are the same between the cal on off pairs
774 //
775 Vector<Float> specCalon, specCaloff;
776 // to store scalar (mean) tsys
777 Vector<Float> tsysout(1);
778 uInt tcalId, polno;
779 Double offint, onint;
780 outpolCol.get(i, polno);
781 outspecCol.get(i, specCaloff);
782 onspecCol.get(i, specCalon);
783 Vector<uChar> flagCaloff, flagCalon;
784 outflagCol.get(i, flagCaloff);
785 onflagCol.get(i, flagCalon);
786 outtcalIdCol.get(i, tcalId);
787 outintCol.get(i, offint);
788 onintCol.get(i, onint);
789 // caluculate mean Tsys
790 uInt nchan = specCaloff.nelements();
791 // percentage of edge cut off
792 uInt pc = 10;
793 uInt bchan = nchan/pc;
794 uInt echan = nchan-bchan;
795
796 Slicer chansl(IPosition(1,bchan-1), IPosition(1,echan-1), IPosition(1,1),Slicer::endIsLast);
797 Vector<Float> testsubsp = specCaloff(chansl);
798 MaskedArray<Float> spoff = maskedArray( specCaloff(chansl),flagCaloff(chansl) );
799 MaskedArray<Float> spon = maskedArray( specCalon(chansl),flagCalon(chansl) );
800 MaskedArray<Float> spdiff = spon-spoff;
801 uInt noff = spoff.nelementsValid();
802 //uInt non = spon.nelementsValid();
803 uInt ndiff = spdiff.nelementsValid();
804 Float meantsys;
805
806/**
807 Double subspec, subdiff;
808 uInt usednchan;
809 subspec = 0;
810 subdiff = 0;
811 usednchan = 0;
812 for(uInt k=(bchan-1); k<echan; k++) {
813 subspec += specCaloff[k];
814 subdiff += static_cast<Double>(specCalon[k]-specCaloff[k]);
815 ++usednchan;
816 }
817**/
818 // get tcal if input tcal <= 0
819 String tcalt;
820 Float tcalUsed;
821 tcalUsed = tcal;
822 if ( tcal <= 0.0 ) {
823 caloff->tcal().getEntry(tcalt, tcalout, tcalId);
824 if (polno<=3) {
825 tcalUsed = tcalout[polno];
826 }
827 else {
828 tcalUsed = tcalout[0];
829 }
830 }
831
832 Float meanoff;
833 Float meandiff;
834 if (noff && ndiff) {
835 //Debug
836 //if(noff!=ndiff) cerr<<"noff and ndiff is not equal"<<endl;
837 meanoff = sum(spoff)/noff;
838 meandiff = sum(spdiff)/ndiff;
839 meantsys= (meanoff/meandiff )*tcalUsed + tcalUsed/2;
840 }
841 else {
842 meantsys=1;
843 }
844
845 tsysout[0] = Float(meantsys);
846 MaskedArray<Float> mcaloff = maskedArray(specCaloff, flagCaloff);
847 MaskedArray<Float> mcalon = maskedArray(specCalon, flagCalon);
848 MaskedArray<Float> sig = Float(0.5) * (mcaloff + mcalon);
849 //uInt ncaloff = mcaloff.nelementsValid();
850 //uInt ncalon = mcalon.nelementsValid();
851
852 outintCol.put(i, offint+onint);
853 outspecCol.put(i, sig.getArray());
854 outflagCol.put(i, flagsFromMA(sig));
855 outtsysCol.put(i, tsysout);
856 }
857 ++sit;
858 ++s2it;
859 }
860 return out;
861}
862
863//dosigref - migrated from GBT IDL's dosigref.pro, do calibration of position switch
864// observatiions.
865// input: sig and ref scantables, and an optional boxcar smoothing width(default width=0,
866// no smoothing).
867// output: resultant scantable [= (sig-ref/ref)*tsys]
868CountedPtr< Scantable > STMath::dosigref( const CountedPtr < Scantable >& sig,
869 const CountedPtr < Scantable >& ref,
870 int smoothref,
871 casa::Float tsysv,
872 casa::Float tau )
873{
874if ( ! ref->conformant(*sig) ) {
875 throw(AipsError("'sig' and 'ref' scantables are not conformant."));
876 }
877 setInsitu(false);
878 CountedPtr< Scantable > out = getScantable(sig, false);
879 CountedPtr< Scantable > smref;
880 if ( smoothref > 1 ) {
881 float fsmoothref = static_cast<float>(smoothref);
882 std::string inkernel = "boxcar";
883 smref = smooth(ref, inkernel, fsmoothref );
884 ostringstream oss;
885 oss<<"Applied smoothing of "<<fsmoothref<<" on the reference."<<endl;
886 pushLog(String(oss));
887 }
888 else {
889 smref = ref;
890 }
891 Table& tout = out->table();
892 const Table& tref = smref->table();
893 if ( tout.nrow() != tref.nrow() ) {
894 throw(AipsError("Mismatch in number of rows to form on-source and reference pair."));
895 }
896 // iteration by scanno? or cycle no.
897 TableIterator sit(tout, "SCANNO");
898 TableIterator s2it(tref, "SCANNO");
899 while ( !sit.pastEnd() ) {
900 Table ton = sit.table();
901 Table t = s2it.table();
902 ScalarColumn<Double> outintCol(ton, "INTERVAL");
903 ArrayColumn<Float> outspecCol(ton, "SPECTRA");
904 ArrayColumn<Float> outtsysCol(ton, "TSYS");
905 ArrayColumn<uChar> outflagCol(ton, "FLAGTRA");
906 ArrayColumn<Float> refspecCol(t, "SPECTRA");
907 ROScalarColumn<Double> refintCol(t, "INTERVAL");
908 ROArrayColumn<Float> reftsysCol(t, "TSYS");
909 ArrayColumn<uChar> refflagCol(t, "FLAGTRA");
910 ROScalarColumn<Float> refelevCol(t, "ELEVATION");
911 for (uInt i=0; i < ton.nrow(); ++i) {
912
913 Double onint, refint;
914 Vector<Float> specon, specref;
915 // to store scalar (mean) tsys
916 Vector<Float> tsysref;
917 outintCol.get(i, onint);
918 refintCol.get(i, refint);
919 outspecCol.get(i, specon);
920 refspecCol.get(i, specref);
921 Vector<uChar> flagref, flagon;
922 outflagCol.get(i, flagon);
923 refflagCol.get(i, flagref);
924 reftsysCol.get(i, tsysref);
925
926 Float tsysrefscalar;
927 if ( tsysv > 0.0 ) {
928 ostringstream oss;
929 Float elev;
930 refelevCol.get(i, elev);
931 oss << "user specified Tsys = " << tsysv;
932 // do recalc elevation if EL = 0
933 if ( elev == 0 ) {
934 throw(AipsError("EL=0, elevation data is missing."));
935 } else {
936 if ( tau <= 0.0 ) {
937 throw(AipsError("Valid tau is not supplied."));
938 } else {
939 tsysrefscalar = tsysv * exp(tau/elev);
940 }
941 }
942 oss << ", corrected (for El) tsys= "<<tsysrefscalar;
943 pushLog(String(oss));
944 }
945 else {
946 tsysrefscalar = tsysref[0];
947 }
948 //get quotient spectrum
949 MaskedArray<Float> mref = maskedArray(specref, flagref);
950 MaskedArray<Float> mon = maskedArray(specon, flagon);
951 MaskedArray<Float> specres = tsysrefscalar*((mon - mref)/mref);
952 Double resint = onint*refint*smoothref/(onint+refint*smoothref);
953
954 //Debug
955 //cerr<<"Tsys used="<<tsysrefscalar<<endl;
956 // fill the result, replay signal tsys by reference tsys
957 outintCol.put(i, resint);
958 outspecCol.put(i, specres.getArray());
959 outflagCol.put(i, flagsFromMA(specres));
960 outtsysCol.put(i, tsysref);
961 }
962 ++sit;
963 ++s2it;
964 }
965 return out;
966}
967
968CountedPtr< Scantable > STMath::donod(const casa::CountedPtr<Scantable>& s,
969 const std::vector<int>& scans,
970 int smoothref,
971 casa::Float tsysv,
972 casa::Float tau,
973 casa::Float tcal )
974
975{
976 setInsitu(false);
977 STSelector sel;
978 std::vector<int> scan1, scan2, beams;
979 std::vector< vector<int> > scanpair;
980 std::vector<string> calstate;
981 String msg;
982
983 CountedPtr< Scantable > s1b1on, s1b1off, s1b2on, s1b2off;
984 CountedPtr< Scantable > s2b1on, s2b1off, s2b2on, s2b2off;
985
986 std::vector< CountedPtr< Scantable > > sctables;
987 sctables.push_back(s1b1on);
988 sctables.push_back(s1b1off);
989 sctables.push_back(s1b2on);
990 sctables.push_back(s1b2off);
991 sctables.push_back(s2b1on);
992 sctables.push_back(s2b1off);
993 sctables.push_back(s2b2on);
994 sctables.push_back(s2b2off);
995
996 //check scanlist
997 int n=s->checkScanInfo(scans);
998 if (n==1) {
999 throw(AipsError("Incorrect scan pairs. "));
1000 }
1001
1002 // Assume scans contain only a pair of consecutive scan numbers.
1003 // It is assumed that first beam, b1, is on target.
1004 // There is no check if the first beam is on or not.
1005 if ( scans.size()==1 ) {
1006 scan1.push_back(scans[0]);
1007 scan2.push_back(scans[0]+1);
1008 } else if ( scans.size()==2 ) {
1009 scan1.push_back(scans[0]);
1010 scan2.push_back(scans[1]);
1011 } else {
1012 if ( scans.size()%2 == 0 ) {
1013 for (uInt i=0; i<scans.size(); i++) {
1014 if (i%2 == 0) {
1015 scan1.push_back(scans[i]);
1016 }
1017 else {
1018 scan2.push_back(scans[i]);
1019 }
1020 }
1021 } else {
1022 throw(AipsError("Odd numbers of scans, cannot form pairs."));
1023 }
1024 }
1025 scanpair.push_back(scan1);
1026 scanpair.push_back(scan2);
1027 calstate.push_back("*calon");
1028 calstate.push_back("*[^calon]");
1029 CountedPtr< Scantable > ws = getScantable(s, false);
1030 uInt l=0;
1031 while ( l < sctables.size() ) {
1032 for (uInt i=0; i < 2; i++) {
1033 for (uInt j=0; j < 2; j++) {
1034 for (uInt k=0; k < 2; k++) {
1035 sel.reset();
1036 sel.setScans(scanpair[i]);
1037 sel.setName(calstate[k]);
1038 beams.clear();
1039 beams.push_back(j);
1040 sel.setBeams(beams);
1041 ws->setSelection(sel);
1042 sctables[l]= getScantable(ws, false);
1043 l++;
1044 }
1045 }
1046 }
1047 }
1048
1049 // replace here by splitData or getData functionality
1050 CountedPtr< Scantable > sig1;
1051 CountedPtr< Scantable > ref1;
1052 CountedPtr< Scantable > sig2;
1053 CountedPtr< Scantable > ref2;
1054 CountedPtr< Scantable > calb1;
1055 CountedPtr< Scantable > calb2;
1056
1057 msg=String("Processing dototalpower for subset of the data");
1058 ostringstream oss1;
1059 oss1 << msg << endl;
1060 pushLog(String(oss1));
1061 // Debug for IRC CS data
1062 //float tcal1=7.0;
1063 //float tcal2=4.0;
1064 sig1 = dototalpower(sctables[0], sctables[1], tcal=tcal);
1065 ref1 = dototalpower(sctables[2], sctables[3], tcal=tcal);
1066 ref2 = dototalpower(sctables[4], sctables[5], tcal=tcal);
1067 sig2 = dototalpower(sctables[6], sctables[7], tcal=tcal);
1068
1069 // correction of user-specified tsys for elevation here
1070
1071 // dosigref calibration
1072 msg=String("Processing dosigref for subset of the data");
1073 ostringstream oss2;
1074 oss2 << msg << endl;
1075 pushLog(String(oss2));
1076 calb1=dosigref(sig1,ref2,smoothref,tsysv,tau);
1077 calb2=dosigref(sig2,ref1,smoothref,tsysv,tau);
1078
1079 // iteration by scanno or cycle no.
1080 Table& tcalb1 = calb1->table();
1081 Table& tcalb2 = calb2->table();
1082 TableIterator sit(tcalb1, "SCANNO");
1083 TableIterator s2it(tcalb2, "SCANNO");
1084 while ( !sit.pastEnd() ) {
1085 Table t1 = sit.table();
1086 Table t2= s2it.table();
1087 ArrayColumn<Float> outspecCol(t1, "SPECTRA");
1088 ArrayColumn<Float> outtsysCol(t1, "TSYS");
1089 ArrayColumn<uChar> outflagCol(t1, "FLAGTRA");
1090 ScalarColumn<Double> outintCol(t1, "INTERVAL");
1091 ArrayColumn<Float> t2specCol(t2, "SPECTRA");
1092 ROArrayColumn<Float> t2tsysCol(t2, "TSYS");
1093 ArrayColumn<uChar> t2flagCol(t2, "FLAGTRA");
1094 ROScalarColumn<Double> t2intCol(t2, "INTERVAL");
1095 for (uInt i=0; i < t1.nrow(); ++i) {
1096 Vector<Float> spec1, spec2;
1097 // to store scalar (mean) tsys
1098 Vector<Float> tsys1, tsys2;
1099 Vector<uChar> flag1, flag2;
1100 Double tint1, tint2;
1101 outspecCol.get(i, spec1);
1102 t2specCol.get(i, spec2);
1103 outflagCol.get(i, flag1);
1104 t2flagCol.get(i, flag2);
1105 outtsysCol.get(i, tsys1);
1106 t2tsysCol.get(i, tsys2);
1107 outintCol.get(i, tint1);
1108 t2intCol.get(i, tint2);
1109 // average
1110 // assume scalar tsys for weights
1111 Float wt1, wt2, tsyssq1, tsyssq2;
1112 tsyssq1 = tsys1[0]*tsys1[0];
1113 tsyssq2 = tsys2[0]*tsys2[0];
1114 wt1 = Float(tint1)/tsyssq1;
1115 wt2 = Float(tint2)/tsyssq2;
1116 Float invsumwt=1/(wt1+wt2);
1117 MaskedArray<Float> mspec1 = maskedArray(spec1, flag1);
1118 MaskedArray<Float> mspec2 = maskedArray(spec2, flag2);
1119 MaskedArray<Float> avspec = invsumwt * (wt1*mspec1 + wt2*mspec2);
1120 //Array<Float> avtsys = Float(0.5) * (tsys1 + tsys2);
1121 // cerr<< "Tsys1="<<tsys1<<" Tsys2="<<tsys2<<endl;
1122 tsys1[0] = sqrt(tsyssq1 + tsyssq2);
1123 Array<Float> avtsys = tsys1;
1124
1125 outspecCol.put(i, avspec.getArray());
1126 outflagCol.put(i, flagsFromMA(avspec));
1127 outtsysCol.put(i, avtsys);
1128 }
1129 ++sit;
1130 ++s2it;
1131 }
1132 return calb1;
1133}
1134
1135//GBTIDL version of frequency switched data calibration
1136CountedPtr< Scantable > STMath::dofs( const CountedPtr< Scantable >& s,
1137 const std::vector<int>& scans,
1138 int smoothref,
1139 casa::Float tsysv,
1140 casa::Float tau,
1141 casa::Float tcal )
1142{
1143
1144
1145 STSelector sel;
1146 CountedPtr< Scantable > ws = getScantable(s, false);
1147 CountedPtr< Scantable > sig, sigwcal, ref, refwcal;
1148 CountedPtr< Scantable > calsig, calref, out, out1, out2;
1149 Bool nofold=False;
1150
1151 //split the data
1152 sel.setName("*_fs");
1153 ws->setSelection(sel);
1154 sig = getScantable(ws,false);
1155 sel.reset();
1156 sel.setName("*_fs_calon");
1157 ws->setSelection(sel);
1158 sigwcal = getScantable(ws,false);
1159 sel.reset();
1160 sel.setName("*_fsr");
1161 ws->setSelection(sel);
1162 ref = getScantable(ws,false);
1163 sel.reset();
1164 sel.setName("*_fsr_calon");
1165 ws->setSelection(sel);
1166 refwcal = getScantable(ws,false);
1167
1168 calsig = dototalpower(sigwcal, sig, tcal=tcal);
1169 calref = dototalpower(refwcal, ref, tcal=tcal);
1170
1171 out1=dosigref(calsig,calref,smoothref,tsysv,tau);
1172 out2=dosigref(calref,calsig,smoothref,tsysv,tau);
1173
1174 Table& tabout1=out1->table();
1175 Table& tabout2=out2->table();
1176 ROScalarColumn<uInt> freqidCol1(tabout1, "FREQ_ID");
1177 ScalarColumn<uInt> freqidCol2(tabout2, "FREQ_ID");
1178 ROArrayColumn<Float> specCol(tabout2, "SPECTRA");
1179 Vector<Float> spec; specCol.get(0, spec);
1180 uInt nchan = spec.nelements();
1181 uInt freqid1; freqidCol1.get(0,freqid1);
1182 uInt freqid2; freqidCol2.get(0,freqid2);
1183 Double rp1, rp2, rv1, rv2, inc1, inc2;
1184 out1->frequencies().getEntry(rp1, rv1, inc1, freqid1);
1185 out2->frequencies().getEntry(rp2, rv2, inc2, freqid2);
1186 //cerr << out1->frequencies().table().nrow() << " " << out2->frequencies().table().nrow() << endl ;
1187 if (rp1==rp2) {
1188 Double foffset = rv1 - rv2;
1189 uInt choffset = static_cast<uInt>(foffset/abs(inc2));
1190 if (choffset >= nchan) {
1191 cerr<<"out-band frequency switching, no folding"<<endl;
1192 nofold = True;
1193 }
1194 }
1195
1196 if (nofold) {
1197 std::vector< CountedPtr< Scantable > > tabs;
1198 tabs.push_back(out1);
1199 tabs.push_back(out2);
1200 out = merge(tabs);
1201 }
1202 else { //folding is not implemented yet
1203 //out = out1;
1204 Int choffset = static_cast<Int>((rv1-rv2)/inc2) ;
1205 out = dofold( out1, out2, choffset ) ;
1206 }
1207
1208 return out;
1209}
1210
1211CountedPtr<Scantable> STMath::dofold( const CountedPtr<Scantable> &sig,
1212 const CountedPtr<Scantable> &ref,
1213 Int choffset )
1214{
1215 // output scantable
1216 CountedPtr<Scantable> out = getScantable( sig, false ) ;
1217
1218 // get column
1219 ROArrayColumn<Float> specCol1( sig->table(), "SPECTRA" ) ;
1220 ROArrayColumn<Float> specCol2( ref->table(), "SPECTRA" ) ;
1221 ROArrayColumn<Float> tsysCol1( sig->table(), "TSYS" ) ;
1222 ROArrayColumn<Float> tsysCol2( ref->table(), "TSYS" ) ;
1223 ROArrayColumn<uChar> flagCol1( sig->table(), "FLAGTRA" ) ;
1224 ROArrayColumn<uChar> flagCol2( ref->table(), "FLAGTRA" ) ;
1225 ROScalarColumn<Double> mjdCol1( sig->table(), "TIME" ) ;
1226 ROScalarColumn<Double> mjdCol2( ref->table(), "TIME" ) ;
1227 ROScalarColumn<Double> intervalCol1( sig->table(), "INTERVAL" ) ;
1228 ROScalarColumn<Double> intervalCol2( ref->table(), "INTERVAL" ) ;
1229
1230 // check
1231 if ( choffset == 0 ) {
1232 cerr << "channel offset is zero, no folding" << endl ;
1233 return out ;
1234 }
1235 int nchan = ref->nchan() ;
1236 if ( abs(choffset) >= nchan ) {
1237 cerr << "out-band frequency switching, no folding" << endl ;
1238 return out ;
1239 }
1240
1241 // attach column for output scantable
1242 ArrayColumn<Float> specColOut( out->table(), "SPECTRA" ) ;
1243 ArrayColumn<uChar> flagColOut( out->table(), "FLAGTRA" ) ;
1244 ArrayColumn<Float> tsysColOut( out->table(), "TSYS" ) ;
1245 ScalarColumn<Double> mjdColOut( out->table(), "TIME" ) ;
1246 ScalarColumn<Double> intervalColOut( out->table(), "INTERVAL" ) ;
1247
1248 // for each row
1249 // assume that the data order are same between sig and ref
1250 RowAccumulator acc( asap::TINTSYS ) ;
1251 for ( int i = 0 ; i < sig->nrow() ; i++ ) {
1252 // get values
1253 Vector<Float> spsig ;
1254 specCol1.get( i, spsig ) ;
1255 Vector<Float> spref ;
1256 specCol2.get( i, spref ) ;
1257 Vector<Float> tsyssig ;
1258 tsysCol1.get( i, tsyssig ) ;
1259 Vector<Float> tsysref ;
1260 tsysCol2.get( i, tsysref ) ;
1261 Vector<uChar> flagsig ;
1262 flagCol1.get( i, flagsig ) ;
1263 Vector<uChar> flagref ;
1264 flagCol2.get( i, flagref ) ;
1265 Double timesig ;
1266 mjdCol1.get( i, timesig ) ;
1267 Double timeref ;
1268 mjdCol2.get( i, timeref ) ;
1269 Double intsig ;
1270 intervalCol1.get( i, intsig ) ;
1271 Double intref ;
1272 intervalCol2.get( i, intref ) ;
1273
1274 // shift reference spectra
1275 int refchan = spref.nelements() ;
1276 if ( choffset > 0 ) {
1277 for ( int j = 0 ; j < refchan-choffset ; j++ ) {
1278 spref[j] = spref[j+choffset] ;
1279 tsysref[j] = tsysref[j+choffset] ;
1280 flagref[j] = flagref[j+choffset] ;
1281 }
1282 for ( int j = refchan-choffset ; j < refchan ; j++ ) {
1283 spref[j] = spref[j-refchan+choffset] ;
1284 tsysref[j] = tsysref[j-refchan+choffset] ;
1285 flagref[j] = flagref[j-refchan+choffset] ;
1286 }
1287 }
1288 else {
1289 for ( int j = 0 ; j < abs(choffset) ; j++ ) {
1290 spref[j] = spref[refchan+choffset+j] ;
1291 tsysref[j] = tsysref[refchan+choffset+j] ;
1292 flagref[j] = flagref[refchan+choffset+j] ;
1293 }
1294 for ( int j = abs(choffset) ; j < refchan ; j++ ) {
1295 spref[j] = spref[j+choffset] ;
1296 tsysref[j] = tsysref[j+choffset] ;
1297 flagref[j] = flagref[j+choffset] ;
1298 }
1299 }
1300
1301 // folding
1302 acc.add( spsig, !flagsig, tsyssig, intsig, timesig ) ;
1303 acc.add( spref, !flagref, tsysref, intref, timeref ) ;
1304
1305 // put result
1306 specColOut.put( i, acc.getSpectrum() ) ;
1307 const Vector<Bool> &msk = acc.getMask() ;
1308 Vector<uChar> flg( msk.shape() ) ;
1309 convertArray( flg, !msk ) ;
1310 flagColOut.put( i, flg ) ;
1311 tsysColOut.put( i, acc.getTsys() ) ;
1312 intervalColOut.put( i, acc.getInterval() ) ;
1313 mjdColOut.put( i, acc.getTime() ) ;
1314
1315 acc.reset() ;
1316 }
1317
1318 return out ;
1319}
1320
1321
1322CountedPtr< Scantable > STMath::freqSwitch( const CountedPtr< Scantable >& in )
1323{
1324 // make copy or reference
1325 CountedPtr< Scantable > out = getScantable(in, false);
1326 Table& tout = out->table();
1327 Block<String> cols(4);
1328 cols[0] = String("SCANNO");
1329 cols[1] = String("CYCLENO");
1330 cols[2] = String("BEAMNO");
1331 cols[3] = String("POLNO");
1332 TableIterator iter(tout, cols);
1333 while (!iter.pastEnd()) {
1334 Table subt = iter.table();
1335 // this should leave us with two rows for the two IFs....if not ignore
1336 if (subt.nrow() != 2 ) {
1337 continue;
1338 }
1339 ArrayColumn<Float> specCol(subt, "SPECTRA");
1340 ArrayColumn<Float> tsysCol(subt, "TSYS");
1341 ArrayColumn<uChar> flagCol(subt, "FLAGTRA");
1342 Vector<Float> onspec,offspec, ontsys, offtsys;
1343 Vector<uChar> onflag, offflag;
1344 tsysCol.get(0, ontsys); tsysCol.get(1, offtsys);
1345 specCol.get(0, onspec); specCol.get(1, offspec);
1346 flagCol.get(0, onflag); flagCol.get(1, offflag);
1347 MaskedArray<Float> on = maskedArray(onspec, onflag);
1348 MaskedArray<Float> off = maskedArray(offspec, offflag);
1349 MaskedArray<Float> oncopy = on.copy();
1350
1351 on /= off; on -= 1.0f;
1352 on *= ontsys[0];
1353 off /= oncopy; off -= 1.0f;
1354 off *= offtsys[0];
1355 specCol.put(0, on.getArray());
1356 const Vector<Bool>& m0 = on.getMask();
1357 Vector<uChar> flags0(m0.shape());
1358 convertArray(flags0, !m0);
1359 flagCol.put(0, flags0);
1360
1361 specCol.put(1, off.getArray());
1362 const Vector<Bool>& m1 = off.getMask();
1363 Vector<uChar> flags1(m1.shape());
1364 convertArray(flags1, !m1);
1365 flagCol.put(1, flags1);
1366 ++iter;
1367 }
1368
1369 return out;
1370}
1371
1372std::vector< float > STMath::statistic( const CountedPtr< Scantable > & in,
1373 const std::vector< bool > & mask,
1374 const std::string& which )
1375{
1376
1377 Vector<Bool> m(mask);
1378 const Table& tab = in->table();
1379 ROArrayColumn<Float> specCol(tab, "SPECTRA");
1380 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1381 std::vector<float> out;
1382 for (uInt i=0; i < tab.nrow(); ++i ) {
1383 Vector<Float> spec; specCol.get(i, spec);
1384 Vector<uChar> flag; flagCol.get(i, flag);
1385 MaskedArray<Float> ma = maskedArray(spec, flag);
1386 float outstat = 0.0;
1387 if ( spec.nelements() == m.nelements() ) {
1388 outstat = mathutil::statistics(which, ma(m));
1389 } else {
1390 outstat = mathutil::statistics(which, ma);
1391 }
1392 out.push_back(outstat);
1393 }
1394 return out;
1395}
1396
1397std::vector< int > STMath::minMaxChan( const CountedPtr< Scantable > & in,
1398 const std::vector< bool > & mask,
1399 const std::string& which )
1400{
1401
1402 Vector<Bool> m(mask);
1403 const Table& tab = in->table();
1404 ROArrayColumn<Float> specCol(tab, "SPECTRA");
1405 ROArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1406 std::vector<int> out;
1407 for (uInt i=0; i < tab.nrow(); ++i ) {
1408 Vector<Float> spec; specCol.get(i, spec);
1409 Vector<uChar> flag; flagCol.get(i, flag);
1410 MaskedArray<Float> ma = maskedArray(spec, flag);
1411 if (ma.ndim() != 1) {
1412 throw (ArrayError(
1413 "std::vector<int> STMath::minMaxChan("
1414 "ContedPtr<Scantable> &in, std::vector<bool> &mask, "
1415 " std::string &which)"
1416 " - MaskedArray is not 1D"));
1417 }
1418 IPosition outpos(1,0);
1419 if ( spec.nelements() == m.nelements() ) {
1420 outpos = mathutil::minMaxPos(which, ma(m));
1421 } else {
1422 outpos = mathutil::minMaxPos(which, ma);
1423 }
1424 out.push_back(outpos[0]);
1425 }
1426 return out;
1427}
1428
1429CountedPtr< Scantable > STMath::bin( const CountedPtr< Scantable > & in,
1430 int width )
1431{
1432 if ( !in->getSelection().empty() ) throw(AipsError("Can't bin subset of the data."));
1433 CountedPtr< Scantable > out = getScantable(in, false);
1434 Table& tout = out->table();
1435 out->frequencies().rescale(width, "BIN");
1436 ArrayColumn<Float> specCol(tout, "SPECTRA");
1437 ArrayColumn<uChar> flagCol(tout, "FLAGTRA");
1438 for (uInt i=0; i < tout.nrow(); ++i ) {
1439 MaskedArray<Float> main = maskedArray(specCol(i), flagCol(i));
1440 MaskedArray<Float> maout;
1441 LatticeUtilities::bin(maout, main, 0, Int(width));
1442 /// @todo implement channel based tsys binning
1443 specCol.put(i, maout.getArray());
1444 flagCol.put(i, flagsFromMA(maout));
1445 // take only the first binned spectrum's length for the deprecated
1446 // global header item nChan
1447 if (i==0) tout.rwKeywordSet().define(String("nChan"),
1448 Int(maout.getArray().nelements()));
1449 }
1450 return out;
1451}
1452
1453CountedPtr< Scantable > STMath::resample( const CountedPtr< Scantable >& in,
1454 const std::string& method,
1455 float width )
1456//
1457// Should add the possibility of width being specified in km/s. This means
1458// that for each freqID (SpectralCoordinate) we will need to convert to an
1459// average channel width (say at the reference pixel). Then we would need
1460// to be careful to make sure each spectrum (of different freqID)
1461// is the same length.
1462//
1463{
1464 //InterpolateArray1D<Double,Float>::InterpolationMethod interp;
1465 Int interpMethod(stringToIMethod(method));
1466
1467 CountedPtr< Scantable > out = getScantable(in, false);
1468 Table& tout = out->table();
1469
1470// Resample SpectralCoordinates (one per freqID)
1471 out->frequencies().rescale(width, "RESAMPLE");
1472 TableIterator iter(tout, "IFNO");
1473 TableRow row(tout);
1474 while ( !iter.pastEnd() ) {
1475 Table tab = iter.table();
1476 ArrayColumn<Float> specCol(tab, "SPECTRA");
1477 //ArrayColumn<Float> tsysCol(tout, "TSYS");
1478 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1479 Vector<Float> spec;
1480 Vector<uChar> flag;
1481 specCol.get(0,spec); // the number of channels should be constant per IF
1482 uInt nChanIn = spec.nelements();
1483 Vector<Float> xIn(nChanIn); indgen(xIn);
1484 Int fac = Int(nChanIn/width);
1485 Vector<Float> xOut(fac+10); // 10 to be safe - resize later
1486 uInt k = 0;
1487 Float x = 0.0;
1488 while (x < Float(nChanIn) ) {
1489 xOut(k) = x;
1490 k++;
1491 x += width;
1492 }
1493 uInt nChanOut = k;
1494 xOut.resize(nChanOut, True);
1495 // process all rows for this IFNO
1496 Vector<Float> specOut;
1497 Vector<Bool> maskOut;
1498 Vector<uChar> flagOut;
1499 for (uInt i=0; i < tab.nrow(); ++i) {
1500 specCol.get(i, spec);
1501 flagCol.get(i, flag);
1502 Vector<Bool> mask(flag.nelements());
1503 convertArray(mask, flag);
1504
1505 IPosition shapeIn(spec.shape());
1506 //sh.nchan = nChanOut;
1507 InterpolateArray1D<Float,Float>::interpolate(specOut, maskOut, xOut,
1508 xIn, spec, mask,
1509 interpMethod, True, True);
1510 /// @todo do the same for channel based Tsys
1511 flagOut.resize(maskOut.nelements());
1512 convertArray(flagOut, maskOut);
1513 specCol.put(i, specOut);
1514 flagCol.put(i, flagOut);
1515 }
1516 ++iter;
1517 }
1518
1519 return out;
1520}
1521
1522STMath::imethod STMath::stringToIMethod(const std::string& in)
1523{
1524 static STMath::imap lookup;
1525
1526 // initialize the lookup table if necessary
1527 if ( lookup.empty() ) {
1528 lookup["nearest"] = InterpolateArray1D<Double,Float>::nearestNeighbour;
1529 lookup["linear"] = InterpolateArray1D<Double,Float>::linear;
1530 lookup["cubic"] = InterpolateArray1D<Double,Float>::cubic;
1531 lookup["spline"] = InterpolateArray1D<Double,Float>::spline;
1532 }
1533
1534 STMath::imap::const_iterator iter = lookup.find(in);
1535
1536 if ( lookup.end() == iter ) {
1537 std::string message = in;
1538 message += " is not a valid interpolation mode";
1539 throw(AipsError(message));
1540 }
1541 return iter->second;
1542}
1543
1544WeightType STMath::stringToWeight(const std::string& in)
1545{
1546 static std::map<std::string, WeightType> lookup;
1547
1548 // initialize the lookup table if necessary
1549 if ( lookup.empty() ) {
1550 lookup["NONE"] = asap::NONE;
1551 lookup["TINT"] = asap::TINT;
1552 lookup["TINTSYS"] = asap::TINTSYS;
1553 lookup["TSYS"] = asap::TSYS;
1554 lookup["VAR"] = asap::VAR;
1555 }
1556
1557 std::map<std::string, WeightType>::const_iterator iter = lookup.find(in);
1558
1559 if ( lookup.end() == iter ) {
1560 std::string message = in;
1561 message += " is not a valid weighting mode";
1562 throw(AipsError(message));
1563 }
1564 return iter->second;
1565}
1566
1567CountedPtr< Scantable > STMath::gainElevation( const CountedPtr< Scantable >& in,
1568 const vector< float > & coeff,
1569 const std::string & filename,
1570 const std::string& method)
1571{
1572 // Get elevation data from Scantable and convert to degrees
1573 CountedPtr< Scantable > out = getScantable(in, false);
1574 Table& tab = out->table();
1575 ROScalarColumn<Float> elev(tab, "ELEVATION");
1576 Vector<Float> x = elev.getColumn();
1577 x *= Float(180 / C::pi); // Degrees
1578
1579 Vector<Float> coeffs(coeff);
1580 const uInt nc = coeffs.nelements();
1581 if ( filename.length() > 0 && nc > 0 ) {
1582 throw(AipsError("You must choose either polynomial coefficients or an ascii file, not both"));
1583 }
1584
1585 // Correct
1586 if ( nc > 0 || filename.length() == 0 ) {
1587 // Find instrument
1588 Bool throwit = True;
1589 Instrument inst =
1590 STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
1591 throwit);
1592
1593 // Set polynomial
1594 Polynomial<Float>* ppoly = 0;
1595 Vector<Float> coeff;
1596 String msg;
1597 if ( nc > 0 ) {
1598 ppoly = new Polynomial<Float>(nc);
1599 coeff = coeffs;
1600 msg = String("user");
1601 } else {
1602 STAttr sdAttr;
1603 coeff = sdAttr.gainElevationPoly(inst);
1604 ppoly = new Polynomial<Float>(3);
1605 msg = String("built in");
1606 }
1607
1608 if ( coeff.nelements() > 0 ) {
1609 ppoly->setCoefficients(coeff);
1610 } else {
1611 delete ppoly;
1612 throw(AipsError("There is no known gain-elevation polynomial known for this instrument"));
1613 }
1614 ostringstream oss;
1615 oss << "Making polynomial correction with " << msg << " coefficients:" << endl;
1616 oss << " " << coeff;
1617 pushLog(String(oss));
1618 const uInt nrow = tab.nrow();
1619 Vector<Float> factor(nrow);
1620 for ( uInt i=0; i < nrow; ++i ) {
1621 factor[i] = 1.0 / (*ppoly)(x[i]);
1622 }
1623 delete ppoly;
1624 scaleByVector(tab, factor, true);
1625
1626 } else {
1627 // Read and correct
1628 pushLog("Making correction from ascii Table");
1629 scaleFromAsciiTable(tab, filename, method, x, true);
1630 }
1631 return out;
1632}
1633
1634void STMath::scaleFromAsciiTable(Table& in, const std::string& filename,
1635 const std::string& method,
1636 const Vector<Float>& xout, bool dotsys)
1637{
1638
1639// Read gain-elevation ascii file data into a Table.
1640
1641 String formatString;
1642 Table tbl = readAsciiTable(formatString, Table::Memory, filename, "", "", False);
1643 scaleFromTable(in, tbl, method, xout, dotsys);
1644}
1645
1646void STMath::scaleFromTable(Table& in,
1647 const Table& table,
1648 const std::string& method,
1649 const Vector<Float>& xout, bool dotsys)
1650{
1651
1652 ROScalarColumn<Float> geElCol(table, "ELEVATION");
1653 ROScalarColumn<Float> geFacCol(table, "FACTOR");
1654 Vector<Float> xin = geElCol.getColumn();
1655 Vector<Float> yin = geFacCol.getColumn();
1656 Vector<Bool> maskin(xin.nelements(),True);
1657
1658 // Interpolate (and extrapolate) with desired method
1659
1660 InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
1661
1662 Vector<Float> yout;
1663 Vector<Bool> maskout;
1664 InterpolateArray1D<Float,Float>::interpolate(yout, maskout, xout,
1665 xin, yin, maskin, interp,
1666 True, True);
1667
1668 scaleByVector(in, Float(1.0)/yout, dotsys);
1669}
1670
1671void STMath::scaleByVector( Table& in,
1672 const Vector< Float >& factor,
1673 bool dotsys )
1674{
1675 uInt nrow = in.nrow();
1676 if ( factor.nelements() != nrow ) {
1677 throw(AipsError("factors.nelements() != table.nelements()"));
1678 }
1679 ArrayColumn<Float> specCol(in, "SPECTRA");
1680 ArrayColumn<uChar> flagCol(in, "FLAGTRA");
1681 ArrayColumn<Float> tsysCol(in, "TSYS");
1682 for (uInt i=0; i < nrow; ++i) {
1683 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
1684 ma *= factor[i];
1685 specCol.put(i, ma.getArray());
1686 flagCol.put(i, flagsFromMA(ma));
1687 if ( dotsys ) {
1688 Vector<Float> tsys = tsysCol(i);
1689 tsys *= factor[i];
1690 tsysCol.put(i,tsys);
1691 }
1692 }
1693}
1694
1695CountedPtr< Scantable > STMath::convertFlux( const CountedPtr< Scantable >& in,
1696 float d, float etaap,
1697 float jyperk )
1698{
1699 CountedPtr< Scantable > out = getScantable(in, false);
1700 Table& tab = in->table();
1701 Unit fluxUnit(tab.keywordSet().asString("FluxUnit"));
1702 Unit K(String("K"));
1703 Unit JY(String("Jy"));
1704
1705 bool tokelvin = true;
1706 Double cfac = 1.0;
1707
1708 if ( fluxUnit == JY ) {
1709 pushLog("Converting to K");
1710 Quantum<Double> t(1.0,fluxUnit);
1711 Quantum<Double> t2 = t.get(JY);
1712 cfac = (t2 / t).getValue(); // value to Jy
1713
1714 tokelvin = true;
1715 out->setFluxUnit("K");
1716 } else if ( fluxUnit == K ) {
1717 pushLog("Converting to Jy");
1718 Quantum<Double> t(1.0,fluxUnit);
1719 Quantum<Double> t2 = t.get(K);
1720 cfac = (t2 / t).getValue(); // value to K
1721
1722 tokelvin = false;
1723 out->setFluxUnit("Jy");
1724 } else {
1725 throw(AipsError("Unrecognized brightness units in Table - must be consistent with Jy or K"));
1726 }
1727 // Make sure input values are converted to either Jy or K first...
1728 Float factor = cfac;
1729
1730 // Select method
1731 if (jyperk > 0.0) {
1732 factor *= jyperk;
1733 if ( tokelvin ) factor = 1.0 / jyperk;
1734 ostringstream oss;
1735 oss << "Jy/K = " << jyperk;
1736 pushLog(String(oss));
1737 Vector<Float> factors(tab.nrow(), factor);
1738 scaleByVector(tab,factors, false);
1739 } else if ( etaap > 0.0) {
1740 if (d < 0) {
1741 Instrument inst =
1742 STAttr::convertInstrument(tab.keywordSet().asString("AntennaName"),
1743 True);
1744 STAttr sda;
1745 d = sda.diameter(inst);
1746 }
1747 jyperk = STAttr::findJyPerK(etaap, d);
1748 ostringstream oss;
1749 oss << "Jy/K = " << jyperk;
1750 pushLog(String(oss));
1751 factor *= jyperk;
1752 if ( tokelvin ) {
1753 factor = 1.0 / factor;
1754 }
1755 Vector<Float> factors(tab.nrow(), factor);
1756 scaleByVector(tab, factors, False);
1757 } else {
1758
1759 // OK now we must deal with automatic look up of values.
1760 // We must also deal with the fact that the factors need
1761 // to be computed per IF and may be different and may
1762 // change per integration.
1763
1764 pushLog("Looking up conversion factors");
1765 convertBrightnessUnits(out, tokelvin, cfac);
1766 }
1767
1768 return out;
1769}
1770
1771void STMath::convertBrightnessUnits( CountedPtr<Scantable>& in,
1772 bool tokelvin, float cfac )
1773{
1774 Table& table = in->table();
1775 Instrument inst =
1776 STAttr::convertInstrument(table.keywordSet().asString("AntennaName"), True);
1777 TableIterator iter(table, "FREQ_ID");
1778 STFrequencies stfreqs = in->frequencies();
1779 STAttr sdAtt;
1780 while (!iter.pastEnd()) {
1781 Table tab = iter.table();
1782 ArrayColumn<Float> specCol(tab, "SPECTRA");
1783 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1784 ROScalarColumn<uInt> freqidCol(tab, "FREQ_ID");
1785 MEpoch::ROScalarColumn timeCol(tab, "TIME");
1786
1787 uInt freqid; freqidCol.get(0, freqid);
1788 Vector<Float> tmpspec; specCol.get(0, tmpspec);
1789 // STAttr.JyPerK has a Vector interface... change sometime.
1790 Vector<Float> freqs(1,stfreqs.getRefFreq(freqid, tmpspec.nelements()));
1791 for ( uInt i=0; i<tab.nrow(); ++i) {
1792 Float jyperk = (sdAtt.JyPerK(inst, timeCol(i), freqs))[0];
1793 Float factor = cfac * jyperk;
1794 if ( tokelvin ) factor = Float(1.0) / factor;
1795 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
1796 ma *= factor;
1797 specCol.put(i, ma.getArray());
1798 flagCol.put(i, flagsFromMA(ma));
1799 }
1800 ++iter;
1801 }
1802}
1803
1804CountedPtr< Scantable > STMath::opacity( const CountedPtr< Scantable > & in,
1805 float tau )
1806{
1807 CountedPtr< Scantable > out = getScantable(in, false);
1808
1809 Table tab = out->table();
1810 ROScalarColumn<Float> elev(tab, "ELEVATION");
1811 ArrayColumn<Float> specCol(tab, "SPECTRA");
1812 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1813 ArrayColumn<Float> tsysCol(tab, "TSYS");
1814 for ( uInt i=0; i<tab.nrow(); ++i) {
1815 Float zdist = Float(C::pi_2) - elev(i);
1816 Float factor = exp(tau/cos(zdist));
1817 MaskedArray<Float> ma = maskedArray(specCol(i), flagCol(i));
1818 ma *= factor;
1819 specCol.put(i, ma.getArray());
1820 flagCol.put(i, flagsFromMA(ma));
1821 Vector<Float> tsys;
1822 tsysCol.get(i, tsys);
1823 tsys *= factor;
1824 tsysCol.put(i, tsys);
1825 }
1826 return out;
1827}
1828
1829CountedPtr< Scantable > STMath::smoothOther( const CountedPtr< Scantable >& in,
1830 const std::string& kernel,
1831 float width )
1832{
1833 CountedPtr< Scantable > out = getScantable(in, false);
1834 Table& table = out->table();
1835 ArrayColumn<Float> specCol(table, "SPECTRA");
1836 ArrayColumn<uChar> flagCol(table, "FLAGTRA");
1837 Vector<Float> spec;
1838 Vector<uChar> flag;
1839 for ( uInt i=0; i<table.nrow(); ++i) {
1840 specCol.get(i, spec);
1841 flagCol.get(i, flag);
1842 Vector<Bool> mask(flag.nelements());
1843 convertArray(mask, flag);
1844 Vector<Float> specout;
1845 Vector<Bool> maskout;
1846 if ( kernel == "hanning" ) {
1847 mathutil::hanning(specout, maskout, spec , !mask);
1848 convertArray(flag, !maskout);
1849 } else if ( kernel == "rmedian" ) {
1850 mathutil::runningMedian(specout, maskout, spec , mask, width);
1851 convertArray(flag, maskout);
1852 }
1853 flagCol.put(i, flag);
1854 specCol.put(i, specout);
1855 }
1856 return out;
1857}
1858
1859CountedPtr< Scantable > STMath::smooth( const CountedPtr< Scantable >& in,
1860 const std::string& kernel, float width )
1861{
1862 if (kernel == "rmedian" || kernel == "hanning") {
1863 return smoothOther(in, kernel, width);
1864 }
1865 CountedPtr< Scantable > out = getScantable(in, false);
1866 Table& table = out->table();
1867 VectorKernel::KernelTypes type = VectorKernel::toKernelType(kernel);
1868 // same IFNO should have same no of channels
1869 // this saves overhead
1870 TableIterator iter(table, "IFNO");
1871 while (!iter.pastEnd()) {
1872 Table tab = iter.table();
1873 ArrayColumn<Float> specCol(tab, "SPECTRA");
1874 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
1875 Vector<Float> tmpspec; specCol.get(0, tmpspec);
1876 uInt nchan = tmpspec.nelements();
1877 Vector<Float> kvec = VectorKernel::make(type, width, nchan, True, False);
1878 Convolver<Float> conv(kvec, IPosition(1,nchan));
1879 Vector<Float> spec;
1880 Vector<uChar> flag;
1881 for ( uInt i=0; i<tab.nrow(); ++i) {
1882 specCol.get(i, spec);
1883 flagCol.get(i, flag);
1884 Vector<Bool> mask(flag.nelements());
1885 convertArray(mask, flag);
1886 Vector<Float> specout;
1887 mathutil::replaceMaskByZero(specout, mask);
1888 conv.linearConv(specout, spec);
1889 specCol.put(i, specout);
1890 }
1891 ++iter;
1892 }
1893 return out;
1894}
1895
1896CountedPtr< Scantable >
1897 STMath::merge( const std::vector< CountedPtr < Scantable > >& in )
1898{
1899 if ( in.size() < 2 ) {
1900 throw(AipsError("Need at least two scantables to perform a merge."));
1901 }
1902 std::vector<CountedPtr < Scantable > >::const_iterator it = in.begin();
1903 bool insitu = insitu_;
1904 setInsitu(false);
1905 CountedPtr< Scantable > out = getScantable(*it, false);
1906 setInsitu(insitu);
1907 Table& tout = out->table();
1908 ScalarColumn<uInt> freqidcol(tout,"FREQ_ID"), molidcol(tout, "MOLECULE_ID");
1909 ScalarColumn<uInt> scannocol(tout,"SCANNO"), focusidcol(tout,"FOCUS_ID");
1910 // Renumber SCANNO to be 0-based
1911 Vector<uInt> scannos = scannocol.getColumn();
1912 uInt offset = min(scannos);
1913 scannos -= offset;
1914 scannocol.putColumn(scannos);
1915 uInt newscanno = max(scannos)+1;
1916 ++it;
1917 while ( it != in.end() ){
1918 if ( ! (*it)->conformant(*out) ) {
1919 // non conformant.
1920 pushLog(String("Warning: Can't merge scantables as header info differs."));
1921 }
1922 out->appendToHistoryTable((*it)->history());
1923 const Table& tab = (*it)->table();
1924 TableIterator scanit(tab, "SCANNO");
1925 while (!scanit.pastEnd()) {
1926 TableIterator freqit(scanit.table(), "FREQ_ID");
1927 while ( !freqit.pastEnd() ) {
1928 Table thetab = freqit.table();
1929 uInt nrow = tout.nrow();
1930 tout.addRow(thetab.nrow());
1931 TableCopy::copyRows(tout, thetab, nrow, 0, thetab.nrow());
1932 ROTableRow row(thetab);
1933 for ( uInt i=0; i<thetab.nrow(); ++i) {
1934 uInt k = nrow+i;
1935 scannocol.put(k, newscanno);
1936 const TableRecord& rec = row.get(i);
1937 Double rv,rp,inc;
1938 (*it)->frequencies().getEntry(rp, rv, inc, rec.asuInt("FREQ_ID"));
1939 uInt id;
1940 id = out->frequencies().addEntry(rp, rv, inc);
1941 freqidcol.put(k,id);
1942 //String name,fname;Double rf;
1943 Vector<String> name,fname;Vector<Double> rf;
1944 (*it)->molecules().getEntry(rf, name, fname, rec.asuInt("MOLECULE_ID"));
1945 id = out->molecules().addEntry(rf, name, fname);
1946 molidcol.put(k, id);
1947 Float frot,fax,ftan,fhand,fmount,fuser, fxy, fxyp;
1948 (*it)->focus().getEntry(fax, ftan, frot, fhand,
1949 fmount,fuser, fxy, fxyp,
1950 rec.asuInt("FOCUS_ID"));
1951 id = out->focus().addEntry(fax, ftan, frot, fhand,
1952 fmount,fuser, fxy, fxyp);
1953 focusidcol.put(k, id);
1954 }
1955 ++freqit;
1956 }
1957 ++newscanno;
1958 ++scanit;
1959 }
1960 ++it;
1961 }
1962 return out;
1963}
1964
1965CountedPtr< Scantable >
1966 STMath::invertPhase( const CountedPtr < Scantable >& in )
1967{
1968 return applyToPol(in, &STPol::invertPhase, Float(0.0));
1969}
1970
1971CountedPtr< Scantable >
1972 STMath::rotateXYPhase( const CountedPtr < Scantable >& in, float phase )
1973{
1974 return applyToPol(in, &STPol::rotatePhase, Float(phase));
1975}
1976
1977CountedPtr< Scantable >
1978 STMath::rotateLinPolPhase( const CountedPtr < Scantable >& in, float phase )
1979{
1980 return applyToPol(in, &STPol::rotateLinPolPhase, Float(phase));
1981}
1982
1983CountedPtr< Scantable > STMath::applyToPol( const CountedPtr<Scantable>& in,
1984 STPol::polOperation fptr,
1985 Float phase )
1986{
1987 CountedPtr< Scantable > out = getScantable(in, false);
1988 Table& tout = out->table();
1989 Block<String> cols(4);
1990 cols[0] = String("SCANNO");
1991 cols[1] = String("BEAMNO");
1992 cols[2] = String("IFNO");
1993 cols[3] = String("CYCLENO");
1994 TableIterator iter(tout, cols);
1995 CountedPtr<STPol> stpol = STPol::getPolClass(out->factories_,
1996 out->getPolType() );
1997 while (!iter.pastEnd()) {
1998 Table t = iter.table();
1999 ArrayColumn<Float> speccol(t, "SPECTRA");
2000 ScalarColumn<uInt> focidcol(t, "FOCUS_ID");
2001 ScalarColumn<Float> parancol(t, "PARANGLE");
2002 Matrix<Float> pols(speccol.getColumn());
2003 try {
2004 stpol->setSpectra(pols);
2005 Float fang,fhand,parang;
2006 fang = in->focusTable_.getTotalFeedAngle(focidcol(0));
2007 fhand = in->focusTable_.getFeedHand(focidcol(0));
2008 parang = parancol(0);
2009 /// @todo re-enable this
2010 // disable total feed angle to support paralactifying Caswell style
2011 stpol->setPhaseCorrections(parang, -parang, fhand);
2012 // use a member function pointer in STPol. This only works on
2013 // the STPol pointer itself, not the Counted Pointer so
2014 // derefernce it.
2015 (&(*(stpol))->*fptr)(phase);
2016 speccol.putColumn(stpol->getSpectra());
2017 } catch (AipsError& e) {
2018 //delete stpol;stpol=0;
2019 throw(e);
2020 }
2021 ++iter;
2022 }
2023 //delete stpol;stpol=0;
2024 return out;
2025}
2026
2027CountedPtr< Scantable >
2028 STMath::swapPolarisations( const CountedPtr< Scantable > & in )
2029{
2030 CountedPtr< Scantable > out = getScantable(in, false);
2031 Table& tout = out->table();
2032 Table t0 = tout(tout.col("POLNO") == 0);
2033 Table t1 = tout(tout.col("POLNO") == 1);
2034 if ( t0.nrow() != t1.nrow() )
2035 throw(AipsError("Inconsistent number of polarisations"));
2036 ArrayColumn<Float> speccol0(t0, "SPECTRA");
2037 ArrayColumn<uChar> flagcol0(t0, "FLAGTRA");
2038 ArrayColumn<Float> speccol1(t1, "SPECTRA");
2039 ArrayColumn<uChar> flagcol1(t1, "FLAGTRA");
2040 Matrix<Float> s0 = speccol0.getColumn();
2041 Matrix<uChar> f0 = flagcol0.getColumn();
2042 speccol0.putColumn(speccol1.getColumn());
2043 flagcol0.putColumn(flagcol1.getColumn());
2044 speccol1.putColumn(s0);
2045 flagcol1.putColumn(f0);
2046 return out;
2047}
2048
2049CountedPtr< Scantable >
2050 STMath::averagePolarisations( const CountedPtr< Scantable > & in,
2051 const std::vector<bool>& mask,
2052 const std::string& weight )
2053{
2054 if (in->npol() < 2 )
2055 throw(AipsError("averagePolarisations can only be applied to two or more"
2056 "polarisations"));
2057 bool insitu = insitu_;
2058 setInsitu(false);
2059 CountedPtr< Scantable > pols = getScantable(in, true);
2060 setInsitu(insitu);
2061 Table& tout = pols->table();
2062 std::string taql = "SELECT FROM $1 WHERE POLNO IN [0,1]";
2063 Table tab = tableCommand(taql, in->table());
2064 if (tab.nrow() == 0 )
2065 throw(AipsError("Could not find any rows with POLNO==0 and POLNO==1"));
2066 TableCopy::copyRows(tout, tab);
2067 TableVector<uInt> vec(tout, "POLNO");
2068 vec = 0;
2069 pols->table_.rwKeywordSet().define("nPol", Int(1));
2070 //pols->table_.rwKeywordSet().define("POLTYPE", String("stokes"));
2071 pols->table_.rwKeywordSet().define("POLTYPE", in->getPolType());
2072 std::vector<CountedPtr<Scantable> > vpols;
2073 vpols.push_back(pols);
2074 CountedPtr< Scantable > out = average(vpols, mask, weight, "SCAN");
2075 return out;
2076}
2077
2078CountedPtr< Scantable >
2079 STMath::averageBeams( const CountedPtr< Scantable > & in,
2080 const std::vector<bool>& mask,
2081 const std::string& weight )
2082{
2083 bool insitu = insitu_;
2084 setInsitu(false);
2085 CountedPtr< Scantable > beams = getScantable(in, false);
2086 setInsitu(insitu);
2087 Table& tout = beams->table();
2088 // give all rows the same BEAMNO
2089 TableVector<uInt> vec(tout, "BEAMNO");
2090 vec = 0;
2091 beams->table_.rwKeywordSet().define("nBeam", Int(1));
2092 std::vector<CountedPtr<Scantable> > vbeams;
2093 vbeams.push_back(beams);
2094 CountedPtr< Scantable > out = average(vbeams, mask, weight, "SCAN");
2095 return out;
2096}
2097
2098
2099CountedPtr< Scantable >
2100 asap::STMath::frequencyAlign( const CountedPtr< Scantable > & in,
2101 const std::string & refTime,
2102 const std::string & method)
2103{
2104 // clone as this is not working insitu
2105 bool insitu = insitu_;
2106 setInsitu(false);
2107 CountedPtr< Scantable > out = getScantable(in, false);
2108 setInsitu(insitu);
2109 Table& tout = out->table();
2110 // Get reference Epoch to time of first row or given String
2111 Unit DAY(String("d"));
2112 MEpoch::Ref epochRef(in->getTimeReference());
2113 MEpoch refEpoch;
2114 if (refTime.length()>0) {
2115 Quantum<Double> qt;
2116 if (MVTime::read(qt,refTime)) {
2117 MVEpoch mv(qt);
2118 refEpoch = MEpoch(mv, epochRef);
2119 } else {
2120 throw(AipsError("Invalid format for Epoch string"));
2121 }
2122 } else {
2123 refEpoch = in->timeCol_(0);
2124 }
2125 MPosition refPos = in->getAntennaPosition();
2126
2127 InterpolateArray1D<Double,Float>::InterpolationMethod interp = stringToIMethod(method);
2128 /*
2129 // Comment from MV.
2130 // the following code has been commented out because different FREQ_IDs have to be aligned together even
2131 // if the frame doesn't change. So far, lack of this check didn't cause any problems.
2132 // test if user frame is different to base frame
2133 if ( in->frequencies().getFrameString(true)
2134 == in->frequencies().getFrameString(false) ) {
2135 throw(AipsError("Can't convert as no output frame has been set"
2136 " (use set_freqframe) or it is aligned already."));
2137 }
2138 */
2139 MFrequency::Types system = in->frequencies().getFrame();
2140 MVTime mvt(refEpoch.getValue());
2141 String epochout = mvt.string(MVTime::YMD) + String(" (") + refEpoch.getRefString() + String(")");
2142 ostringstream oss;
2143 oss << "Aligned at reference Epoch " << epochout
2144 << " in frame " << MFrequency::showType(system);
2145 pushLog(String(oss));
2146 // set up the iterator
2147 Block<String> cols(4);
2148 // select by constant direction
2149 cols[0] = String("SRCNAME");
2150 cols[1] = String("BEAMNO");
2151 // select by IF ( no of channels varies over this )
2152 cols[2] = String("IFNO");
2153 // select by restfrequency
2154 cols[3] = String("MOLECULE_ID");
2155 TableIterator iter(tout, cols);
2156 while ( !iter.pastEnd() ) {
2157 Table t = iter.table();
2158 MDirection::ROScalarColumn dirCol(t, "DIRECTION");
2159 TableIterator fiter(t, "FREQ_ID");
2160 // determine nchan from the first row. This should work as
2161 // we are iterating over BEAMNO and IFNO // we should have constant direction
2162
2163 ROArrayColumn<Float> sCol(t, "SPECTRA");
2164 const MDirection direction = dirCol(0);
2165 const uInt nchan = sCol(0).nelements();
2166
2167 // skip operations if there is nothing to align
2168 if (fiter.pastEnd()) {
2169 continue;
2170 }
2171
2172 Table ftab = fiter.table();
2173 // align all frequency ids with respect to the first encountered id
2174 ScalarColumn<uInt> freqidCol(ftab, "FREQ_ID");
2175 // get the SpectralCoordinate for the freqid, which we are iterating over
2176 SpectralCoordinate sC = in->frequencies().getSpectralCoordinate(freqidCol(0));
2177 FrequencyAligner<Float> fa( sC, nchan, refEpoch,
2178 direction, refPos, system );
2179 // realign the SpectralCoordinate and put into the output Scantable
2180 Vector<String> units(1);
2181 units = String("Hz");
2182 Bool linear=True;
2183 SpectralCoordinate sc2 = fa.alignedSpectralCoordinate(linear);
2184 sc2.setWorldAxisUnits(units);
2185 const uInt id = out->frequencies().addEntry(sc2.referencePixel()[0],
2186 sc2.referenceValue()[0],
2187 sc2.increment()[0]);
2188 while ( !fiter.pastEnd() ) {
2189 ftab = fiter.table();
2190 // spectral coordinate for the current FREQ_ID
2191 ScalarColumn<uInt> freqidCol2(ftab, "FREQ_ID");
2192 sC = in->frequencies().getSpectralCoordinate(freqidCol2(0));
2193 // create the "global" abcissa for alignment with same FREQ_ID
2194 Vector<Double> abc(nchan);
2195 for (uInt i=0; i<nchan; i++) {
2196 Double w;
2197 sC.toWorld(w,Double(i));
2198 abc[i] = w;
2199 }
2200 TableVector<uInt> tvec(ftab, "FREQ_ID");
2201 // assign new frequency id to all rows
2202 tvec = id;
2203 // cache abcissa for same time stamps, so iterate over those
2204 TableIterator timeiter(ftab, "TIME");
2205 while ( !timeiter.pastEnd() ) {
2206 Table tab = timeiter.table();
2207 ArrayColumn<Float> specCol(tab, "SPECTRA");
2208 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2209 MEpoch::ROScalarColumn timeCol(tab, "TIME");
2210 // use align abcissa cache after the first row
2211 // these rows should be just be POLNO
2212 bool first = true;
2213 for (int i=0; i<int(tab.nrow()); ++i) {
2214 // input values
2215 Vector<uChar> flag = flagCol(i);
2216 Vector<Bool> mask(flag.shape());
2217 Vector<Float> specOut, spec;
2218 spec = specCol(i);
2219 Vector<Bool> maskOut;Vector<uChar> flagOut;
2220 convertArray(mask, flag);
2221 // alignment
2222 Bool ok = fa.align(specOut, maskOut, abc, spec,
2223 mask, timeCol(i), !first,
2224 interp, False);
2225 // back into scantable
2226 flagOut.resize(maskOut.nelements());
2227 convertArray(flagOut, maskOut);
2228 flagCol.put(i, flagOut);
2229 specCol.put(i, specOut);
2230 // start abcissa caching
2231 first = false;
2232 }
2233 // next timestamp
2234 ++timeiter;
2235 }
2236 // next FREQ_ID
2237 ++fiter;
2238 }
2239 // next aligner
2240 ++iter;
2241 }
2242 // set this afterwards to ensure we are doing insitu correctly.
2243 out->frequencies().setFrame(system, true);
2244 return out;
2245}
2246
2247CountedPtr<Scantable>
2248 asap::STMath::convertPolarisation( const CountedPtr<Scantable>& in,
2249 const std::string & newtype )
2250{
2251 if (in->npol() != 2 && in->npol() != 4)
2252 throw(AipsError("Can only convert two or four polarisations."));
2253 if ( in->getPolType() == newtype )
2254 throw(AipsError("No need to convert."));
2255 if ( ! in->selector_.empty() )
2256 throw(AipsError("Can only convert whole scantable. Unset the selection."));
2257 bool insitu = insitu_;
2258 setInsitu(false);
2259 CountedPtr< Scantable > out = getScantable(in, true);
2260 setInsitu(insitu);
2261 Table& tout = out->table();
2262 tout.rwKeywordSet().define("POLTYPE", String(newtype));
2263
2264 Block<String> cols(4);
2265 cols[0] = "SCANNO";
2266 cols[1] = "CYCLENO";
2267 cols[2] = "BEAMNO";
2268 cols[3] = "IFNO";
2269 TableIterator it(in->originalTable_, cols);
2270 String basetype = in->getPolType();
2271 STPol* stpol = STPol::getPolClass(in->factories_, basetype);
2272 try {
2273 while ( !it.pastEnd() ) {
2274 Table tab = it.table();
2275 uInt row = tab.rowNumbers()[0];
2276 stpol->setSpectra(in->getPolMatrix(row));
2277 Float fang,fhand,parang;
2278 fang = in->focusTable_.getTotalFeedAngle(in->mfocusidCol_(row));
2279 fhand = in->focusTable_.getFeedHand(in->mfocusidCol_(row));
2280 parang = in->paraCol_(row);
2281 /// @todo re-enable this
2282 // disable total feed angle to support paralactifying Caswell style
2283 stpol->setPhaseCorrections(parang, -parang, fhand);
2284 Int npolout = 0;
2285 for (uInt i=0; i<tab.nrow(); ++i) {
2286 Vector<Float> outvec = stpol->getSpectrum(i, newtype);
2287 if ( outvec.nelements() > 0 ) {
2288 tout.addRow();
2289 TableCopy::copyRows(tout, tab, tout.nrow()-1, 0, 1);
2290 ArrayColumn<Float> sCol(tout,"SPECTRA");
2291 ScalarColumn<uInt> pCol(tout,"POLNO");
2292 sCol.put(tout.nrow()-1 ,outvec);
2293 pCol.put(tout.nrow()-1 ,uInt(npolout));
2294 npolout++;
2295 }
2296 }
2297 tout.rwKeywordSet().define("nPol", npolout);
2298 ++it;
2299 }
2300 } catch (AipsError& e) {
2301 delete stpol;
2302 throw(e);
2303 }
2304 delete stpol;
2305 return out;
2306}
2307
2308CountedPtr< Scantable >
2309 asap::STMath::mxExtract( const CountedPtr< Scantable > & in,
2310 const std::string & scantype )
2311{
2312 bool insitu = insitu_;
2313 setInsitu(false);
2314 CountedPtr< Scantable > out = getScantable(in, true);
2315 setInsitu(insitu);
2316 Table& tout = out->table();
2317 std::string taql = "SELECT FROM $1 WHERE BEAMNO != REFBEAMNO";
2318 if (scantype == "on") {
2319 taql = "SELECT FROM $1 WHERE BEAMNO == REFBEAMNO";
2320 }
2321 Table tab = tableCommand(taql, in->table());
2322 TableCopy::copyRows(tout, tab);
2323 if (scantype == "on") {
2324 // re-index SCANNO to 0
2325 TableVector<uInt> vec(tout, "SCANNO");
2326 vec = 0;
2327 }
2328 return out;
2329}
2330
2331CountedPtr< Scantable >
2332 asap::STMath::lagFlag( const CountedPtr< Scantable > & in,
2333 double frequency, double width )
2334{
2335 CountedPtr< Scantable > out = getScantable(in, false);
2336 Table& tout = out->table();
2337 TableIterator iter(tout, "FREQ_ID");
2338 FFTServer<Float,Complex> ffts;
2339 while ( !iter.pastEnd() ) {
2340 Table tab = iter.table();
2341 Double rp,rv,inc;
2342 ROTableRow row(tab);
2343 const TableRecord& rec = row.get(0);
2344 uInt freqid = rec.asuInt("FREQ_ID");
2345 out->frequencies().getEntry(rp, rv, inc, freqid);
2346 ArrayColumn<Float> specCol(tab, "SPECTRA");
2347 ArrayColumn<uChar> flagCol(tab, "FLAGTRA");
2348 for (int i=0; i<int(tab.nrow()); ++i) {
2349 Vector<Float> spec = specCol(i);
2350 Vector<uChar> flag = flagCol(i);
2351 Int lag0 = Int(spec.nelements()*abs(inc)/(frequency+width)+0.5);
2352 Int lag1 = Int(spec.nelements()*abs(inc)/(frequency-width)+0.5);
2353 for (int k=0; k < flag.nelements(); ++k ) {
2354 if (flag[k] > 0) {
2355 spec[k] = 0.0;
2356 }
2357 }
2358 Vector<Complex> lags;
2359 ffts.fft0(lags, spec);
2360 Int start = max(0, lag0);
2361 Int end = min(Int(lags.nelements()-1), lag1);
2362 if (start == end) {
2363 lags[start] = Complex(0.0);
2364 } else {
2365 for (int j=start; j <=end ;++j) {
2366 lags[j] = Complex(0.0);
2367 }
2368 }
2369 ffts.fft0(spec, lags);
2370 specCol.put(i, spec);
2371 }
2372 ++iter;
2373 }
2374 return out;
2375}
2376
2377// Averaging spectra with different channel/resolution
2378CountedPtr<Scantable>
2379STMath::new_average( const std::vector<CountedPtr<Scantable> >& in,
2380 const bool& compel,
2381 const std::vector<bool>& mask,
2382 const std::string& weight,
2383 const std::string& avmode )
2384 throw ( casa::AipsError )
2385{
2386 if ( avmode == "SCAN" && in.size() != 1 )
2387 throw(AipsError("Can't perform 'SCAN' averaging on multiple tables.\n"
2388 "Use merge first."));
2389
2390 // check if OTF observation
2391 String obstype = in[0]->getHeader().obstype ;
2392 //bool otfscan = false ;
2393 Double tol = 0.0 ;
2394 if ( obstype.find( "OTF" ) != String::npos ) {
2395 //cout << "OTF scan" << endl ;
2396 //otfscan = true ;
2397 tol = TOL_OTF ;
2398 }
2399 else {
2400 tol = TOL_POINT ;
2401 }
2402
2403 CountedPtr<Scantable> out ; // processed result
2404 if ( compel ) {
2405 std::vector< CountedPtr<Scantable> > newin ; // input for average process
2406 uInt insize = in.size() ; // number of input scantables
2407
2408 // TEST: do normal average in each table before IF grouping
2409 cout << "Do preliminary averaging" << endl ;
2410 vector< CountedPtr<Scantable> > tmpin( insize ) ;
2411 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2412 vector< CountedPtr<Scantable> > v( 1, in[itable] ) ;
2413 tmpin[itable] = average( v, mask, weight, avmode ) ;
2414 }
2415
2416 // warning
2417 cout << "Average spectra with different spectral resolution" << endl ;
2418 cout << endl ;
2419
2420 // temporarily set coordinfo
2421 vector<string> oldinfo( insize ) ;
2422 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2423 vector<string> coordinfo = in[itable]->getCoordInfo() ;
2424 oldinfo[itable] = coordinfo[0] ;
2425 coordinfo[0] = "Hz" ;
2426 tmpin[itable]->setCoordInfo( coordinfo ) ;
2427 }
2428
2429 // columns
2430 ScalarColumn<uInt> freqIDCol ;
2431 ScalarColumn<uInt> ifnoCol ;
2432 ScalarColumn<uInt> scannoCol ;
2433
2434
2435 // check IF frequency coverage
2436 // freqid: list of FREQ_ID, which is used, in each table
2437 // iffreq: list of minimum and maximum frequency for each FREQ_ID in
2438 // each table
2439 // freqid[insize][numIF]
2440 // freqid: [[id00, id01, ...],
2441 // [id10, id11, ...],
2442 // ...
2443 // [idn0, idn1, ...]]
2444 // iffreq[insize][numIF*2]
2445 // iffreq: [[min_id00, max_id00, min_id01, max_id01, ...],
2446 // [min_id10, max_id10, min_id11, max_id11, ...],
2447 // ...
2448 // [min_idn0, max_idn0, min_idn1, max_idn1, ...]]
2449 //cout << "Check IF settings in each table" << endl ;
2450 vector< vector<uInt> > freqid( insize );
2451 vector< vector<double> > iffreq( insize ) ;
2452 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2453 uInt rows = tmpin[itable]->nrow() ;
2454 uInt freqnrows = tmpin[itable]->frequencies().table().nrow() ;
2455 for ( uInt irow = 0 ; irow < rows ; irow++ ) {
2456 if ( freqid[itable].size() == freqnrows ) {
2457 break ;
2458 }
2459 else {
2460 freqIDCol.attach( tmpin[itable]->table(), "FREQ_ID" ) ;
2461 ifnoCol.attach( tmpin[itable]->table(), "IFNO" ) ;
2462 uInt id = freqIDCol( irow ) ;
2463 if ( freqid[itable].size() == 0 || count( freqid[itable].begin(), freqid[itable].end(), id ) == 0 ) {
2464 //cout << "itable = " << itable << ": IF " << id << " is included in the list" << endl ;
2465 vector<double> abcissa = tmpin[itable]->getAbcissa( irow ) ;
2466 freqid[itable].push_back( id ) ;
2467 iffreq[itable].push_back( abcissa[0] - 0.5 * ( abcissa[1] - abcissa[0] ) ) ;
2468 iffreq[itable].push_back( abcissa[abcissa.size()-1] + 0.5 * ( abcissa[1] - abcissa[0] ) ) ;
2469 }
2470 }
2471 }
2472 }
2473
2474 // debug
2475 //cout << "IF settings summary:" << endl ;
2476 //for ( uInt i = 0 ; i < freqid.size() ; i++ ) {
2477 //cout << " Table" << i << endl ;
2478 //for ( uInt j = 0 ; j < freqid[i].size() ; j++ ) {
2479 //cout << " id = " << freqid[i][j] << " (min,max) = (" << iffreq[i][2*j] << "," << iffreq[i][2*j+1] << ")" << endl ;
2480 //}
2481 //}
2482 //cout << endl ;
2483
2484 // IF grouping based on their frequency coverage
2485 // ifgrp: list of table index and FREQ_ID for all members in each IF group
2486 // ifgfreq: list of minimum and maximum frequency in each IF group
2487 // ifgrp[numgrp][nummember*2]
2488 // ifgrp: [[table00, freqrow00, table01, freqrow01, ...],
2489 // [table10, freqrow10, table11, freqrow11, ...],
2490 // ...
2491 // [tablen0, freqrown0, tablen1, freqrown1, ...]]
2492 // ifgfreq[numgrp*2]
2493 // ifgfreq: [min0_grp0, max0_grp0, min1_grp1, max1_grp1, ...]
2494 //cout << "IF grouping based on their frequency coverage" << endl ;
2495 vector< vector<uInt> > ifgrp ;
2496 vector<double> ifgfreq ;
2497
2498 // parameter for IF grouping
2499 // groupmode = OR retrieve all region
2500 // AND only retrieve overlaped region
2501 //string groupmode = "AND" ;
2502 string groupmode = "OR" ;
2503 uInt sizecr = 0 ;
2504 if ( groupmode == "AND" )
2505 sizecr = 2 ;
2506 else if ( groupmode == "OR" )
2507 sizecr = 0 ;
2508
2509 vector<double> sortedfreq ;
2510 for ( uInt i = 0 ; i < iffreq.size() ; i++ ) {
2511 for ( uInt j = 0 ; j < iffreq[i].size() ; j++ ) {
2512 if ( count( sortedfreq.begin(), sortedfreq.end(), iffreq[i][j] ) == 0 )
2513 sortedfreq.push_back( iffreq[i][j] ) ;
2514 }
2515 }
2516 sort( sortedfreq.begin(), sortedfreq.end() ) ;
2517 for ( vector<double>::iterator i = sortedfreq.begin() ; i != sortedfreq.end()-1 ; i++ ) {
2518 ifgfreq.push_back( *i ) ;
2519 ifgfreq.push_back( *(i+1) ) ;
2520 }
2521 ifgrp.resize( ifgfreq.size()/2 ) ;
2522 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2523 for ( uInt iif = 0 ; iif < freqid[itable].size() ; iif++ ) {
2524 double range0 = iffreq[itable][2*iif] ;
2525 double range1 = iffreq[itable][2*iif+1] ;
2526 for ( uInt j = 0 ; j < ifgrp.size() ; j++ ) {
2527 double fmin = max( range0, ifgfreq[2*j] ) ;
2528 double fmax = min( range1, ifgfreq[2*j+1] ) ;
2529 if ( fmin < fmax ) {
2530 ifgrp[j].push_back( itable ) ;
2531 ifgrp[j].push_back( freqid[itable][iif] ) ;
2532 }
2533 }
2534 }
2535 }
2536 vector< vector<uInt> >::iterator fiter = ifgrp.begin() ;
2537 vector<double>::iterator giter = ifgfreq.begin() ;
2538 while( fiter != ifgrp.end() ) {
2539 if ( fiter->size() <= sizecr ) {
2540 fiter = ifgrp.erase( fiter ) ;
2541 giter = ifgfreq.erase( giter ) ;
2542 giter = ifgfreq.erase( giter ) ;
2543 }
2544 else {
2545 fiter++ ;
2546 advance( giter, 2 ) ;
2547 }
2548 }
2549
2550 // Grouping continuous IF groups (without frequency gap)
2551 // freqgrp: list of IF group indexes in each frequency group
2552 // freqrange: list of minimum and maximum frequency in each frequency group
2553 // freqgrp[numgrp][nummember]
2554 // freqgrp: [[ifgrp00, ifgrp01, ifgrp02, ...],
2555 // [ifgrp10, ifgrp11, ifgrp12, ...],
2556 // ...
2557 // [ifgrpn0, ifgrpn1, ifgrpn2, ...]]
2558 // freqrange[numgrp*2]
2559 // freqrange: [min_grp0, max_grp0, min_grp1, max_grp1, ...]
2560 vector< vector<uInt> > freqgrp ;
2561 double freqrange = 0.0 ;
2562 uInt grpnum = 0 ;
2563 for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
2564 // Assumed that ifgfreq was sorted
2565 if ( grpnum != 0 && freqrange == ifgfreq[2*i] ) {
2566 freqgrp[grpnum-1].push_back( i ) ;
2567 }
2568 else {
2569 vector<uInt> grp0( 1, i ) ;
2570 freqgrp.push_back( grp0 ) ;
2571 grpnum++ ;
2572 }
2573 freqrange = ifgfreq[2*i+1] ;
2574 }
2575
2576
2577 // print IF groups
2578 cout << "IF Group summary: " << endl ;
2579 cout << " GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
2580 for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
2581 cout << " GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
2582 for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
2583 cout << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ;
2584 }
2585 cout << endl ;
2586 }
2587 cout << endl ;
2588
2589 // print frequency group
2590 cout << "Frequency Group summary: " << endl ;
2591 cout << " GROUP_ID [FREQ_MIN, FREQ_MAX]: IF_GROUP_ID" << endl ;
2592 for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
2593 cout << " GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*freqgrp[i][0]] << "," << ifgfreq[2*freqgrp[i][freqgrp[i].size()-1]+1] << "]: " ;
2594 for ( uInt j = 0 ; j < freqgrp[i].size() ; j++ ) {
2595 cout << freqgrp[i][j] << " " ;
2596 }
2597 cout << endl ;
2598 }
2599 cout << endl ;
2600
2601 // membership check
2602 // groups: list of IF group indexes whose frequency range overlaps with
2603 // that of each table and IF
2604 // groups[numtable][numIF][nummembership]
2605 // groups: [[[grp, grp,...], [grp, grp,...],...],
2606 // [[grp, grp,...], [grp, grp,...],...],
2607 // ...
2608 // [[grp, grp,...], [grp, grp,...],...]]
2609 vector< vector< vector<uInt> > > groups( insize ) ;
2610 for ( uInt i = 0 ; i < insize ; i++ ) {
2611 groups[i].resize( freqid[i].size() ) ;
2612 }
2613 for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
2614 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
2615 uInt tableid = ifgrp[igrp][2*imem] ;
2616 vector<uInt>::iterator iter = find( freqid[tableid].begin(), freqid[tableid].end(), ifgrp[igrp][2*imem+1] ) ;
2617 if ( iter != freqid[tableid].end() ) {
2618 uInt rowid = distance( freqid[tableid].begin(), iter ) ;
2619 groups[tableid][rowid].push_back( igrp ) ;
2620 }
2621 }
2622 }
2623
2624 // print membership
2625 //for ( uInt i = 0 ; i < insize ; i++ ) {
2626 //cout << "Table " << i << endl ;
2627 //for ( uInt j = 0 ; j < groups[i].size() ; j++ ) {
2628 //cout << " FREQ_ID " << setw( 2 ) << freqid[i][j] << ": " ;
2629 //for ( uInt k = 0 ; k < groups[i][j].size() ; k++ ) {
2630 //cout << setw( 2 ) << groups[i][j][k] << " " ;
2631 //}
2632 //cout << endl ;
2633 //}
2634 //}
2635
2636 // set back coordinfo
2637 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2638 vector<string> coordinfo = tmpin[itable]->getCoordInfo() ;
2639 coordinfo[0] = oldinfo[itable] ;
2640 tmpin[itable]->setCoordInfo( coordinfo ) ;
2641 }
2642
2643 // Create additional table if needed
2644 bool oldInsitu = insitu_ ;
2645 setInsitu( false ) ;
2646 vector< vector<uInt> > addrow( insize ) ;
2647 vector<uInt> addtable( insize, 0 ) ;
2648 vector<uInt> newtableids( insize ) ;
2649 vector<uInt> newifids( insize, 0 ) ;
2650 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2651 //cout << "Table " << setw(2) << itable << ": " ;
2652 for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
2653 addrow[itable].push_back( groups[itable][ifrow].size()-1 ) ;
2654 //cout << addrow[itable][ifrow] << " " ;
2655 }
2656 addtable[itable] = *max_element( addrow[itable].begin(), addrow[itable].end() ) ;
2657 //cout << "(" << addtable[itable] << ")" << endl ;
2658 }
2659 newin.resize( insize ) ;
2660 copy( tmpin.begin(), tmpin.end(), newin.begin() ) ;
2661 for ( uInt i = 0 ; i < insize ; i++ ) {
2662 newtableids[i] = i ;
2663 }
2664 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2665 for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
2666 CountedPtr<Scantable> add = getScantable( newin[itable], false ) ;
2667 vector<int> freqidlist ;
2668 for ( uInt i = 0 ; i < groups[itable].size() ; i++ ) {
2669 if ( groups[itable][i].size() > iadd + 1 ) {
2670 freqidlist.push_back( freqid[itable][i] ) ;
2671 }
2672 }
2673 stringstream taqlstream ;
2674 taqlstream << "SELECT FROM $1 WHERE FREQ_ID IN [" ;
2675 for ( uInt i = 0 ; i < freqidlist.size() ; i++ ) {
2676 taqlstream << i ;
2677 if ( i < freqidlist.size() - 1 )
2678 taqlstream << "," ;
2679 else
2680 taqlstream << "]" ;
2681 }
2682 string taql = taqlstream.str() ;
2683 //cout << "taql = " << taql << endl ;
2684 STSelector selector = STSelector() ;
2685 selector.setTaQL( taql ) ;
2686 add->setSelection( selector ) ;
2687 newin.push_back( add ) ;
2688 newtableids.push_back( itable ) ;
2689 newifids.push_back( iadd + 1 ) ;
2690 }
2691 }
2692
2693 // udpate ifgrp
2694 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2695 for ( uInt iadd = 0 ; iadd < addtable[itable] ; iadd++ ) {
2696 for ( uInt ifrow = 0 ; ifrow < groups[itable].size() ; ifrow++ ) {
2697 if ( groups[itable][ifrow].size() > iadd + 1 ) {
2698 uInt igrp = groups[itable][ifrow][iadd+1] ;
2699 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
2700 if ( ifgrp[igrp][2*imem] == newtableids[iadd+insize] && ifgrp[igrp][2*imem+1] == freqid[newtableids[iadd+insize]][ifrow] ) {
2701 ifgrp[igrp][2*imem] = insize + iadd ;
2702 }
2703 }
2704 }
2705 }
2706 }
2707 }
2708
2709 // print IF groups again for debug
2710 //cout << "IF Group summary: " << endl ;
2711 //cout << " GROUP_ID [FREQ_MIN, FREQ_MAX]: (TABLE_ID, FREQ_ID)" << endl ;
2712 //for ( uInt i = 0 ; i < ifgrp.size() ; i++ ) {
2713 //cout << " GROUP " << setw( 2 ) << i << " [" << ifgfreq[2*i] << "," << ifgfreq[2*i+1] << "]: " ;
2714 //for ( uInt j = 0 ; j < ifgrp[i].size()/2 ; j++ ) {
2715 //cout << "(" << ifgrp[i][2*j] << "," << ifgrp[i][2*j+1] << ") " ;
2716 //}
2717 //cout << endl ;
2718 //}
2719 //cout << endl ;
2720
2721 // reset SCANNO and IFNO/FREQ_ID: IF is reset by the result of sortation
2722 cout << "All scan number is set to 0" << endl ;
2723 //cout << "All IF number is set to IF group index" << endl ;
2724 cout << endl ;
2725 insize = newin.size() ;
2726 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2727 uInt rows = newin[itable]->nrow() ;
2728 Table &tmpt = newin[itable]->table() ;
2729 freqIDCol.attach( tmpt, "FREQ_ID" ) ;
2730 scannoCol.attach( tmpt, "SCANNO" ) ;
2731 ifnoCol.attach( tmpt, "IFNO" ) ;
2732 for ( uInt irow=0 ; irow < rows ; irow++ ) {
2733 scannoCol.put( irow, 0 ) ;
2734 uInt freqID = freqIDCol( irow ) ;
2735 vector<uInt>::iterator iter = find( freqid[newtableids[itable]].begin(), freqid[newtableids[itable]].end(), freqID ) ;
2736 if ( iter != freqid[newtableids[itable]].end() ) {
2737 uInt index = distance( freqid[newtableids[itable]].begin(), iter ) ;
2738 ifnoCol.put( irow, groups[newtableids[itable]][index][newifids[itable]] ) ;
2739 }
2740 else {
2741 throw(AipsError("IF grouping was wrong in additional tables.")) ;
2742 }
2743 }
2744 }
2745 oldinfo.resize( insize ) ;
2746 setInsitu( oldInsitu ) ;
2747
2748 // temporarily set coordinfo
2749 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2750 vector<string> coordinfo = newin[itable]->getCoordInfo() ;
2751 oldinfo[itable] = coordinfo[0] ;
2752 coordinfo[0] = "Hz" ;
2753 newin[itable]->setCoordInfo( coordinfo ) ;
2754 }
2755
2756 // save column values in the vector
2757 vector< vector<uInt> > freqTableIdVec( insize ) ;
2758 vector< vector<uInt> > freqIdVec( insize ) ;
2759 vector< vector<uInt> > ifNoVec( insize ) ;
2760 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2761 ScalarColumn<uInt> freqIDs ;
2762 freqIDs.attach( newin[itable]->frequencies().table(), "ID" ) ;
2763 ifnoCol.attach( newin[itable]->table(), "IFNO" ) ;
2764 freqIDCol.attach( newin[itable]->table(), "FREQ_ID" ) ;
2765 for ( uInt irow = 0 ; irow < newin[itable]->frequencies().table().nrow() ; irow++ ) {
2766 freqTableIdVec[itable].push_back( freqIDs( irow ) ) ;
2767 }
2768 for ( uInt irow = 0 ; irow < newin[itable]->table().nrow() ; irow++ ) {
2769 freqIdVec[itable].push_back( freqIDCol( irow ) ) ;
2770 ifNoVec[itable].push_back( ifnoCol( irow ) ) ;
2771 }
2772 }
2773
2774 // reset spectra and flagtra: pick up common part of frequency coverage
2775 //cout << "Pick common frequency range and align resolution" << endl ;
2776 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2777 uInt rows = newin[itable]->nrow() ;
2778 int nminchan = -1 ;
2779 int nmaxchan = -1 ;
2780 vector<uInt> freqIdUpdate ;
2781 for ( uInt irow = 0 ; irow < rows ; irow++ ) {
2782 uInt ifno = ifNoVec[itable][irow] ; // IFNO is reset by group index
2783 double minfreq = ifgfreq[2*ifno] ;
2784 double maxfreq = ifgfreq[2*ifno+1] ;
2785 //cout << "frequency range: [" << minfreq << "," << maxfreq << "]" << endl ;
2786 vector<double> abcissa = newin[itable]->getAbcissa( irow ) ;
2787 int nchan = abcissa.size() ;
2788 double resol = abcissa[1] - abcissa[0] ;
2789 //cout << "abcissa range : [" << abcissa[0] << "," << abcissa[nchan-1] << "]" << endl ;
2790 if ( minfreq <= abcissa[0] )
2791 nminchan = 0 ;
2792 else {
2793 //double cfreq = ( minfreq - abcissa[0] ) / resol ;
2794 double cfreq = ( minfreq - abcissa[0] + 0.5 * resol ) / resol ;
2795 nminchan = int(cfreq) + ( ( cfreq - int(cfreq) <= 0.5 ) ? 0 : 1 ) ;
2796 }
2797 if ( maxfreq >= abcissa[abcissa.size()-1] )
2798 nmaxchan = abcissa.size() - 1 ;
2799 else {
2800 //double cfreq = ( abcissa[abcissa.size()-1] - maxfreq ) / resol ;
2801 double cfreq = ( abcissa[abcissa.size()-1] - maxfreq + 0.5 * resol ) / resol ;
2802 nmaxchan = abcissa.size() - 1 - int(cfreq) - ( ( cfreq - int(cfreq) >= 0.5 ) ? 1 : 0 ) ;
2803 }
2804 //cout << "channel range (" << irow << "): [" << nminchan << "," << nmaxchan << "]" << endl ;
2805 if ( nmaxchan > nminchan ) {
2806 newin[itable]->reshapeSpectrum( nminchan, nmaxchan, irow ) ;
2807 int newchan = nmaxchan - nminchan + 1 ;
2808 if ( count( freqIdUpdate.begin(), freqIdUpdate.end(), freqIdVec[itable][irow] ) == 0 && newchan < nchan )
2809 freqIdUpdate.push_back( freqIdVec[itable][irow] ) ;
2810 }
2811 else {
2812 throw(AipsError("Failed to pick up common part of frequency range.")) ;
2813 }
2814 }
2815 for ( uInt i = 0 ; i < freqIdUpdate.size() ; i++ ) {
2816 uInt freqId = freqIdUpdate[i] ;
2817 Double refpix ;
2818 Double refval ;
2819 Double increment ;
2820
2821 // update row
2822 newin[itable]->frequencies().getEntry( refpix, refval, increment, freqId ) ;
2823 refval = refval - ( refpix - nminchan ) * increment ;
2824 refpix = 0 ;
2825 newin[itable]->frequencies().setEntry( refpix, refval, increment, freqId ) ;
2826 }
2827 }
2828
2829
2830 // reset spectra and flagtra: align spectral resolution
2831 //cout << "Align spectral resolution" << endl ;
2832 // gmaxdnu: the coarsest frequency resolution in the frequency group
2833 // gmemid: member index that have a resolution equal to gmaxdnu
2834 // gmaxdnu[numfreqgrp]
2835 // gmaxdnu: [dnu0, dnu1, ...]
2836 // gmemid[numfreqgrp]
2837 // gmemid: [id0, id1, ...]
2838 vector<double> gmaxdnu( freqgrp.size(), 0.0 ) ;
2839 vector<uInt> gmemid( freqgrp.size(), 0 ) ;
2840 for ( uInt igrp = 0 ; igrp < ifgrp.size() ; igrp++ ) {
2841 double maxdnu = 0.0 ; // maximum (coarsest) frequency resolution
2842 int minchan = INT_MAX ; // minimum channel number
2843 Double refpixref = -1 ; // reference of 'reference pixel'
2844 Double refvalref = -1 ; // reference of 'reference frequency'
2845 Double refinc = -1 ; // reference frequency resolution
2846 uInt refreqid ;
2847 uInt reftable = INT_MAX;
2848 // process only if group member > 1
2849 if ( ifgrp[igrp].size() > 2 ) {
2850 // find minchan and maxdnu in each group
2851 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
2852 uInt tableid = ifgrp[igrp][2*imem] ;
2853 uInt rowid = ifgrp[igrp][2*imem+1] ;
2854 vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
2855 if ( iter != freqIdVec[tableid].end() ) {
2856 uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
2857 vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
2858 int nchan = abcissa.size() ;
2859 double dnu = abcissa[1] - abcissa[0] ;
2860 //cout << "GROUP " << igrp << " (" << tableid << "," << rowid << "): nchan = " << nchan << " (minchan = " << minchan << ")" << endl ;
2861 if ( nchan < minchan ) {
2862 minchan = nchan ;
2863 maxdnu = dnu ;
2864 newin[tableid]->frequencies().getEntry( refpixref, refvalref, refinc, rowid ) ;
2865 refreqid = rowid ;
2866 reftable = tableid ;
2867 }
2868 }
2869 }
2870 // regrid spectra in each group
2871 cout << "GROUP " << igrp << endl ;
2872 cout << " Channel number is adjusted to " << minchan << endl ;
2873 cout << " Corresponding frequency resolution is " << maxdnu << "Hz" << endl ;
2874 for ( uInt imem = 0 ; imem < ifgrp[igrp].size()/2 ; imem++ ) {
2875 uInt tableid = ifgrp[igrp][2*imem] ;
2876 uInt rowid = ifgrp[igrp][2*imem+1] ;
2877 freqIDCol.attach( newin[tableid]->table(), "FREQ_ID" ) ;
2878 //cout << "tableid = " << tableid << " rowid = " << rowid << ": " << endl ;
2879 //cout << " regridChannel applied to " ;
2880 if ( tableid != reftable )
2881 refreqid = newin[tableid]->frequencies().addEntry( refpixref, refvalref, refinc ) ;
2882 for ( uInt irow = 0 ; irow < newin[tableid]->table().nrow() ; irow++ ) {
2883 uInt tfreqid = freqIdVec[tableid][irow] ;
2884 if ( tfreqid == rowid ) {
2885 //cout << irow << " " ;
2886 newin[tableid]->regridChannel( minchan, maxdnu, irow ) ;
2887 freqIDCol.put( irow, refreqid ) ;
2888 freqIdVec[tableid][irow] = refreqid ;
2889 }
2890 }
2891 //cout << endl ;
2892 }
2893 }
2894 else {
2895 uInt tableid = ifgrp[igrp][0] ;
2896 uInt rowid = ifgrp[igrp][1] ;
2897 vector<uInt>::iterator iter = find( freqIdVec[tableid].begin(), freqIdVec[tableid].end(), rowid ) ;
2898 if ( iter != freqIdVec[tableid].end() ) {
2899 uInt index = distance( freqIdVec[tableid].begin(), iter ) ;
2900 vector<double> abcissa = newin[tableid]->getAbcissa( index ) ;
2901 minchan = abcissa.size() ;
2902 maxdnu = abcissa[1] - abcissa[0] ;
2903 }
2904 }
2905 for ( uInt i = 0 ; i < freqgrp.size() ; i++ ) {
2906 if ( count( freqgrp[i].begin(), freqgrp[i].end(), igrp ) > 0 ) {
2907 if ( maxdnu > gmaxdnu[i] ) {
2908 gmaxdnu[i] = maxdnu ;
2909 gmemid[i] = igrp ;
2910 }
2911 break ;
2912 }
2913 }
2914 }
2915 cout << endl ;
2916
2917 // set back coordinfo
2918 for ( uInt itable = 0 ; itable < insize ; itable++ ) {
2919 vector<string> coordinfo = newin[itable]->getCoordInfo() ;
2920 coordinfo[0] = oldinfo[itable] ;
2921 newin[itable]->setCoordInfo( coordinfo ) ;
2922 }
2923
2924 // accumulate all rows into the first table
2925 // NOTE: assumed in.size() = 1
2926 vector< CountedPtr<Scantable> > tmp( 1 ) ;
2927 if ( newin.size() == 1 )
2928 tmp[0] = newin[0] ;
2929 else
2930 tmp[0] = merge( newin ) ;
2931
2932 //return tmp[0] ;
2933
2934 // average
2935 CountedPtr<Scantable> tmpout = average( tmp, mask, weight, avmode ) ;
2936
2937 //return tmpout ;
2938
2939 // combine frequency group
2940 cout << "Combine spectra based on frequency grouping" << endl ;
2941 cout << "IFNO is renumbered as frequency group ID (see above)" << endl ;
2942 vector<string> coordinfo = tmpout->getCoordInfo() ;
2943 oldinfo[0] = coordinfo[0] ;
2944 coordinfo[0] = "Hz" ;
2945 tmpout->setCoordInfo( coordinfo ) ;
2946 // create proformas of output table
2947 stringstream taqlstream ;
2948 taqlstream << "SELECT FROM $1 WHERE IFNO IN [" ;
2949 for ( uInt i = 0 ; i < gmemid.size() ; i++ ) {
2950 taqlstream << gmemid[i] ;
2951 if ( i < gmemid.size() - 1 )
2952 taqlstream << "," ;
2953 else
2954 taqlstream << "]" ;
2955 }
2956 string taql = taqlstream.str() ;
2957 //cout << "taql = " << taql << endl ;
2958 STSelector selector = STSelector() ;
2959 selector.setTaQL( taql ) ;
2960 oldInsitu = insitu_ ;
2961 setInsitu( false ) ;
2962 out = getScantable( tmpout, false ) ;
2963 setInsitu( oldInsitu ) ;
2964 out->setSelection( selector ) ;
2965 // regrid rows
2966 ifnoCol.attach( tmpout->table(), "IFNO" ) ;
2967 for ( uInt irow = 0 ; irow < tmpout->table().nrow() ; irow++ ) {
2968 uInt ifno = ifnoCol( irow ) ;
2969 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
2970 if ( count( freqgrp[igrp].begin(), freqgrp[igrp].end(), ifno ) > 0 ) {
2971 vector<double> abcissa = tmpout->getAbcissa( irow ) ;
2972 double bw = ( abcissa[1] - abcissa[0] ) * abcissa.size() ;
2973 int nchan = (int)( bw / gmaxdnu[igrp] ) ;
2974 tmpout->regridChannel( nchan, gmaxdnu[igrp], irow ) ;
2975 break ;
2976 }
2977 }
2978 }
2979 // combine spectra
2980 ArrayColumn<Float> specColOut ;
2981 specColOut.attach( out->table(), "SPECTRA" ) ;
2982 ArrayColumn<uChar> flagColOut ;
2983 flagColOut.attach( out->table(), "FLAGTRA" ) ;
2984 ScalarColumn<uInt> ifnoColOut ;
2985 ifnoColOut.attach( out->table(), "IFNO" ) ;
2986 ScalarColumn<uInt> polnoColOut ;
2987 polnoColOut.attach( out->table(), "POLNO" ) ;
2988 ScalarColumn<uInt> freqidColOut ;
2989 freqidColOut.attach( out->table(), "FREQ_ID" ) ;
2990 MDirection::ScalarColumn dirColOut ;
2991 dirColOut.attach( out->table(), "DIRECTION" ) ;
2992 Table &tab = tmpout->table() ;
2993 Block<String> cols(1);
2994 cols[0] = String("POLNO") ;
2995 TableIterator iter( tab, cols ) ;
2996 bool done = false ;
2997 vector< vector<uInt> > sizes( freqgrp.size() ) ;
2998 while( !iter.pastEnd() ) {
2999 vector< vector<Float> > specout( freqgrp.size() ) ;
3000 vector< vector<uChar> > flagout( freqgrp.size() ) ;
3001 ArrayColumn<Float> specCols ;
3002 specCols.attach( iter.table(), "SPECTRA" ) ;
3003 ArrayColumn<uChar> flagCols ;
3004 flagCols.attach( iter.table(), "FLAGTRA" ) ;
3005 ifnoCol.attach( iter.table(), "IFNO" ) ;
3006 ScalarColumn<uInt> polnos ;
3007 polnos.attach( iter.table(), "POLNO" ) ;
3008 MDirection::ScalarColumn dircol ;
3009 dircol.attach( iter.table(), "DIRECTION" ) ;
3010 uInt polno = polnos( 0 ) ;
3011 //cout << "POLNO iteration: " << polno << endl ;
3012// for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3013// sizes[igrp].resize( freqgrp[igrp].size() ) ;
3014// for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
3015// for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
3016// uInt ifno = ifnoCol( irow ) ;
3017// if ( ifno == freqgrp[igrp][imem] ) {
3018// Vector<Float> spec = specCols( irow ) ;
3019// Vector<uChar> flag = flagCols( irow ) ;
3020// vector<Float> svec ;
3021// spec.tovector( svec ) ;
3022// vector<uChar> fvec ;
3023// flag.tovector( fvec ) ;
3024// //cout << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << endl ;
3025// specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ;
3026// flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ;
3027// //cout << "specout[" << igrp << "].size() = " << specout[igrp].size() << endl ;
3028// sizes[igrp][imem] = spec.nelements() ;
3029// }
3030// }
3031// }
3032// for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
3033// uInt ifout = ifnoColOut( irow ) ;
3034// uInt polout = polnoColOut( irow ) ;
3035// if ( ifout == gmemid[igrp] && polout == polno ) {
3036// // set SPECTRA and FRAGTRA
3037// Vector<Float> newspec( specout[igrp] ) ;
3038// Vector<uChar> newflag( flagout[igrp] ) ;
3039// specColOut.put( irow, newspec ) ;
3040// flagColOut.put( irow, newflag ) ;
3041// // IFNO renumbering
3042// ifnoColOut.put( irow, igrp ) ;
3043// }
3044// }
3045// }
3046 // get a list of number of channels for each frequency group member
3047 if ( !done ) {
3048 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3049 sizes[igrp].resize( freqgrp[igrp].size() ) ;
3050 for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
3051 for ( uInt irow = 0 ; irow < iter.table().nrow() ; irow++ ) {
3052 uInt ifno = ifnoCol( irow ) ;
3053 if ( ifno == freqgrp[igrp][imem] ) {
3054 Vector<Float> spec = specCols( irow ) ;
3055 sizes[igrp][imem] = spec.nelements() ;
3056 break ;
3057 }
3058 }
3059 }
3060 }
3061 done = true ;
3062 }
3063 // combine spectra
3064 for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
3065 uInt polout = polnoColOut( irow ) ;
3066 if ( polout == polno ) {
3067 uInt ifout = ifnoColOut( irow ) ;
3068 Vector<Double> direction = dirColOut(irow).getAngle(Unit(String("rad"))).getValue() ;
3069 uInt igrp ;
3070 for ( uInt jgrp = 0 ; jgrp < freqgrp.size() ; jgrp++ ) {
3071 if ( ifout == gmemid[jgrp] ) {
3072 igrp = jgrp ;
3073 break ;
3074 }
3075 }
3076 for ( uInt imem = 0 ; imem < freqgrp[igrp].size() ; imem++ ) {
3077 for ( uInt jrow = 0 ; jrow < iter.table().nrow() ; jrow++ ) {
3078 uInt ifno = ifnoCol( jrow ) ;
3079 Vector<Double> tdir = dircol(jrow).getAngle(Unit(String("rad"))).getValue() ;
3080 //if ( ifno == freqgrp[igrp][imem] && allTrue( tdir == direction ) ) {
3081 Double dx = tdir[0] - direction[0] ;
3082 Double dy = tdir[1] - direction[1] ;
3083 Double dd = sqrt( dx * dx + dy * dy ) ;
3084 //if ( ifno == freqgrp[igrp][imem] && allNearAbs( tdir, direction, tol ) ) {
3085 if ( ifno == freqgrp[igrp][imem] && dd <= tol ) {
3086 Vector<Float> spec = specCols( jrow ) ;
3087 Vector<uChar> flag = flagCols( jrow ) ;
3088 vector<Float> svec ;
3089 spec.tovector( svec ) ;
3090 vector<uChar> fvec ;
3091 flag.tovector( fvec ) ;
3092 //cout << "spec.size() = " << svec.size() << " fvec.size() = " << fvec.size() << endl ;
3093 specout[igrp].insert( specout[igrp].end(), svec.begin(), svec.end() ) ;
3094 flagout[igrp].insert( flagout[igrp].end(), fvec.begin(), fvec.end() ) ;
3095 //cout << "specout[" << igrp << "].size() = " << specout[igrp].size() << endl ;
3096 }
3097 }
3098 }
3099 // set SPECTRA and FRAGTRA
3100 Vector<Float> newspec( specout[igrp] ) ;
3101 Vector<uChar> newflag( flagout[igrp] ) ;
3102 specColOut.put( irow, newspec ) ;
3103 flagColOut.put( irow, newflag ) ;
3104 // IFNO renumbering
3105 ifnoColOut.put( irow, igrp ) ;
3106 }
3107 }
3108 iter++ ;
3109 }
3110 // update FREQUENCIES subtable
3111 vector<bool> updated( freqgrp.size(), false ) ;
3112 for ( uInt igrp = 0 ; igrp < freqgrp.size() ; igrp++ ) {
3113 uInt index = 0 ;
3114 uInt pixShift = 0 ;
3115 while ( freqgrp[igrp][index] != gmemid[igrp] ) {
3116 pixShift += sizes[igrp][index++] ;
3117 }
3118 for ( uInt irow = 0 ; irow < out->table().nrow() ; irow++ ) {
3119 if ( ifnoColOut( irow ) == gmemid[igrp] && !updated[igrp] ) {
3120 uInt freqidOut = freqidColOut( irow ) ;
3121 //cout << "freqgrp " << igrp << " freqidOut = " << freqidOut << endl ;
3122 double refpix ;
3123 double refval ;
3124 double increm ;
3125 out->frequencies().getEntry( refpix, refval, increm, freqidOut ) ;
3126 refpix += pixShift ;
3127 out->frequencies().setEntry( refpix, refval, increm, freqidOut ) ;
3128 updated[igrp] = true ;
3129 }
3130 }
3131 }
3132
3133 //out = tmpout ;
3134
3135 coordinfo = tmpout->getCoordInfo() ;
3136 coordinfo[0] = oldinfo[0] ;
3137 tmpout->setCoordInfo( coordinfo ) ;
3138 }
3139 else {
3140 // simple average
3141 out = average( in, mask, weight, avmode ) ;
3142 }
3143
3144 return out ;
3145}
Note: See TracBrowser for help on using the repository browser.