source: branches/polybatch/src/STMath.cpp@ 2383

Last change on this file since 2383 was 1907, checked in by WataruKawasaki, 14 years ago

New Development: No

JIRA Issue: Yes CAS-1937,CAS-2373

Ready to Release: Yes

Interface Changes: Yes

What Interface Changed: A new parameter 'batch' was added to

sd.scantable.poly_baseline(), while
'uselin' was removed.

Test Programs:

Put in Release Notes: Yes

Module(s): sdbaseline

Description: A faster version of sd.scantable.poly_baseline().


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