source: branches/parallel/src/STMath.cpp@ 2047

Last change on this file since 2047 was 2014, checked in by ShinnosukeKawakami, 14 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): STMath.cpp

Description: calsig and calref were to be parallelized


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