source: branches/alma/src/STMath.cpp@ 1682

Last change on this file since 1682 was 1682, checked in by Takeshi Nakazato, 15 years ago

New Development: No

JIRA Issue: Yes CAS-1823

Ready to Release: Yes/No

Interface Changes: Yes/No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Bug fix.
Supported operation with 1D list with length of 1, and with 2D list with
shape of (1,1) or (1,nchan).


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