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

Last change on this file since 2538 was 2534, checked in by Takeshi Nakazato, 12 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

calibrateALMA will do nothing when input ON scantable is empty row.


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