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

Last change on this file since 2362 was 2345, checked in by WataruKawasaki, 13 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): SD

Description: bugfix.


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