source: trunk/src/STMath.cpp@ 2285

Last change on this file since 2285 was 2278, checked in by Malte Marquarding, 13 years ago

Interestingly this fix nver madi it into trunk. Need to iterate over subtable

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