source: branches/parallel/src/STMath.cpp@ 2291

Last change on this file since 2291 was 2288, checked in by KohjiNakamura, 13 years ago

dototalpower fixed

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