source: trunk/src/STMath.cpp@ 2203

Last change on this file since 2203 was 2186, checked in by WataruKawasaki, 13 years ago

New Development: Yes

JIRA Issue: Yes CAS-3149

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: scantable.*sinusoid_baseline() params

Test Programs:

Put in Release Notes: Yes

Module(s):

Description: (1) Implemented an automated sinusoidal fitting functionality

(2) FFT available with scantable.fft()
(3) fixed a bug of parsing 'edge' param used by linefinder.
(4) a function to show progress status for row-based iterations.


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