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

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

New Development: No

JIRA Issue: N

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...


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