source: trunk/src/STMath.cpp@ 2526

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

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

bug fixes in STMath::new_average

  • proper handling of scan average
  • proper handling of multi-beam data


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