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

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

New Development: No

JIRA Issue: Yes CAS-1799

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: srcType is added to PKSrecord class

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): casa/atnf

Description: Describe your changes here...

Calibration scheme is updated to use SRCTYPE as an identifier of scan intent
(ON, OFF, ...). Currently, updated calibration scheme is implemented only for
ALMA data. Calibration of other data is still use SRCNAME as the identifier.

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