source: branches/casa-prerelease/pre-asap/src/STMath.cpp@ 2307

Last change on this file since 2307 was 2289, checked in by ShinnosukeKawakami, 13 years ago

merged parallel branch to trunk

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