source: branches/hpc33/src/STMath.cpp@ 2446

Last change on this file since 2446 was 2443, checked in by Takeshi Nakazato, 13 years ago

New Development: No

JIRA Issue: No

Ready for Test: No

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs:

Put in Release Notes: No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Test tuning ALMA calibration code.


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